Multi-Objective Optimizations of a Serpentine Micromixer with Crossing Channels at Low and High Reynolds Numbers

In order to maximize the mixing performance of a micromixer with an integrated three-dimensional serpentine and split-and-recombination configuration, multi-objective optimizations were performed at two different Reynolds numbers, 1 and 120, based on numerical simulation. Numerical analyses of fluid flow and mixing in the micromixer were performed using three-dimensional Navier-Stokes equations and convection-diffusion equation. Three dimensionless design variables that were related to the geometry of the micromixer were selected as design variables for optimization. Mixing index at the exit and pressure drop through the micromixer were employed as two objective functions. A parametric study was carried out to explore the effects of the design variables on the objective functions. Latin hypercube sampling method as a design-of-experiment technique has been used to select design points in the design space. Surrogate modeling of the objective functions was performed by using radial basis neural network. Concave Pareto-optimal curves comprising of Pareto-optimal solutions that represents the trade-off between the objective functions were obtained using a multi-objective genetic algorithm at Re = 1 and 120. Through the optimizations, maximum enhancements of 18.8% and 6.0% in mixing index were achieved at Re = 1 and 120, respectively.


Introduction
Microfluidics is related to an expeditiously emerging technology enabling manipulation and control of minute volumes of fluids with high accuracy in a miniaturized system for various fluidic functions, such as transporting, metering, valving, mixing, reacting, and separating [1,2]. Micromixer is an integral component of the microfluidic systems that have promising impact in the fields of biomedical diagnostics, drug development, and chemical industry [3,4]. Efficient mixing of liquid samples is a challenging task for successful operation of different processes in the microfluidic systems. The flow nature in the microfluidic systems is laminar, due to low Reynolds number. Thus, the mixing of fluid species depends mainly on mass diffusion in the absence of turbulence. However, the diffusion-dependent mixing is relatively slow and ineffective. In order to enhance the mixing performance, numerous methods have been proposed during the last two decades [3,[5][6][7].
Depending upon the working principle, micromixers are categorized either as a passive or as an active type. Active micromixers employ external energy sources, such as electrokinetic, ultrasonic vibration and magnetic field to generate flow perturbations inside the microchannel. Although active micromixers generally show excellent mixing capability and control during the mixing, high fabrication cost, and difficulty in integration with microfluidic systems make them less practical. In contrast, In the present work, multi-objective optimizations have been performed to further enhance the performance of the micromixer proposed by Raza et al. [25], at both low and high Reynolds numbers. Mixing index at the exit of the micromixer and pressure drop though the micromixer were employed as two objective functions. The Pareto-optimal fronts compromising these two objectives were generated using RBNN surrogate model and MOGA. The optimizations were carried out with the two objective functions at the two different Reynolds numbers, Re = 1 and 120. The corresponding Peclet numbers are 1 × 10 4 and 120 × 10 4 , respectively.

Micromixer Geometry
The micromixer configuration, as proposed by Raza et al. [25], has been used for the optimization. The micromixer consists of two repeated OX-shaped mixing units of 3D serpentine SAR structure as shown in Figure 1. The number of mixing units was reduced from five of previous work [25] to two in this work. Chaotic advection through stretching and folding of fluid interface at the X-junction, where the fluids in both the crossing channels are exchanged, enhances mixing at low Reynolds numbers. On the other hand, 3D serpentine path that was created by O-structure promotes mixing by producing chaotic advection at high Reynolds numbers, as explained in the previous study [25].
In the present work, multi-objective optimizations have been performed to further enhance the performance of the micromixer proposed by Raza et al. [25], at both low and high Reynolds numbers. Mixing index at the exit of the micromixer and pressure drop though the micromixer were employed as two objective functions. The Pareto-optimal fronts compromising these two objectives were generated using RBNN surrogate model and MOGA. The optimizations were carried out with the two objective functions at the two different Reynolds numbers, Re = 1 and 120. The corresponding Peclet numbers are 1 × 10 4 and 120 × 10 4 , respectively.

Micromixer Geometry
The micromixer configuration, as proposed by Raza et al. [25], has been used for the optimization. The micromixer consists of two repeated OX-shaped mixing units of 3D serpentine SAR structure as shown in Figure 1. The number of mixing units was reduced from five of previous work [25] to two in this work. Chaotic advection through stretching and folding of fluid interface at the X-junction, where the fluids in both the crossing channels are exchanged, enhances mixing at low Reynolds numbers. On the other hand, 3D serpentine path that was created by O-structure promotes mixing by producing chaotic advection at high Reynolds numbers, as explained in the previous study [25].
Values of the geometric parameters for the reference micromixer are listed in Table 1. The traditional method of stacking the polydimethylsiloxane (PDMS) layers can be used to fabricate the proposed micromixer. As described in the previous works [17,26], the micro-molding technology, i.e., soft lithography technique using SU-8 master molds, can be employed to fabricate these PDMS layers.
(a) (b) Figure 1. Geometry of micromixer with OX-shaped units: (a) planar view [25]; and, (b) three-dimensional (3D) view.  Values of the geometric parameters for the reference micromixer are listed in Table 1. The traditional method of stacking the polydimethylsiloxane (PDMS) layers can be used to fabricate the proposed micromixer. As described in the previous works [17,26], the micro-molding technology, i.e., soft lithography technique using SU-8 master molds, can be employed to fabricate these PDMS layers.

Numerical Analysis
The analyses of fluid flow and mixing inside the micromixer were performed numerically using a commercial CFD code, ANSYS CFX 15.0 (ANSYS, Inc., Canonsburg, PA, USA) [27], based on finite volume approximations. The numerical analysis of steady, incompressible, 3D laminar flow was carried out by solving the continuity and Navier-Stokes equations: where the symbols ν, ρ, and → V denote the fluid kinematic viscosity, density, and velocity, respectively. A solution of dye in water and water at a temperature of 20 • C were considered as two working fluids. To model the mass transport of the fluids having constant viscosity and density in the mixing process, a scalar transport equation of advection-diffusion type [28] was used, as follows: where the symbols α and C denote the diffusivity coefficient and dye concentration, respectively. In order to model mixing of two fluids, the scalar transport equation has been used and validated previously for different micromixers [29,30]. The computational domain was discretized by a combination of tetrahedral and hexahedral elements to reduce the total number of computational nodes. The mixing units were meshed with tetrahedral element due to their complex geometry, while the remaining part in the microchannel was meshed with hexahedral elements.
Uniform velocity profiles were assigned at the inlets, while the atmospheric pressure was used at the outlet. At the walls, no-slip condition was used. A solution of dye in water (mass fraction equal to 1) and water (mass fraction equal to 0) were introduced at the inlet 1 and inlet 2, respectively. The diffusivity coefficient value of the water-dye mixture was 1.0 × 10 −10 m 2 /s. The values of density and dynamic viscosity of water (and also water-dye mixture) were 1000 kg/m 3 and 1.0 × 10 −3 kg·m −1 ·s −1 , respectively [31]. Truncation errors that were associated with numerical discretization for the advection terms in the governing differential equations give rise to numerical diffusion. The extent of numerical diffusion depends on the numerical scheme used. Numerical diffusion can be decreased by using higher order approximation schemes, such as second-order upwind and third order QUICK [32] scheme instead of first-order upwind scheme [33]. A high-resolution scheme of second-order approximation was applied for discretization of advection terms. The solution convergence criterion was set as root-mean-square residual value of less than 1.0 × 10 −6 .
A statistical method [34] based on concept of intensity of segregation was used to define mixing index. The mixing index at a plane perpendicular to streamwise direction is represented as follows: where σ and σ max are the standard deviation of the concentration at the cross-sectional plane, and maximum standard deviation over the entire data range, respectively. The standard deviation of dye mass fraction at the cross-sectional plane can be written as: where c i , c m , and N are the mass fraction at sampling point i, the optimal mixing mass fraction, and the number of sampling points on the plane, respectively. A mixing index value of 0 indicates completely unmixed fluids, while a value of 1 indicates completely mixed fluids. Mixing index at the exit (M o ) was defined as the mixing index 700 µm downstream of the main channel starting position. The Reynolds Number was defined using hydraulic diameter (D h ) of the main channel as follows: where ρ, µ, and V denote the density of water, dynamic viscosity of water, and fluid average velocity, respectively. Hydraulic diameter was defined as follows: where W and H denote the width and depth of the microchannel.

Design Variables and Objective Functions
For the optimization, the ratios of distance between O-structures to pitch (d/Pi), crossing channel width to pitch (w 2 /Pi), and total depth of the micromixer to main channel width (H/W), as shown in Figure 1, were selected as design variables among various geometric parameters through a preliminary parametric study. Effects of these design variables on mixing index at the exit and pressure drop through the micromixer at Re = 1 for the reference micromixer, are shown in Figure 2. For all of the design variables, the mixing index shows maximum values in the tested ranges of the variables. Pressure drop increases with increase in d/Pi, while it decreases with increase in w 2 /Pi and H/W in the tested ranges. Design ranges for the design variables are summarized in Table 2.
where ρ, µ, and V denote the density of water, dynamic viscosity of water, and fluid average velocity, respectively. Hydraulic diameter was defined as follows: where W and H denote the width and depth of the microchannel.

Design Variables and Objective Functions
For the optimization, the ratios of distance between O-structures to pitch (d/Pi), crossing channel width to pitch (w2/Pi), and total depth of the micromixer to main channel width (H/W), as shown in Figure 1, were selected as design variables among various geometric parameters through a preliminary parametric study. Effects of these design variables on mixing index at the exit and pressure drop through the micromixer at Re = 1 for the reference micromixer, are shown in Figure 2. For all of the design variables, the mixing index shows maximum values in the tested ranges of the variables. Pressure drop increases with increase in d/Pi, while it decreases with increase in w2/Pi and H/W in the tested ranges. Design ranges for the design variables are summarized in Table 2.
In this work, mixing index at the exit and pressure drop in the micromixer were used to define two objective functions. The pressure drop that affects the power required to drive the fluids through the micromixer, was calculated as the difference in area-weighted average of the total pressure between the planes at the main channel inlet and exit. To make the optimization problem be defined as minimization of the objective functions, one of the objective functions (FM) was taken as the negative of mixing index, and the other objective function (F∆P) was defined as pressure drop. Two pairs of objective functions were selected as follows: mixing index at Re = 1 (FM at Re = 1)-pressure drop at Re = 1 (F∆P at Re = 1) and mixing index at Re = 120 (FM at Re = 120)-pressure drop at Re = 120 (F∆P at Re = 120). Values of the objective functions for the reference design are listed in Table 3.  In this work, mixing index at the exit and pressure drop in the micromixer were used to define two objective functions. The pressure drop that affects the power required to drive the fluids through the micromixer, was calculated as the difference in area-weighted average of the total pressure between the planes at the main channel inlet and exit. To make the optimization problem be defined as minimization of the objective functions, one of the objective functions (F M ) was taken as the negative of mixing index, and the other objective function (F ∆P ) was defined as pressure drop. Two pairs of objective functions were selected as follows: mixing index at Re = 1 (F M at Re = 1 )-pressure drop at Re = 1 (F ∆P at Re = 1 ) and mixing index at Re = 120 (F M at Re = 120 )-pressure drop at Re = 120 (F ∆P at Re = 120 ). Values of the objective functions for the reference design are listed in Table 3.

Surrogate Model and Multi-Objective Optimization
Procedure of multi-objective optimization used in the present work is represented in Figure 3. The first step is selecting design variables and their ranges through a previous parametric study, and the objective functions considering the design goals. The next step is to build the design space using the design of experiments (DOE). In the present work, Latin Hypercube Sampling (LHS) [35] is used as DOE. LHS is an effective sampling technique that uses an m × n simulation matrix where m is the number of levels (sampling points) to be examined and n is the number of design parameters. Each of the n columns of the matrix containing the levels, 1, 2, . . . , m, is randomly paired to form a Latin hypercube. This approach produces random sample points, ensuring that all of the portions of the design space are represented. Objective function values are calculated by Navier-Stokes analysis at these design points. A surrogate model is constructed to approximate the objective functions based on these objective function values, and MOGA is used to explore global Pareto-optimal solutions in the design space. If the global optimum solution exists in the design space and a termination criterion is satisfied, the multi-objective optimization procedure terminates.
The multi objective optimization problem was formulated, as follows: where F(x) is an objective function, x is a vector of design variables, and LB and UB denote the vectors of the lower and upper bounds, respectively.
is satisfied, the multi-objective optimization procedure terminates. The multi objective optimization problem was formulated, as follows: Design variable bound: LB ≤ x ≤ UB, x ∈ R  In order to obtain LHS design points, MATLAB (MathWorks, Inc., Natick, MA, USA) function 'lhsdesign' [36] was used with the criterion 'maxmin' (maximize minimum distance between adjacent design point). As a result, uniformly distributed 27 design points were selected for the three design variables in order to construct a surrogate model. The objective function values calculated at these design points by Navier-Stokes analysis are listed in Table 4. A RBNN model [37] was used to construct surrogate models of the objective functions. The RBNN model contains a two-layered network comprised of a hidden layer of radial neuron and an output layer of linear neuron. The hidden layer executes a nonlinear alteration of the input space to a middle space by using a set of radial basis elements. The output layer then implements a linear combiner to yield the desired targets. The linear model g(x) for the function can be expressed as a linear combination of a set of M radially-symmetric functions, as follows: where w j is weights and Φ j are radial basis functions. Benefit of this surrogate modeling is an ability to reduce computational time owing to the linear nature of the radial basis functions. In the present work, the RBNN function, newrb [36] was used to construct the models. The network training was performed by varying spread constant (SC) to adjust the cross-validation error. In this study, SC 1 is related to mixing index, and SC 2 is related to pressure drop. The SC values of each objective function were selected through a k-fold cross-validation [38] error test in which the errors were minimum at SC 1 = 0.9 and SC 2 = 0.5, as shown in Figure 4. The k-fold cross-validation is a validation method to estimate how the results of a statistical analysis generalize to an independent data set.
Micromachines 2018, 9, x 8 of 18 MOGA was used to find optimal solution on the constructed surrogate model using the MATLAB optimization tool box [36]. MOGA is a randomized global search method that solves functions by imitating progressions observed from natural evolution [37]. Based on the survival and growth of the fittest, MOGA finds new and improved results repeatedly. MOGA describes an initial population of individuals, which represent a part of the solution to the functions. Before the search begins, a series of chromosomes is randomly selected from the design space to obtain the initial population. Through computations, the individuals are selected in a competitive way, based on their fitness functions as measured by each specific objective function. The genetic search operators (i.e., "selection", "mutation", and "crossover") were applied to obtain a next generation of chromosomes, for which the predicted quality of all the chromosomes is better than that of the previous generation. This process is repeated until the termination criterion, which is function tolerance, is met. The following parameters were used: population size = 400, cross over fraction = 0.7, generations = 800, and function tolerance =10 −8 .

Results and Discussion
A grid-dependency test was performed to find out an optimal number of nodes for the spatial discretization of computational space. Four different grid systems with 4.78 × 10 5 to 1.43 × 10 6 nodes were tested for development of mixing index along the channel length at two different Reynolds numbers (Re = 1 and 120), as shown in Figure 5. The mixing indices were calculated on the planes perpendicular to the axial direction at four different locations (i.e., start of the main channel, two successive intersection nodes of the crossing channels in X-structure, and the exit). Almost similar profiles of the mixing index development are observed for grid systems with 1.21 × 10 6 and 1.43 × 10 6 nodes at both Re = 1 and 120. Hence, grid system with 1.21 × 10 6 nodes was selected as an optimal MOGA was used to find optimal solution on the constructed surrogate model using the MATLAB optimization tool box [36]. MOGA is a randomized global search method that solves functions by imitating progressions observed from natural evolution [37]. Based on the survival and growth of the fittest, MOGA finds new and improved results repeatedly. MOGA describes an initial population of individuals, which represent a part of the solution to the functions. Before the search begins, a series of chromosomes is randomly selected from the design space to obtain the initial population. Through computations, the individuals are selected in a competitive way, based on their fitness functions as measured by each specific objective function. The genetic search operators (i.e., "selection", "mutation", and "crossover") were applied to obtain a next generation of chromosomes, for which the predicted quality of all the chromosomes is better than that of the previous generation. This process is repeated until the termination criterion, which is function tolerance, is met. The following parameters were used: population size = 400, cross over fraction = 0.7, generations = 800, and function tolerance =10 −8 .

Results and Discussion
A grid-dependency test was performed to find out an optimal number of nodes for the spatial discretization of computational space. Four different grid systems with 4.78 × 10 5 to 1.43 × 10 6 nodes were tested for development of mixing index along the channel length at two different Reynolds numbers (Re = 1 and 120), as shown in Figure 5. The mixing indices were calculated on the planes perpendicular to the axial direction at four different locations (i.e., start of the main channel, two successive intersection nodes of the crossing channels in X-structure, and the exit). Almost similar profiles of the mixing index development are observed for grid systems with 1.21 × 10 6 and 1.43 × 10 6 nodes at both Re = 1 and 120. Hence, grid system with 1.21 × 10 6 nodes was selected as an optimal grid, commonly at Re = 1 and 120. The present numerical model was validated in the previous study [25] by both qualitative and quantitative comparisons of the numerical results with the experimental results of Hossain et al. [39], as shown in Figure 6. The optical image of fluid mixing in the first mixing unit at Re = 60 is compared with the numerical result of dye mass fraction distribution on x-y plane located halfway along the channel depth, as shown in Figure 6a. The numerical values of mixing indices at the exit of the micromixer in a Reynolds number range of 0.2-120 are also compared with the experimental data, as shown in Figure 6b. The numerically predicted mixing indices are slightly higher than the experimental data over the whole range. The uncertainties that are involved in the experimental procedures, such as capturing and analyzing experimental image, geometrical variation in fabrication, and wall roughness, can be attributed as the causes for the differences in the mixing indices, as shown in Figure 6b. However, the qualitative and quantitative comparisons between the numerical and experimental results, show acceptable agreements.  The present numerical model was validated in the previous study [25] by both qualitative and quantitative comparisons of the numerical results with the experimental results of Hossain et al. [39], as shown in Figure 6. The optical image of fluid mixing in the first mixing unit at Re = 60 is compared with the numerical result of dye mass fraction distribution on x-y plane located halfway along the channel depth, as shown in Figure 6a. The numerical values of mixing indices at the exit of the micromixer in a Reynolds number range of 0.2-120 are also compared with the experimental data, as shown in Figure 6b. The numerically predicted mixing indices are slightly higher than the experimental data over the whole range. The uncertainties that are involved in the experimental procedures, such as capturing and analyzing experimental image, geometrical variation in fabrication, and wall roughness, can be attributed as the causes for the differences in the mixing indices, as shown in Figure 6b. However, the qualitative and quantitative comparisons between the numerical and experimental results, show acceptable agreements.  Figure 7 shows effect of number of the mixing units on the mixing performance through a quantitative comparison among the mixing indices at the exits of the micromixers with two, three, and five mixing units at different Reynolds numbers. It is observed that the micromixer with five mixing units exhibits almost complete mixing in most of the Reynolds number range, while the micromixers with two and three mixing units show largely reduced mixing index values over the whole Re range. It is also observed that all the micromixers show minimum mixing indices at Re = 1. Hence, the micromixer with two mixing units was selected in the present work to have more space for enhancement of the mixing performance by optimization. Figure 6. Validation of numerical results compared with experimental data [25]: (a) qualitative comparison of dye mass fraction distribution on x-y plane located halfway along the channel depth between experiment and numerical analysis at Re = 60; and, (b) quantitative comparison of mixing indices at the exit for different Reynolds numbers between experiment and numerical analysis. Figure 7 shows effect of number of the mixing units on the mixing performance through a quantitative comparison among the mixing indices at the exits of the micromixers with two, three, and five mixing units at different Reynolds numbers. It is observed that the micromixer with five mixing units exhibits almost complete mixing in most of the Reynolds number range, while the micromixers with two and three mixing units show largely reduced mixing index values over the whole Re range. It is also observed that all the micromixers show minimum mixing indices at Re = 1. Hence, the micromixer with two mixing units was selected in the present work to have more space for enhancement of the mixing performance by optimization.
Effect of Reynolds number on development of mixing index along the axial length of the micromixer (reference design) is shown in Figure 8. The mixing indices were calculated on four y-z planes (at the start of the main channel, the two successive intersection nodes of the crossing channels in X-structure, and the exit) for Reynolds numbers, 1, 20, 40, 80, and 120. The results show that developing rate of mixing index depends on Reynolds number, as shown in Figure 8. The rate of development of mixing index is found to be highest between the first and second intersection nodes at Re = 40. mixing units exhibits almost complete mixing in most of the Reynolds number range, while the micromixers with two and three mixing units show largely reduced mixing index values over the whole Re range. It is also observed that all the micromixers show minimum mixing indices at Re = 1. Hence, the micromixer with two mixing units was selected in the present work to have more space for enhancement of the mixing performance by optimization. planes (at the start of the main channel, the two successive intersection nodes of the crossing channels in X-structure, and the exit) for Reynolds numbers, 1, 20, 40, 80, and 120. The results show that developing rate of mixing index depends on Reynolds number, as shown in Figure 8. The rate of development of mixing index is found to be highest between the first and second intersection nodes at Re = 40. Following the procedure outlined in Figure 3, Pareto-optimal fronts presenting the optimum trade-off between the two conflicting objectives are plotted, as shown in Figure 9. Pareto-optimal fronts, termed as POF-1 and POF-2 in the Figure 9a,b, respectively, represent two pairs of the objective functions at two different Reynolds numbers: FM at Re = 1 − F∆P at Re = 1 and FM at Re = 120 − F∆P at Re = 120. Concave shape of the Pareto-optimal fronts indicates that an improvement in mixing index occurs with simultaneous increase in pressure drop. Selection of Pareto-optimal solution depends upon the choice of the designer, since each solution is a global Pareto-optimal solution and none of these Pareto-optimal solutions is superior to the others for both objective functions. Following the procedure outlined in Figure 3, Pareto-optimal fronts presenting the optimum trade-off between the two conflicting objectives are plotted, as shown in Figure 9. Pareto-optimal fronts, termed as POF-1 and POF-2 in the Figure 9a,b, respectively, represent two pairs of the objective functions at two different Reynolds numbers: F M at Re = 1 − F ∆P at Re = 1 and F M at Re = 120 − F ∆P at Re = 120 . Concave shape of the Pareto-optimal fronts indicates that an improvement in mixing index occurs with simultaneous increase in pressure drop. Selection of Pareto-optimal solution depends upon the choice of the designer, since each solution is a global Pareto-optimal solution and none of these Pareto-optimal solutions is superior to the others for both objective functions. objective functions at two different Reynolds numbers: FM at Re = 1 − F∆P at Re = 1 and FM at Re = 120 − F∆P at Re = 120. Concave shape of the Pareto-optimal fronts indicates that an improvement in mixing index occurs with simultaneous increase in pressure drop. Selection of Pareto-optimal solution depends upon the choice of the designer, since each solution is a global Pareto-optimal solution and none of these Pareto-optimal solutions is superior to the others for both objective functions.  In order to analyze the Pareto-optimal solutions, five representative Pareto-optimal designs (PODs) were selected by K-means clustering on each Pareto-optimal front, as shown in Figure 9. PODs A and E at extreme ends of each Pareto-optimal front represent pressure drop-oriented and mixing index-oriented designs, respectively. Accomplishment of one objective function leads to forfeit of the other objective function. At Re = 1, as compared to POD A on POF-1, the mixing index-oriented design, POD E shows a relative enhancement of 34.5 % in the mixing index at the exit. The pressure drop-oriented design, POD A shows a 64.5 % reduction in pressure drop as compared to POD E. Similarly, at Re = 120, POD E shows a relative enhancement of 2.8 % in mixing index at the exit as compared to POD A on POF-2, and POD A shows a 65.7% reduction in pressure drop as compared to POD E. These results reveal that the relative enhancement of mixing index at the exit is much more pronounced at Re = 1 as compared to Re = 120. However, the relative percentage reductions in pressure drop are not much different at both the Reynolds numbers. The higher relative percentage changes in pressure drop than those in mixing index at the exit indicate that mixing index is less sensitive to the design variables as compared to pressure drop.
The results at the representative PODs for POF-1 and POF-2 are listed in Tables 5 and 6, respectively. It is observed that mixing index at the exit and pressure drop are most sensitive to the design variable, H/W, while w 2 /Pi remains nearly invariant on POF-1 with the exception of POD A. A high value of mixing index is observed for d/Pi, w 2 /Pi, and H/W values close to middle of the range, upper bound, and lower bound, respectively. In case of POF-2, the objective functions become more sensitive to the design variables when compared to POF-1. The results of numerical analysis at the PODs shows that maximum enhancements of 18.8% (POD E on POF-1) and 6.0% (POD E on POF-2) in mixing index at the exit are achieved as compared to the reference design at Re = 1 and 120, respectively, by the optimization. Maximum reductions of 5.8% (POD A on POF-1) and 11.1% (POD A on POF-2) in pressure drop are obtained as compared to the reference design at Re = 1 and 120, respectively. The surrogate model predictions of the objective functions values are also compared with the numerical results at the same designs in Tables 5 and 6. The relative errors of the surrogate predictions for mixing index at the exit and pressure drop are less than 2% and 10%, respectively, at Re = 1, as shown in Table 5. However, these relative errors are increased at Re = 120 to less than 5% and 22%, respectively, as shown in Table 6. For qualitative comparison of mixings at different PODs on POF-1 (Re = 1) and POF-2 (Re = 120), dye mass fraction contours are plotted on y-z planes at the beginning of the main channel (x/L t = 0), two successive intersection nodes of the crossing channel in X-structure (x/L t = 0.13 and 0.26), and the exit (x/L t = 0.32), respectively, as shown in Figure 10. Two Pareto-optimal designs, i.e. POD A and POD E at the extreme ends of each Pareto-optimal front are selected for the comparison. Number of two-fluid interfaces increases along the channel length for both the PODs. However, there is a distinct difference in the dye mass fraction distribution at the exit of the micromixer between POD A and POD E on POF-1. A more uniform dye mass fraction distribution is observed at the exit for POD E as compared to POD A on POF-1, as shown in Figure 10a, whereas almost similar distributions are observed for POD A and POD E on POF-2 ( Figure 10b). This is related to the fact that the change in mixing index along POF-1 (0.439-0.581 in Table 5) is much larger than that along POF-2 (0.936-0.957 in Table 6).
Micromachines 2018, 9, x 13 of 18 difference in the dye mass fraction distribution at the exit of the micromixer between POD A and POD E on POF-1. A more uniform dye mass fraction distribution is observed at the exit for POD E as compared to POD A on POF-1, as shown in Figure 10a, whereas almost similar distributions are observed for POD A and POD E on POF-2 ( Figure 10b). This is related to the fact that the change in mixing index along POF-1 (0.439-0.581 in Table 5) is much larger than that along POF-2 (0.936-0.957 in Table 6). In order to analyze the mixing mechanism, velocity vectors superimposed on the dye mass fraction contours are plotted for POD A and POD E on POF-1 (Re = 1), in Figure 11. It is observed that a saddle-shaped flow structure, which promotes chaotic advection [25], is formed at each intersection node (x/Lt = 0.13 and 0.26) of the crossing microchannels in both of the Pareto-optimal designs, POD A and POD E. However, there is a difference in the fraction of total depth of the microchannel cross-section covered by saddle-shaped pattern between the two designs. In case of POD A, the saddle-shaped flow structure exists only in the middle height of the y-z plane, as shown in Figure 11a. Whereas, in case of POD E, the saddle-shaped flow pattern covers more height of the y-z In order to analyze the mixing mechanism, velocity vectors superimposed on the dye mass fraction contours are plotted for POD A and POD E on POF-1 (Re = 1), in Figure 11. It is observed that a saddle-shaped flow structure, which promotes chaotic advection [25], is formed at each intersection node (x/L t = 0.13 and 0.26) of the crossing microchannels in both of the Pareto-optimal designs, POD A and POD E. However, there is a difference in the fraction of total depth of the microchannel cross-section covered by saddle-shaped pattern between the two designs. In case of POD A, the saddle-shaped flow structure exists only in the middle height of the y-z plane, as shown in Figure 11a. Whereas, in case of POD E, the saddle-shaped flow pattern covers more height of the y-z plane, as shown in Figure 11b. This results in more uniform mass fraction distribution ( Figure 11b) and consequently, higher mixing index in POD E, as compared to POD A. Figure 10. Dye mass fraction distributions on y-z planes for POD A and POD E: (a) Pareto-optimal fronts (POF)-1 (Re = 1); and (b) POF-2 (Re = 120).
In order to analyze the mixing mechanism, velocity vectors superimposed on the dye mass fraction contours are plotted for POD A and POD E on POF-1 (Re = 1), in Figure 11. It is observed that a saddle-shaped flow structure, which promotes chaotic advection [25], is formed at each intersection node (x/Lt = 0.13 and 0.26) of the crossing microchannels in both of the Pareto-optimal designs, POD A and POD E. However, there is a difference in the fraction of total depth of the microchannel cross-section covered by saddle-shaped pattern between the two designs. In case of POD A, the saddle-shaped flow structure exists only in the middle height of the y-z plane, as shown in Figure 11a. Whereas, in case of POD E, the saddle-shaped flow pattern covers more height of the y-z plane, as shown in Figure 11b. This results in more uniform mass fraction distribution ( Figure 11b) and consequently, higher mixing index in POD E, as compared to POD A. In order to have an idea about the mixing performance of a Pareto-optimal design at Reynolds numbers other than the Reynolds number at which it was obtained, mixing performances of PODs E on POF-1 and POF-2 are compared with the reference design at Reynolds numbers, 1, 30, 60, 90, and In order to have an idea about the mixing performance of a Pareto-optimal design at Reynolds numbers other than the Reynolds number at which it was obtained, mixing performances of PODs E on POF-1 and POF-2 are compared with the reference design at Reynolds numbers, 1, 30, 60, 90, and 120, as shown in Figure 12. It is observed that POD E on POF-1 shows values of mixing index larger than those of the reference design except at intermediate Reynolds number of 30, but the mixing index values are smaller than those of POD E on POF-2 except at Re = 1. POD E on POF-2 shows improvements in the mixing index at Reynolds numbers larger than 60 when compared to the other designs.  In order to compare the merit of PODs at Re = 1 with previous micromixers, as shown in Table 7, two different mixing performance parameters named, "mixing cost" [40], and "mixing energy cost" (mec) [41] have been used. The terms, "mixing cost" and "mixing energy cost" are defined as: In order to compare the merit of PODs at Re = 1 with previous micromixers, as shown in Table 7, two different mixing performance parameters named, "mixing cost" [40], and "mixing energy cost" (mec) [41] have been used. The terms, "mixing cost" and "mixing energy cost" are defined as: where η, ∆P, and C p denote mixing efficiency, pressure drop and mean input power coefficient, respectively. The detail description about these terms can be found in the works of Ortega-Casanova [41,42].
Mixing efficiency (η) and mean input power coefficient (C p ) are defined as follows: where σ, σ max , ∆P, q, ρ, and V are the standard deviation of the concentration, maximum standard deviation, pressure drop, dimensionless flow rate, density of fluid, and average flow velocity, respectively. A micromixer with higher mixing cost and lower mec indicates a more efficient micromixer. From Table 7, it is evident that PODs have neither the best nor the worst values of these parameters. Higher mixing cost values are shown by all of the PODs as compared to the micromixer proposed by Cheri et al. [43]. PODs show higher mec values as compared to the rhombic micromixer with asymmetrical flow [40]. The highest mec value is shown by the micromixer configuration of a rectangular chamber with obstruction proposed by Cheri et al. [43]. Although the micromixer proposed by Cheri et al. [43] shows the highest pressure drop, due to the highest mixing efficiency and the lowest mean input power coefficient, the mec value becomes least.

Conclusions
Multi-objective optimizations of a micromixer with 3D serpentine and SAR configuration have been performed at Reynolds numbers, 1 and 120, based on flow and mixing analyses using 3D Navier-Stokes equations and convection-diffusion equation. Three design variables, i.e., d/Pi, w 2 /Pi, and H/W, were selected, and two objective functions were defined in terms of mixing index at the exit of the micromixer and pressure drop through the micromixer for the optimization. In a parametric study, the mixing index shows maxima for all of the design variables, but the pressure drop shows monotonic variations for all the design variables in the tested ranges. Two concave Pareto-optimal fronts (POF-1 and POF-2) representing trade-off between the two objective functions at Re = 1 and 120, respectively, were obtained by RBNN surrogate model and MOGA. In applying RBNN model, it was found that the k-hold cross-validation errors for mixing index and the pressure drop were minimized with the spread constants, 0.9 and 0.5, respectively. On POF-1 (Re = 1), the preference of a mixing index-oriented design, POD E over a pressure drop-oriented design, POD A leads to 34.5% relative increase in mixing index at the exit, and the preference of POD A over POD E showed 64.5% reduction in pressure drop. On POF-2 (Re = 120), the preference of POD E over POD A leads to only 2.8% relative increase in mixing index at the exit, and the preference of POD A over POD E showed 65.7% reduction in pressure drop. These results indicate that mixing index is less sensitive to the design variables as compared to pressure drop, on the Pareto-optimal fronts. It was found from the numerical analysis at the PODs that the maximum enhancements of 18.8% at POD E on POF-1 and 6.0% at POD E on POF-2 in mixing index at the exit were obtained when compared to the reference design. And, maximum reductions of 5.8% at POD A on POF-1 and 11.1% at POD A on POF-2 in pressure drop were achieved compared to the reference design. Maximum relative error of the surrogate prediction compared to the numerical analysis was smaller for mixing index than pressure drop, and increased with the Reynolds number. In a range of Re = 1-120, POD E on POF-1 showed values of mixing index larger than those of the reference design except at Re = 30, but POD E on POF-2 showed higher mixing performance than the reference design, only for the Reynolds numbers larger than 60.