Single and Multi-Objective Optimization of a Three-Dimensional Unbalanced Split-and-Recombine Micromixer

The three-dimensional geometry of a micromixer with an asymmetrical split-and-recombine mechanism was optimized to enhance the fluid-mixing capability at a Reynolds number of 20. Single and multi-objective optimizations were carried out by using particle swarm optimization and a genetic algorithm on a modeled surrogate surface. Surrogate modeling was performed using the computational results for the mixing. Mixing and flow analyses were carried out by solving the convection–diffusion equation in combination with the three-dimensional continuity and momentum equations. The optimization was carried out with two design variables related to dimensionless geometric parameters. The mixing effectiveness was chosen as the objective function for the single-objective optimization, and the pressure drop and mixing index at the outlet were chosen for the multi-objective optimization. The sampling points in the design space were determined using a design of experiment technique called Latin hypercube sampling. The surrogates for the objective functions were developed using a Kriging model. The single-objective optimization resulted in 58.9% enhancement of the mixing effectiveness compared to the reference design. The multi-objective optimization provided Pareto-optimal solutions that showed a maximum increase of 48.5% in the mixing index and a maximum decrease of 55.0% in the pressure drop in comparison to the reference design.


Introduction
Microfluidics is rapidly emerging for precise control and manipulating fluids in microscale channels, and has found a wide range of applications in bioengineering, the chemical industry, and environmental monitoring [1][2][3]. Most applications of microfluidic systems involve the mixing of two or more species as a fundamental process. Successful completion of the applications requires rapid and efficient mixing, but microfluidic systems exhibit laminar flow characteristics due to their operation at low Reynolds numbers. The mixing achieved through molecular diffusion in the absence of turbulence is very slow and inefficient for microfluidic applications. Hence, it is crucial to attain efficient mixing to enhance the performance of microfluidic devices.
Computational fluid dynamics (CFD) is extensively employed as a tool to evaluate the performance and the flow structure in micromixers. Over the past two decades, various passive micromixer designs [4,5] have been proposed based on mixing strategies such as parallel/sequential lamination [6,7], hydrodynamic focusing [8], and chaotic advection [9]. Combinations of these strategies can also be employed [10,11]. The sequential processes of splitting, rejoining, and rearranging the flow in The dimensions of the reference design are the same as in previous work [17]: The axial length of the initial part of the main channel (L0), exit channel length (Le), total length (Lt), pitch length (Pi), width of the main channel (W), width of the major sub-channel (w1), width of the minor sub-channel (w2), width of the dislocation (w3), depth of the major sub-channel (h1), depth of the minor sub-channel (h2), and total depth of the main channel (H) are 500 µ m, 2950 µ m, 6750 µ m, 1200 µ m, 300 µ m, 200 µ m,100 µ m, 100 µ m, 100 µ m, 100 µ m, and 200 µ m, respectively. To analyze the mixing performance, the CFD code ANSYS CFX 15.0 ® [23] was used for numerical analyses of the mass transport and fluid flow inside the micromixer model. The CFD code is based on finite-volume approximations and finds the solutions to the momentum and continuity equations for the Newtonian, incompressible, and steady 3D laminar flow: where ⃗ , ν, and ρ are the fluid velocity, kinematic viscosity, and density, respectively. A dye-water solution and water enter inlets 1 and 2, respectively. Mass transport occurs in mixing due to diffusion and advection. Therefore, an advection-diffusion equation [24] was used to simulate the mass transport of fluids with constant properties during mixing: The dimensions of the reference design are the same as in previous work [17]: The axial length of the initial part of the main channel (L 0 ), exit channel length (L e ), total length (L t ), pitch length (Pi), width of the main channel (W), width of the major sub-channel (w 1 ), width of the minor sub-channel (w 2 ), width of the dislocation (w 3 ), depth of the major sub-channel (h 1 ), depth of the minor sub-channel (h 2 ), and total depth of the main channel (H) are 500 µm, 2950 µm, 6750 µm, 1200 µm, 300 µm, 200 µm,100 µm, 100 µm, 100 µm, 100 µm, and 200 µm, respectively. To analyze the mixing performance, the CFD code ANSYS CFX 15.0 ® [23] was used for numerical analyses of the mass transport and fluid flow inside the micromixer model. The CFD code is based on finite-volume approximations and finds the solutions to the momentum and continuity equations for the Newtonian, incompressible, and steady 3D laminar flow: ∇· where → V , ν, and ρ are the fluid velocity, kinematic viscosity, and density, respectively. A dye-water solution and water enter inlets 1 and 2, respectively. Mass transport occurs in mixing due to diffusion and advection. Therefore, an advection-diffusion equation [24] was used to simulate the mass transport of fluids with constant properties during mixing: where C and D are the concentration and diffusivity coefficients, respectively. The diffusion coefficient was 1.2 × 10 −9 m 2 /s. The properties of the two fluids are set to that of water at 25 • C (density ρ = 997 kg/m 3 and viscosity µ = 0.00089 kg·m −1 ·s −1 ). Reynolds number was defined as Re = ρUW/µ where U is average velocity in the inlet channels. The micromixer domain was discretized using a tetrahedral grid system. ANSYS ICEM CFD 15.0 ® [23] was used to generate the mesh. Due to the unstructured tetrahedral elements, it is not possible to control the number of grid points in a specific direction. However, the average number of grid nodes used along H (Figure 1b) was 36. The boundary conditions were specified as follows: identical uniform velocities at the inlets, zero velocity at the walls, and zero static pressure at the outlet. For species transport, the molar concentration fractions of the dye-water solution were assigned a value of 1 at inlet 1 and 0 at inlet 2. Zero mass flux was employed at the channel walls.
The accuracy of the numerical results of mixing in the microchannel is sensitive to the discretization scheme. The discretization of the advection terms usually produces artificial diffusion due to truncation errors. The artificial diffusion cannot be fully ruled out but can be curtailed by employing certain methods [25]. Higher-order schemes such as the third-order QUICK [26] and second-order upwind schemes are less prone to numerical diffusion, in contrast to the first-order upwind scheme. Therefore, the second-order high-resolution scheme [23] was chosen for the discretization of the advection terms. A converged solution was obtained when the root-mean-square residual value of each variable was less than 1.0 × 10 −6 . The distribution of species at a plane perpendicular to the streamwise direction was used to measure the mixing capability. The variance of the mass fraction of species based on the concept of the intensity of segregation [27] was defined as follows: where N, c i , and c m are the number of data points, the mass fraction at data point i, and the optimal mixing mass fraction, respectively. The mixing index at the plane was derived as follows: where σ max is the maximum variance over the entire data range. A higher mixing index means a higher concentration uniformity. Therefore, for completely unmixed fluids, the mixing index is zero, while for completely mixed fluids, the value is unity. The mixing index at the exit (M o ) was calculated on a plane located at 6700 µm from the main channel entry.

