DLDNN: Deterministic Lateral Displacement Design Automation by Neural Networks

A DLD array of circular posts was built in the CFD solver to generate a comprehensive data-set for neural network training. To do so, the incompressible Navier-Stokes and continuity equations were solved since the Reynolds number lies within the laminar flow regime. Figure 1a illustrates a general schematic of the DLD design and the role of critical diameter in the separation of green ($D_g>D_c$) and red ($D_r<D_c$) particles. The dashed segment demonstrates the excerpt used for solving the periodic flow model which holds all the geometrical parameters that influence the critical diameter. These geometrical parameters consist of R, N, $G_x$, and $G_y$. Where, R is pillars radius, N is the pillars arrays period number ($N=1/\alpha$) defining their periodic repetition, and $G=G_x=G_y$ are the identical horizontal and vertical gaps between pillars. Then, in order to feed the velocity fields to the proposed CNN, the decomposed fields in the Cartesian coordinate system were mapped to squared fields of 128$\times$128 structured data as drawn in Figure 1b. To achieve this mapping, the domain was normalized by $L=2R+G$, followed by a shear mapping in the y direction as,

Field-generative neural networks provide a great alternative for CFD solvers with lower computational costs and the capacity of transfer learning for similar applications. Here, a convolutional network was used for reconstructing velocity fields as illustrated in Figure 1c. The neural network architecture used in this work, inspired by Dosovitskiy et al. [dosovitskiy2016learning], generates data as a great asset in DLD cases since the existing methods can be time-consuming and resource intensive. To capture a wide range of operational conditions, three inputs were defined namely, Re, N, and f where,

It is important to take the ratio of G and R instead of their values because the latter does not necessarily represent a unique shape. The inputs were fed into a separate set of dense layers for higher dimensional representation and after concatenation, they went through two sequential dense layers. Next, they were connected to four upsampling layers combined with convolutional. It should be noted that since both the x and y components of the velocity field are required, the mentioned sub-structure is used in parallel in order to reconstruct both fields.

Particle tracing and critical diameter extraction are crucial parts of the performance evaluation of DLD devices. Consequently, particles’ behavior in the flow field and their interaction with obstacles were simulated based on the following models. First, the massless particle tracing scheme was utilized, because the particles and the flow have negligible inertia. Then, particles positions were approximated from velocity fields with Runge-Kutta method. Finally, for simulating the particle-pillar interaction, the wall distance function was numerically calculated to measure the minimum spatial distance from pillars. In case of contact (when the distance between particle and pillar is equal to the particle radius), the reflect function modifies the particle velocity ($V_p$) as,

where the $\overrightarrow{n}$ and $\overrightarrow{t}$ are the unit vectors in the normal and tangential directions of particle and pillar contact point, respectively (see Appendix A). Once the particle’s accurate pathlines are obtained, the critical diameter was extracted via the Newtonian root finding method. In this approach, first, two particles with the highest and lowest diameters possible are simulated. The particle with the bumped mode receives a positive value and the particle with the zigzag mode is assigned a negative value. Thus, the Newtonian method can narrow the diameter interval to detect the root which is the critical diameter (see Appendix B).

The design automation algorithms convert the user-defined specification to geometry properties and flow rate to readily meet the desired criteria. To boost the function evaluation speed, one can utilize neural networks as the surrogate model instead of the whole physics/experimentation. This alternative facilitates the design process of a DLD device to reach the best possible configuration, especially when combined with an evolutionary algorithm.

As presented in Figure 2, design automation inputs are the diameters of particles to be separated ($D_1$ and $D_2$), desired constraints to be applied on f, N, Re ($C_f$, $C_N$, and $C_{Re}$), and finally the trade-off between flexibility and stability of the design ($\varphi$). Flexibility indicates the ability of the proposed design to cover a wide range of critical diameters while stability specifies its immunity to change due to minor flow rate fluctuations (see Appendix C). These two terms were included to consider the possible device behavior after fabrication since it can be altered by adjusting the Reynolds number.