Design Variables and Objective Functions
The selection of appropriate design variables that largely affect the design objectives is an important step in the optimization process. The chosen design variables were the ratio of the width of the dislocation to the width of the major sub-channel (w 3 /w 1 ) and the ratio of the depth of the major sub-channel to the depth of the minor sub-channel (h 1 /h 2 ), as shown in Figure 1. The ranges for these variables were selected based on a preliminary parametric study. The parameters W, H, w 1 , and w 2 were kept constant throughout the optimization. Figure 2 shows the responses of the mixing and pressure drop to the changes of these design variables in the reference design of the micromixer at Re = 20. Variations of 128.5% and 26% in the mixing index are noticed with changes in w 3 /w 1 and h 1 /h 2 , respectively, within the tested ranges of the design variables. Similarly, there were variations of 110.2% and 88.5% in the pressure drop with changes in w 3 /w 1 and h 1 /h 2 , respectively. Table 1 lists the design variables and their ranges.  The objective functions for the optimization are the mixing index at the exit, mixing effectiveness, and pressure drop through the micromixer. The pressure drop is related to the pumping power necessary to drive the flow of mixing species and was determined as the difference between the area-averaged pressures at the inlet and exit of the main channel. The mixing effectiveness takes into account both the mixing index at the exit and the pressure drop in a single parameter as follows: where the pressure loss coefficient (K) is calculated as follows: U, ρ, and ΔP are the average inlet velocity, density, and pressure drop, respectively. The present work suggests two different approaches to deal with two objectives, i.e., enhancing mixing index and reducing pressure drop, in the optimization: one is to combine these two objectives into a single-objective function and to perform a single-objective optimization, and the other is to perform a multi-objective optimization with two objective functions related to these objectives. Therefore, two different optimizations were performed in this work: single-objective optimization for maximum mixing effectiveness and multi-objective optimization to find optimum solutions between the mixing index and the pressure drop.
The single-objective optimization uses an objective function FME based on the mixing efficiency as follows: The negative sign was introduced to define the optimization problem as a minimization of the objective function. The multi-objective optimization was done using objective functions FMI and FΔP, which were defined using the mixing index at the exit (M0) and pressure drop ( P), respectively:  The objective functions for the optimization are the mixing index at the exit, mixing effectiveness, and pressure drop through the micromixer. The pressure drop is related to the pumping power necessary to drive the flow of mixing species and was determined as the difference between the area-averaged pressures at the inlet and exit of the main channel. The mixing effectiveness takes into account both the mixing index at the exit and the pressure drop in a single parameter as follows: where the pressure loss coefficient (K) is calculated as follows: U, ρ, and ∆P are the average inlet velocity, density, and pressure drop, respectively. The present work suggests two different approaches to deal with two objectives, i.e., enhancing mixing index and reducing pressure drop, in the optimization: one is to combine these two objectives into a single-objective function and to perform a single-objective optimization, and the other is to perform a multi-objective optimization with two objective functions related to these objectives. Therefore, two different optimizations were performed in this work: single-objective optimization for maximum mixing effectiveness and multi-objective optimization to find optimum solutions between the mixing index and the pressure drop. The single-objective optimization uses an objective function F ME based on the mixing efficiency as follows: The negative sign was introduced to define the optimization problem as a minimization of the objective function. The multi-objective optimization was done using objective functions F MI and F ∆P , which were defined using the mixing index at the exit (M 0 ) and pressure drop (∆P), respectively: The values of different objective functions attained for the reference geometry are shown in Table 2.  Figure 3 shows the optimization process used for the single and multi-objective optimizations based on surrogate modeling. The optimization problem was formulated as follows:

Single and Multi-objective Optimization Methods
where F(x) represents the objective functions, x is the vector of design variables, and LB and UB are the lower and upper bounds of the design variables, respectively [28].
DOE was used to select arbitrary points in the design space to construct surrogate models. DOE methods can be generally classified into two categories: orthogonal design and random design. The orthogonality of a design means that the model parameters are statistically independent. This indicates that the factors in an experiment are uncorrelated and can be varied independently.
A factorial design, represented by the orthogonality design, has some disadvantages: initially it is usually unclear which factor is important. Since the underlying function is deterministic, there is a possibility that some of the initial design points collapse and one or more of the time-consuming computer experiments become useless. This issue is called the collapse problem [29]. Most classic DOEs, including the factorial design, are only applicable to rectangular design regions. To overcome this problem, a random design is used. A random design means that the model parameter values for the experiments are assigned on the basis of a random process [30]. Thus, the collapse problem does not occur with the random design. This is because if one or more factors appear not to be important, every point in the design still provides some information regarding the influence of the other factors on the response. In this way, none of the time-consuming computer experiments will turn out to be useless [30].
Therefore, in the present work, LHS [31] as the random design was used to construct surrogate models which were created to approximate the objective functions. 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 portions of the design space are represented. A genetic algorithm (GA) [32] and Particle Swarm Optimization (PSO) [33] were used as searching algorithms to find the optimum points on the surrogate models in the design space. The MATLAB function lhsdesign was used to obtain the design points, and maxmin was used to maximize the minimum distance between adjacent design points [34]. LHS selected 12 design points for the two design variables, and the objective functions values at the design points were calculated numerically, as shown in Table 3.  The KRG model [35] was used for surrogate modeling of the objective functions. The KRG model is expressed by the unknown function y(x) as follows: where x is an m-dimensional vector (for m design variables). f(x) is a known function representing the tendency for the design space and is generally expressed as a polynomial (also called a global model). Z(x) denotes the local deviations from the global model. In the KRG model, the local deviation at an unknown point can be expressed using a stochastic process. The sample points are interpolated with a Gaussian random function as a correlation function to estimate the trend of the stochastic process. The covariance matrix of Z(x) is given by:  The KRG model [35] was used for surrogate modeling of the objective functions. The KRG model is expressed by the unknown function y(x) as follows: where x is an m-dimensional vector (for m design variables). f(x) is a known function representing the tendency for the design space and is generally expressed as a polynomial (also called a global model). Z(x) denotes the local deviations from the global model. In the KRG model, the local deviation at an unknown point can be expressed using a stochastic process. The sample points are interpolated with a Gaussian random function as a correlation function to estimate the trend of the stochastic process. The covariance matrix of Z(x) is given by: where R is a correlation matrix with a spatial correlation function (SCF) and R x i , x j as its elements. σ 2 is the process variance, which is a scalar of the SCF that quantifies the correlation between any two n s sampled data points x i and x j to control the smoothness of the KRG model. The Gaussian function is the most desirable SCF when used with gradient-based optimization algorithms because it provides a relatively smooth and infinitely differentiable surface [36]. A previous study [37] showed that PSO can find an optimal point with better performance in a single-objective optimization than GA. On the other hand, GA is suitable for relatively complex optimization problems such as multi-objective optimizations because it assesses various points in the design space. Therefore, PSO and GA were used for the single-and multi-objective optimizations, respectively. In PSO, the velocities and positions of the particles are constantly updated with the particles that previously had the best performance [33]. GA is a random global searching technique that finds an optimum for the fitness function by mimicking the evolutions observed in natural phenomena [32]. Based on growth and survival of the fittest, GA repeatedly searches for improved results.
MOGA was used to obtain Pareto-optimal solutions using the customized MATLAB function, "gamultiobj" [34]. The optimization used the GA to search the Pareto-optimal solutions in a feasible solution space bounded by the lower and upper bounds on the design variables. The function approximations for the two objectives, F MI and F ∆p , were supplied in vector form. GA uses three main types of operation to create future generation from the current: selection, crossover and mutations [34]. In the single-objective optimization, the superiority of a solution over other solutions was easily determined by comparing their objective function values. However, in the multi-objective optimization problem, the goodness of a solution was determined by the dominance. Given the multi-objective optimization results, the non-dominated solution set of the entire feasible space is called the Pareto-optimal set and the boundary from the Pareto-optimal set is called the Pareto-optimal front. A solution is called non-dominated if none of the objective functions can be improved in value without degrading some of the other objective values. The function, gamultiobj, used in the present multi-objective optimization, was also coded to determine the dominance of solutions [38]. To identify among the Pareto-optimal solutions, these solutions were grouped into clusters using K-mean clustering [39].