After defining the inputs, a multi-objective optimization algorithm coupled with a pre-trained FCNN tunes the f, N, Re, and G to reach the optimized design. In this design automation tool, an NSGA3 with five reference directions and population size of 260 was chosen and implemented in Pymoo library [pymoo] in Python 3. The pre-trained FCNN was used for direct critical diameter prediction (direct neural network) since extracting critical diameter from CNN outputs slows down the optimization process. Finally, the tuned parameters in addition to bandwidth (BW) are extracted and presented in the outputs. BW is the difference between the maximum and the minimum critical diameter for that design serving as an index of flexibility. Note that the parameter G is used in optimization to nondimensionalize the critical diameter matching the developed neural network’s $D_c$.

Furthermore, the pre-trained CNN was attached to the end of the automation process to provide the ability to simulate the particle trajectory in the optimum design for a given number of periods. This feature helps visualize the particle behavior at the full length of the device.

The CNN was developed for fast velocity field generation which accelerates particle tracing. In this section, the network architecture, performance, and utility for data augmentation were analyzed. For this purpose, a sufficient amount of data was selected to cover a wide range of DLD device configurations. In order to choose the domain limits, a combination of reported working conditions from the literature was used [tang2022geometric]. The range for parameters f, N, and Re are as follows. f was chosen from 0.25 to 0.75 with step size of 0.02 which resulted in 26 data points. N was chosen to be {3, 4, 5, 6, 7, 8, 9, 10} and the Reynolds number was selected to be {0.01, 0.1, 1, 2.5, 5, 7.5, 10, 15, 20, 25}. In total, 2288 data points were created and split 80%-20% for train-development sets. Furthermore, another data-set was created to test the network, consisting of data points that were not given to the network in the development process. This data-set was build from a combination of f between 0.25 and 0.75 with a step size of 0.06, N of {3, 4, 5, 6}, and Re of {0.05, 1.5, 6.5, 8.5, 12.5, 18.5}, which resulted in a data-set with 239 data points.

Next, a network architecture was designed to effectively understand the relationship between the data. As can be seen in Table 1, nine architectures were tested starting with a base architecture (CNN1). At first, the width of the dense layers was increased (CNN1-3) which caused no significant improvement despite increasing the network trainable parameters dramatically. Then, the width and depth of convolutional layers were examined (CNN4-9) where CNN9 exhibited the best performance and its detailed structure is presented in Table 2. Moreover, other hyper-parameters were set to constant as follows, batch size of 64, a learning rate of $2e-3$ for the first 100 epochs, and $2e-4$ for the other 100 epochs. The chosen network was also gone through another 100 epochs with a learning rate of $2e-5$ and reached the validation loss of $5.1e-6$.

Finally, Figure 3 depicts the performance evaluation of the chosen CNN structure to predict the critical diameter. The accuracy of the network in development and test set is illustrated in Figure 3a, where its mean squared errors ($MSE$) were $2.7e-3$ and $2.6e-3$, respectively. Furthermore, The network predictions’ $MSE$ and standard deviation ($SD$) are also demonstrated based on Re, N, and f on the development set to detect any trends in the errors (Figure 3b). In sum, $MSE$ and $SD$ are less than 3% and 5%, respectively, which is an indication of the satisfactory performance of the network. However, the error trends convey that it will be intensified with the increase in N and f. Yet, there is no strong correlation between the error and Re. Additionally, Figure 3c, depicts the velocity fields with the most critical diameter prediction error from the test set ($f=0.38$, $N=8$, and $Re=15$). By comparing the original field from the numerical study (Ground Truth) and the neural network-generated field (Prediction), it is clear that these fields match well. Moreover, the discrepancies between the x and y components of velocity fields (Difference) quantitatively prove the latter observation. The results reveal that the nominated CNN has a high accuracy ($error<1\%$) in the field generation leaving most of the critical diameter prediction errors in the process of critical diameter computation.

Then, the trained CNN was used to generate new data (data augmentation) to facilitate the training of the direct neural network. This network helped with the augmentation of 2288 to 10400 data points by enhancing the resolution of Re in the new data-set by decreasing the step size to 0.499.