Results and Discussion
A mesh test was performed to find the best grid system that ensures that the numerical results do not change upon further refinement. A mesh with 9.23 × 10 6 nodes was chosen as an optimum mesh in a wide range of node numbers from 0.10 × 10 6 to 1.41 × 10 7 , as shown in Figure 4. The mixing index at the microchannel exit for the selected mesh and a finer mesh show an insignificant difference of 0.52%.
Considering the methods outlined in the ref. [25], numerical diffusion was evaluated for the selected mesh. The numerical diffusion coefficient was 1.04 × 10 −8 m 2 /s. Although this value is higher than the value of the molecular diffusion coefficient, the grid on further refinement shows a negligible change in the mixing index. This indicates that the mixing index changes negligibly with further reduction in the numerical diffusion, because the numerical diffusion coefficient is proportional to the mesh size. This may be caused by the dominance of chaotic advection over diffusion, and mixing is not attributed solely to diffusion. Furthermore, based on the numerical diffusion vs. mesh density graph [40], the required number of nodes to achieve the numerical diffusion coefficient equal or less than the molecular diffusion coefficient is approximately 3.41 × 10 8 for this problem. Considering the computational limit, this is not possible at present. multi-objective optimization, was also coded to determine the dominance of solutions [38]. To identify among the Pareto-optimal solutions, these solutions were grouped into clusters using Kmean clustering [39].

Results and Discussion
A mesh test was performed to find the best grid system that ensures that the numerical results do not change upon further refinement. A mesh with 9.23 × 10 6 nodes was chosen as an optimum mesh in a wide range of node numbers from 0.10 × 10 6 to 1.41 × 10 7 , as shown in Figure 4. The mixing index at the microchannel exit for the selected mesh and a finer mesh show an insignificant difference of 0.52%. Considering the methods outlined in the ref. [25], numerical diffusion was evaluated for the selected mesh. The numerical diffusion coefficient was 1.04 × 10 −8 m 2 /s. Although this value is higher The accuracy of the numerical simulation was verified by comparing the results with experimental data from Li et al. [16], as shown in Figure 5. The numerical plot of the dye mass fraction rendering of the micromixer was compared with the optical image of the dye mass fraction at Re = 80. There is acceptable agreement between the two concentration distributions, which validates the numerical approach.
Micromachines 2019, 10, x 9 of 17 than the value of the molecular diffusion coefficient, the grid on further refinement shows a negligible change in the mixing index. This indicates that the mixing index changes negligibly with further reduction in the numerical diffusion, because the numerical diffusion coefficient is proportional to the mesh size. This may be caused by the dominance of chaotic advection over diffusion, and mixing is not attributed solely to diffusion. Furthermore, based on the numerical diffusion vs. mesh density graph [40], the required number of nodes to achieve the numerical diffusion coefficient equal or less than the molecular diffusion coefficient is approximately 3.41 × 10 8 for this problem. Considering the computational limit, this is not possible at present. The accuracy of the numerical simulation was verified by comparing the results with experimental data from Li et al. [16], as shown in Figure 5. The numerical plot of the dye mass fraction rendering of the micromixer was compared with the optical image of the dye mass fraction at Re = 80. There is acceptable agreement between the two concentration distributions, which validates the numerical approach.  [16] for the dye mass fraction distribution at Re = 80. Table 4 shows the results of the single-objective optimization for the mixing effectiveness. The numerical results at the optimal design indicate an improvement of 58.9% in the mixing effectiveness compared to the reference design. This indicates an enhancement of 3.7% in the mixing index at the exit and a reduction of 34.8% in the pressure drop in comparison to the reference design as a result  [16] for the dye mass fraction distribution at Re = 80. Table 4 shows the results of the single-objective optimization for the mixing effectiveness. The numerical results at the optimal design indicate an improvement of 58.9% in the mixing effectiveness compared to the reference design. This indicates an enhancement of 3.7% in the mixing index at the exit and a reduction of 34.8% in the pressure drop in comparison to the reference design as a result of the optimization. The optimum design was obtained for w 3 /w 1 = 0.532 and h 1 /h 2 = 3.067, which are close to the middle range and upper bound, respectively. The surrogate model predicts the mixing effectiveness with a relative error of 2.5% in comparison to the numerical result with the optimal design.  Figure 6 shows a comparison of the mixing effectiveness in the axial direction for the reference and optimum micromixers. The mixing effectiveness was calculated on the y-z planes located after each mixing unit and at the exit (x/L t = 0.23, 0.41, 0.62, and 1). The mixing effectiveness is almost constant along the length of the optimized design, while it increases in the reference geometry. The optimum design shows much higher mixing effectiveness than the reference design over the entire length of the micromixer.   Figure 6 shows a comparison of the mixing effectiveness in the axial direction for the reference and optimum micromixers. The mixing effectiveness was calculated on the y-z planes located after each mixing unit and at the exit (x/Lt = 0.23, 0.41, 0.62, and 1). The mixing effectiveness is almost constant along the length of the optimized design, while it increases in the reference geometry. The optimum design shows much higher mixing effectiveness than the reference design over the entire length of the micromixer. Following the procedure in Figure 3, a multi-objective optimization was performed to find the Pareto-optimal solutions, as shown in Figure 7. The figure shows optimal compromises between the pair of conflicting objective functions: FMI and FΔP. The concave feature of the Pareto-optimal front signifies that a gain in the mixing index occurs at the expense of a higher pressure drop. Each design corresponding to a solution on the Pareto-optimal front is a universally best design. Therefore, Pareto-optimal solutions provide a choice for the selection of the preferred combination of objective functions. Following the procedure in Figure 3, a multi-objective optimization was performed to find the Pareto-optimal solutions, as shown in Figure 7. The figure shows optimal compromises between the pair of conflicting objective functions: F MI and F ∆P . The concave feature of the Pareto-optimal front signifies that a gain in the mixing index occurs at the expense of a higher pressure drop. Each design corresponding to a solution on the Pareto-optimal front is a universally best design. Therefore, Pareto-optimal solutions provide a choice for the selection of the preferred combination of objective functions.
Pareto-optimal solutions, as shown in Figure 7. The figure shows optimal compromises between the pair of conflicting objective functions: FMI and FΔP. The concave feature of the Pareto-optimal front signifies that a gain in the mixing index occurs at the expense of a higher pressure drop. Each design corresponding to a solution on the Pareto-optimal front is a universally best design. Therefore, Pareto-optimal solutions provide a choice for the selection of the preferred combination of objective functions. Five symbolic Pareto-optimal designs (PODs) were selected for analysis in terms of the two conflicting objective functions. These designs are marked on the Pareto-optimal curve in Figure 7 and listed in Table 5. POD 1 and POD 5 located at the extremes of the Pareto-optimal curve are mixing index-oriented and pressure drop-oriented designs, respectively. The fulfillment of one objective function causes deterioration in the competing objective function. The mixing index-oriented design is closer to the lower and upper limits of the design ranges for w 3 /w 1 and h 1 /h 2 , respectively. The objective function values predicted by the surrogate model are compared with the numerically calculated values for the representative PODs, as shown in Figure 8 and Table 5. The forecasted values of PODs show good agreement with the numerically calculated values for these PODs. The maximum relative errors in the surrogate prediction compared to the numerical results of the representative PODs were 2.7% and 5.9% for the mixing index at the exit and the pressure drop, respectively. The numerical result of POD 1 shows that a maximum improvement of 48.5% in the mixing index is achieved with POD 1 compared to the reference design. Similarly, a maximum reduction of 55.0% in pressure drop is achieved with POD 5 compared to the reference design.
The literature studied indicates that single-objective and multi-objective optimizations have been carried on different micromixers but not compared. However, in this study, the comparison of results between single-and multi-objective optimizations of the micromixer indicates that the mixing index improvements achieved by the multi-objective optimization are higher than the single-objective optimization. Also, the pressure drop achieved by the multi-objective optimization is lower than the single-objective optimization result. Hence, the multi-objective optimization is more suited to the optimization of the micromixer performance.
forecasted values of PODs show good agreement with the numerically calculated values for these PODs. The maximum relative errors in the surrogate prediction compared to the numerical results of the representative PODs were 2.7% and 5.9% for the mixing index at the exit and the pressure drop, respectively. The numerical result of POD 1 shows that a maximum improvement of 48.5% in the mixing index is achieved with POD 1 compared to the reference design. Similarly, a maximum reduction of 55.0% in pressure drop is achieved with POD 5 compared to the reference design. The literature studied indicates that single-objective and multi-objective optimizations have been carried on different micromixers but not compared. However, in this study, the comparison of results between single-and multi-objective optimizations of the micromixer indicates that the mixing index improvements achieved by the multi-objective optimization are higher than the singleobjective optimization. Also, the pressure drop achieved by the multi-objective optimization is lower The dye concentration distributions are plotted on transverse planes along the length of the micromixer (x/L t = 0, 0.23, 0.41, and 1) for a qualitative comparison of the mixing between POD 1 and POD 5, as shown in Figure 9. Greater distortion of the fluid interfaces is observed for POD 1 due to the enhanced stretching and folding process, which extends the interfacial contact area for diffusion and thereby establishes a more uniform concentration distribution at a much shorter axial length than POD 5. than the single-objective optimization result. Hence, the multi-objective optimization is more suited to the optimization of the micromixer performance.
The dye concentration distributions are plotted on transverse planes along the length of the micromixer (x/Lt = 0, 0.23, 0.41, and 1) for a qualitative comparison of the mixing between POD 1 and POD 5, as shown in Figure 9. Greater distortion of the fluid interfaces is observed for POD 1 due to the enhanced stretching and folding process, which extends the interfacial contact area for diffusion and thereby establishes a more uniform concentration distribution at a much shorter axial length than POD 5. The flow structure responsible for the dissimilar mixing qualities in POD 1 and POD 5 were investigated using the 3D streamlines of the fluids entering from the two inlets shown in Figure 10. The fluids from the two inlets come into contact at the T-joint. Distinct disparity in the streamlines is observed at the splitting of the flow. At the second joint, the flow from a single inlet splits into the major and minor sub-channels of the second mixing unit in POD 1, but most of the flow from the first major sub-channel passes through the second major sub-channel only in POD 5. The splitting of the flow into the two sub-channels in POD 1 causes increased stretching and folding of the fluid interface in the collision zone, as shown in Figure 10. The flipping of streamlines and flow circulation are not presented in POD 5, but are observed in the major sub-channels of POD 1, which accelerates the mixing. In addition, the smaller flow area at the dislocations of POD 1 than in POD 5 increases the local flow velocity, which also enhances the flow perturbation and thereby aids in mixing. The flow structure responsible for the dissimilar mixing qualities in POD 1 and POD 5 were investigated using the 3D streamlines of the fluids entering from the two inlets shown in Figure 10. The fluids from the two inlets come into contact at the T-joint. Distinct disparity in the streamlines is observed at the splitting of the flow. At the second joint, the flow from a single inlet splits into the major and minor sub-channels of the second mixing unit in POD 1, but most of the flow from the first major sub-channel passes through the second major sub-channel only in POD 5. The splitting of the flow into the two sub-channels in POD 1 causes increased stretching and folding of the fluid interface in the collision zone, as shown in Figure 10. The flipping of streamlines and flow circulation are not presented in POD 5, but are observed in the major sub-channels of POD 1, which accelerates the mixing. In addition, the smaller flow area at the dislocations of POD 1 than in POD 5 increases the local flow velocity, which also enhances the flow perturbation and thereby aids in mixing.
The fluids from the two inlets come into contact at the T-joint. Distinct disparity in the streamlines is observed at the splitting of the flow. At the second joint, the flow from a single inlet splits into the major and minor sub-channels of the second mixing unit in POD 1, but most of the flow from the first major sub-channel passes through the second major sub-channel only in POD 5. The splitting of the flow into the two sub-channels in POD 1 causes increased stretching and folding of the fluid interface in the collision zone, as shown in Figure 10. The flipping of streamlines and flow circulation are not presented in POD 5, but are observed in the major sub-channels of POD 1, which accelerates the mixing. In addition, the smaller flow area at the dislocations of POD 1 than in POD 5 increases the local flow velocity, which also enhances the flow perturbation and thereby aids in mixing.  The mixing capability of POD 1 is compared with that of the reference design at different Reynolds numbers, as shown in Figure 11. Both the reference micromixer and POD 1 achieve a mixing index over 0.90 at the exit for higher Reynolds numbers (Re = 60, 80, and 100). However, POD 1 achieves this mixing index at a much shorter length. Therefore, the mixing index at the end of the second mixing unit (x/Lt = 0.41) was chosen for comparison instead of the exit to clarify the improvement in the mixing. POD 1 shows better mixing capability than the reference design at all Reynolds numbers. This design shows improvements of 77.8% and 26.2% in the mixing index at Reynolds numbers of 30 and 40, respectively. The mixing capability of POD 1 is compared with that of the reference design at different Reynolds numbers, as shown in Figure 11. Both the reference micromixer and POD 1 achieve a mixing index over 0.90 at the exit for higher Reynolds numbers (Re = 60, 80, and 100). However, POD 1 achieves this mixing index at a much shorter length. Therefore, the mixing index at the end of the second mixing unit (x/Lt = 0.41) was chosen for comparison instead of the exit to clarify the improvement in the mixing. POD 1 shows better mixing capability than the reference design at all Reynolds numbers. This design shows improvements of 77.8% and 26.2% in the mixing index at Reynolds numbers of 30 and 40, respectively. The mixing effectiveness (ME) [41], mixing cost (MC) [42], and mixing energy cost (MEC) [43,44] have been used as performance parameters in the literature to evaluate micromixers. Higher mixing effectiveness and MC and a lower value MEC represent better micromixer performance. The expressions for the MC and MEC can be written as follows: where ΔP, , and η denote pressure drop, mean input power coefficient, and mixing efficiency, respectively. Mixing efficiency and mean input power coefficient ( ) are written as follows: Figure 11. Comparison of mixing index at different Reynolds numbers between POD 1 and reference design.
The mixing effectiveness (ME) [41], mixing cost (MC) [42], and mixing energy cost (MEC) [43,44] have been used as performance parameters in the literature to evaluate micromixers. Higher mixing effectiveness and MC and a lower value MEC represent better micromixer performance. The expressions for the MC and MEC can be written as follows: where ∆P, C p , and η denote pressure drop, mean input power coefficient, and mixing efficiency, respectively. Mixing efficiency and mean input power coefficient (C p ) are written as follows: where σ, σ max , ρ, ∆P, q, and V are the standard deviation of the concentration, maximum standard deviation, density of fluid, pressure drop, dimensionless flow rate, and average flow velocity, respectively. The performance parameters of the representative PODs were compared with those of three previous micromixers involving unbalanced collisions, as shown in Table 6. All the PODs outperform the previous micromixers in terms of the mixing index at the exit, but POD 1 shows ME and MC values that are less than or equal to all micromixers except that proposed by Xia et al. [45] due to higher pressure drop. In the case of MEC, all the PODs except POD 1 and POD 5 show lower values than those of the previous micromixers. POD 4 shows the best values for ME, MC, and MEC among all the micromixer designs. This micromixer design shows 315%, 37%, and 38% higher values of the mixing index, ME, and MC than the lowest values shown by the other micromixer designs, respectively.