Here, the process of finding the best architecture of the FCNN for direct prediction of critical diameter alongside its training procedure is presented. For this particular problem, a systematic search was carried out where networks with a maximum of 10 hidden layers (HL) and nodes (N) of up to 128 were tested. The results are given in Table 3 after training with 1000 epochs and a learning rate of $1e-4$. The results imply that the 8 HL and 128 N architecture display superior performance. Therefore, this structure was used in the optimization algorithm.

The convergence process and performance of the chosen network are depicted in Figure 4. It is shown that the network’s $MSE$ on the training set (training loss) and the development set (validation loss) decreased per episode, Figure 4a. Considering the proximity and low values of training and validation losses ($2.8e-5$ and $3.8e-5$, respectively), it can be concluded that there were no avoidable bias and over-fitting problems. In addition, the network predictions of critical diameter and the actual critical diameter were compared on the entire original and generated data-sets. The network prediction has satisfying performance in both data-sets with the $MSE$ of $4.6e-3$ on the 2288 data-set and $4.0e-3$ on the 10400 data-set, Figure 4b.

After ensuring the accuracy of the direct neural network, the next step was examining the robustness and versatility of the design automation platform consisting of the network and evolutionary algorithm. In order to do this, a set of conditions were used as shown in Table 4. In this table, two conditions of the maximum stability ($\varphi =0$) and the maximum flexibility ($\varphi =1$) were used in each of which f, N, or Re were set to either maximum or minimum. In addition, for all the test subjects $D_1$ and $D_2$ were kept at 5 and 8 $\mu m$, respectively. The output of these tests indicates that the proposed NSGA3 could apply the constraints and find the design parameters efficiently. For example, in case 1, the f reached its minimum value with BW of 2.27 $\mu m$, while its corresponding case with the maximum flexibility (case 7) increased the BW up to 2.61 $\mu m$. This trend for BW occurred similarly where the cases with $\varphi =1$ always had bigger BWs. Moreover, as the last stage of checking the method’s accuracy, the proposed configurations were also simulated in the original software for comparison with the direct neural network’s prediction. The relative differences are shown in Table 4 as errors (E) being acceptably less than 4%.

Furthermore, to visually elaborate the automation outputs, Figure 5a depicts the range of $D_c$ based on $Re$. The results presented in these plots are aligned with the finding from Table 4 and it can be seen that the cases of maximum flexibility have a higher slope than the cases of maximum stability. In addition, the post-processing of case with minimum $N$ and maximum flexibility (case 9) is demonstrated in Figure 5b for 10 periods. In this figure, the particle trajectory plots illustrate how particle with bumped and zigzag modes would behave within the full length of the device and the recurrence maps provide the exact lateral positions of the particles at the input and output of each period.

Wall function was used to simulate particle interactions with obstacles. It stores the shortest distance from any position to pillars in a discretized periodic domain. Then, this information serves as an efficient mapping between the spatial positions of particles to their distance from pillars, since it is only calculated once for each geometry. Moreover, the normal and tangential unit vector fields for velocity modification in the contact point can be calculated by extracting the wall function gradient and its perpendicular direction, respectively. Figure A.1 shows an example of the wall distance function and the normal vector field.

The Newtonian root finding method was used for extracting the critical diameter from particle trajectory simulation. Here, the detailed explanation of this method is provided in the Algorithm 1 where it gets the desired tolerance as input and derives the critical diameter.

Figure C.1 depicts the relationship between all the effective parameters and critical diameter in the original data-set in four f values to provide a better understanding of the DLD physics and criteria that should be used in automation step. The overall trends of the critical diameter changes with an increase of f: at low f values, in almost all N values there is a peak around Re=5, followed by the critical diameter value decreases. Nevertheless, for N=3 there are no critical diameters at all. However, the peak vanishes as the f value grows and it only shows a positive slope. Moreover, at N=3 some area appears to have a critical diameter around a low Reynolds number and this area expands as f value increases. Furthermore, the area with a slope near zero offers the best stability since the critical diameter would not change much, and separation will occur regardless of the fluctuation in Re. On the other hand, it appears that higher f and lower N have a steeper slope and provide a better chance of flexible design.