Conclusions
Single and multi-objective optimizations of a 3D-ASAR micromixer were performed at Reynolds number of 20, using surrogate modeling in conjunction with mixing and flow analyses. The optimizations were carried out with two dimensionless design variables, w 3 /w 1 and h 1 /h 2 . The mixing effectiveness was selected as a single-objective function, while both the mixing index at the exit and pressure drop across the micromixer were chosen for the two-objective optimization. The KRG model was selected to generate the surrogate of the objective functions.
The parametric study demonstrated that the mixing index at the outlet shows non-monotonic behavior and has maxima, but the pressure drop shows a monotonically decreasing behavior for the two design variables. An improvement of 58.9% in the mixing effectiveness was obtained in comparison to the reference design by the single-objective optimization. This enhancement corresponds to a 3.7% rise in the mixing index at the exit and a 34.8% decrease in the pressure drop in comparison to the reference design. The surrogate prediction error for the single-objective optimization was estimated to be 2.5%.
The reciprocity between the two objective functions represented by the Pareto-optimal front for the multi-objective optimization was obtained by employing surrogate modeling and MOGA. The numerical analysis of five representative PODs indicated a maximum improvement of 48.5% in the mixing index for a mixing index-oriented design (POD 1) and a maximum reduction of 55.0% in the pressure loss for a pressure drop-oriented design (POD 5) in comparison to the reference design. On the other hand, POD 4 showed the best values for ME, MC, and MEC among all the micromixer designs, including three previous micromixers involving unbalanced collisions. The surrogate model predicted the mixing index and pressure drop with maximum relative errors of 2.7% and 5.9%, respectively. There were low percentages of improvement in the mixing index and the reduction in the pressure drop by the single-objective optimization in comparison to the multi-objective optimization, which indicates that the multi-objective optimization is much more favorable for improvement of the micromixer performance.