Next Article in Journal
Fabrication of Through via Holes in Ultra-Thin Fused Silica Wafers for Microwave and Millimeter-Wave Applications
Next Article in Special Issue
Numerical and Experimental Study on Mixing Performances of Simple and Vortex Micro T-Mixers
Previous Article in Journal
Optofluidics Refractometers
Previous Article in Special Issue
Multi-Objective Optimizations of a Serpentine Micromixer with Crossing Channels at Low and High Reynolds Numbers
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Topology Optimization of Passive Micromixers Based on Lagrangian Mapping Method

1
Changchun Institute of Optics, Fine Mechanics and Physics (CIOMP), Chinese Academy of Science, Changchun 130033, China
2
University of Chinese Academy of Science, Beijing 100049, China
3
State Key Laboratory of Mechanical System and Vibration, Shanghai Jiao Tong University, Shanghai 200240, China
4
State Key Laboratory of Applied Optics, Changchun Institute of Optics, Fine Mechanics and Physics (CIOMP), Chinese Academy of Science, Changchun 130033, China
*
Author to whom correspondence should be addressed.
Micromachines 2018, 9(3), 137; https://doi.org/10.3390/mi9030137
Submission received: 26 January 2018 / Revised: 17 March 2018 / Accepted: 18 March 2018 / Published: 20 March 2018
(This article belongs to the Special Issue Passive Micromixers)

Abstract

:
This paper presents an optimization-based design method of passive micromixers for immiscible fluids, which means that the Peclet number infinitely large. Based on topology optimization method, an optimization model is constructed to find the optimal layout of the passive micromixers. Being different from the topology optimization methods with Eulerian description of the convection-diffusion dynamics, this proposed method considers the extreme case, where the mixing is dominated completely by the convection with negligible diffusion. In this method, the mixing dynamics is modeled by the mapping method, a Lagrangian description that can deal with the case with convection-dominance. Several numerical examples have been presented to demonstrate the validity of the proposed method.

1. Introduction

Lab-on-a-chip devices have been widely used in the area of the analysis, synthesis, and separations due to the advantages of high efficiency, portability, and low reagent consumption [1]. Injection, mixing, reaction, cleaning, separation, and detection, which are the functions of conventional analytical laboratory, can be achieved on a centimeter-level chip [2]. Various microfluidic devices have been integrated in lab-on-a-chip, such as micropumps, microvalves, micromixers, and microchannel. Rapid and complete mixing can influence the efficiency of a microfluidic system [3]. Therefore, micromixer plays a significant role in the lab-on-a-chip devices. Based on actuation methods, micromixers can be classified into two categories: active micromixers and passive micromixers [4]. Active mixers use external energy to create chaotic convection, such as pressure [5], magnetohydrodynamics [6], electrokinetics [7] and acoustics [8,9,10]. Active mixers have short mixing times and distances and can be controlled to be on and off, according to the needs of the users. Because of the requirement of external energy, however, the fabrication and integration of active mixers is complicated and expensive [11]. Comparatively, passive mixers can achieve fluid mixing solely by the geometries of channels and be integrated in a complex microfluidic system simply and directly. The inconsistent cross sections of a microchannel can cause the fluid to be stretched and folded in the transversal direction [12,13,14,15]. A reasonable layout of passive micromixer can strengthen the chaotic convection to enhance the mixing performance. Immiscible fluids are widely used in chemical industry, where the particle flow can be considered as the immiscible fluids [16]. The mixture of immiscible fluids appears in biochemical experiments, which is usually not discussed compare to the mixers with phenomenon of convection-diffusion. With the trend of miniaturization in recent years, the micromixers of immiscible fluids with efficient mixing performance are also desired in a microfluidic system.
Topology optimization of fluid flows has been proposed by Borrvall and Petersson for Stokes flow [17] and the Navier-Stokes flow [18,19,20,21]. When compared to shape optimization method, the detail topology and shape of the microchannels can be obtained simultaneously by topology optimization. This method has been used to design the microchannel networks, micropumps, no-moving part microvalves, and micromixers [22,23,24,25]. In the topology optimization of micromixers, the convection-diffusion equation is usually used to describe the mixing process of the fluids, and objective function is the variance between the actual obtained concentration at the outlet and the expected concentration. All of the topological optimization methods that are mentioned above used to design the mixers of miscible fluids are based on the Eulerian description of the convection-diffusion dynamics and can not be directly applied to the design of immiscible fluid mixers. The mixture of immiscible fluids is a convective problem in physics and should be described by the Navier–Stokes (NS) equations and convection equation in numerical calculation. However, when using the standard Galerkin finite element method to solve the above equation, the numerical instability will occur. For an optimization process, the numerical instability of solving forward problems is a big challenge. To avoid this problem, the mapping method is used in this paper to describe the convection problem in numerical calculation. The mapping method proposed by Singh et al. [26,27] can describe the mixing performance by calculating the mapping matrix and be integrated into the topology optimization method. For the topology optimization problem, the measure of the mixing in the mapping method needs to be discussed. The coarse grained concentration [28] can change the value of the objective function by replacing the concentration in the objective function with the coarse grained concentration defined in discrete areas. In order to shorten the optimization time, the mixing performance of entire periodic structure is obtained by analyzing only one cycle structure in this paper.
In passive mixers, the mixing efficiency is mainly determined by the layout of mixers. The chaotic convection can be promoted by the various cross-sectional structures along the flow direction. Based on the topology optimization method of fluidic flows and Lagrangian description, this paper is focused on the layout design method of the passive micromixers of immiscible fluid, where the mapping method is used to describe the mixing process of the immiscible fluid. A new mixing measurement is applied in this paper, based on the change of position of the mixed fluids. This paper is organized as follows: the choice of mixing measurement and the mapping method are stated in Section 2; the topology optimization model of the passive micromixers and the corresponding adjoint equations and derivative are derived in Section 3; several numerical results are presented in Section 4; and, the discussion and conclusion are stated in Section 5.

2. Measure of Mixing and Lagrangian Mapping Description

2.1. Measure of Mixing

In the topology optimization method of micromixers based on Eulerian description of the convection-diffusion dynamics, we define the concentration as the volume fraction of a fluid in a mixture of two fluids at a point. To quantify and compare mixing performance, the least squares variance between the actual obtained concentration at the outlet and the expected concentration can be expressed as [29]:
J ( c ) = Γ o u t ( c c ¯ ) 2 d Γ Γ i n ( c 0 c ¯ ) 2 d Γ
where Гin and Гout are the inlet and outlet of the mixer, respectively; c0 is the reference concentration, which is usually the designer specified concentration at the inlet. Ideally, the two fluids are sufficiently mixed in the mixer. Therefore, the expected concentration c ¯ is chosen as the ideal concentration after sufficient mixing at the outlet, which is the average concentration. The concentrations of two fluids are set to be dimensionless value 0 and 1. When the two fluids are miscible, J(c) = 0 means that the two fluids can be considered to be completely mixed and J(c) = 1 means that the two fluids can be considered to be completely separated. However, when the two liquids are immiscible, the concentration c of a fluid is always be 0 or 1, due to the absence of diffusion and J(c) remains constant. The expression in Equation (1) is not suitable to quantify the mixing performance. The coarse grain concentration method that was proposed by Welander et al. [28] can avoid this situation. The coarse grain concentration Ci defined on a finite cell Гi is:
C i = Γ i c ( X ) d Γ
where Гi is the i-th cell obtained discretely on the cross section perpendicular to the flow direction; Ci can vary in interval [0,1]. Therefore, Equation (1) can be rewrite as
J C ¯ ( C ) = 1 A Γ i = 1 N ( C i C ¯ ) 2 ( C i 0 C ¯ ) 2 A Γ i , C ¯ = 1 A Γ i = 1 N C i A Γ i
where A Γ and A Γ i is the area of cross section and the i-th cell, respectively; C ¯ is the average coarse grain concentration; C i 0 is the coarse grain concentration of the i-th cell on the cross section of inlet; and N is the number of cells divided in cross section.
Since the mixed fluids are immiscible, the coarse grain concentration changes only at the contact cells of the two fluids. Due to the low Reynolds numbers of fluid and the limited mixing length of micromixer used in the topology optimization model, the change of the value of mixing performance in Equation (3) is not large enough. When the number of discrete cells increases at the same time, the value of Equation (3) is not sensitive to the change of coarse grain concentration. A new measurement mixing is proposed to amplify the value modification in the mixing performance by using the change of the positions of two fluids at the cross section of outlet. The least squares variance between the actual obtained coarse grain concentration Ci at the outlet and the initial coarse grain concentration C i 0 at the inlet is used:
J C 0 ( C ) = 1 1 + 1 A Γ i = 1 N ( C i C i 0 ) 2 A Γ i
When J C 0 ( C ) = 1 , the two immiscible fluids can be considered to be completely separated; when J C 0 ( C ) < 1 , the two immiscible fluids can be considered to be mixing.

2.2. Mapping Method

A distribution matrix is used in the mapping method to store the information, which describes the changes of fluidic distribution between two specified cross sections [30,31]. To obtain the each coefficient of mapping matrix, the initial cross section is divided into a large number of discrete cells with specified size. The material of fluid transferred to several recipient cells from a donor cell along the fluidic flow. The fraction of material in the recipient cell Ωj in section at X = X0 + ΔX, which is found in the donor cell Ωi in section at X = X0 is the coefficient of mapping matrix. The cross section with N cells can construct a distribution matrix of the order N × N:
φ i j = Ω j | x = x 0 + Δ x Ω i | x = x 0 d A Ω j | x = x 0 d A
To describe the detail of flow, a lot of sections are traced which are set along the flow direction. Tracing all cells in all sections during a flow over a distance ΔX, the detail of flow can be described by the complex deformation of cells. Although this tracking method is feasible, it is too time-consuming to apply into topology optimization. Singh et al. proposed a convenient method for calculating the coefficients of the mapping matrix [26,27]. To approximate the coefficients of mapping matrix, K markers are filled into each cell uniformly in donor section at X = X0, and then, tracing these markers can obtain the information about the distribution in recipient section at X = X0 + ΔX. When the number of markers in the cell Ωj in donor section is Mj and Mij markers are traced in the cell Ωi in recipient section, the coefficient of mapping matrix can be calculated as
Φ i j = M i j M j
The convection dynamics of fluids can be analyzed the following procedure. Using the same cell to describe the coarse grain concentration, the coarse grain concentration distribution of cross-section can construct a vector C R N × 1 (N is the number of cells). After passing through the structure that is described by the mapping method, the coarse grain concentration distribution of recipient section Cr can be calculated from the coarse grain concentration distribution of donor section Cd as:
C r = Φ C d
The coarse grain concentration distribution at outlet can be obtained as:
C o = Φ a l l C i
where Ci and Co are the coarse grain concentration distribution of inlet and outlet cross-section, respectively; and Φall is the mapping matrix of whole mixer.
Since the construction of the mapping matrix has no correlation with the initial cross-section concentration distribution, the mixer of periodic layout, with a known initial mapping matrix Φ1 of a single cycle, has the following relation:
C i + 1 = Φ 1 C i , C n = ( Φ 1 ( Φ 1 ( ( Φ 1 n   times C 0 ) ) ) )
where Ci is the coarse grain concentration distribution after the i-th mixing. Therefore, for the periodic mixers in any cycle, the computational cost is saved that the concentration distribution. Cn can be obtained by simply multiplying the single cycle matrix in corresponding times.
Assuming that the length direction of mixing channel is the x-axis, the cross section that is shown in Figure 1 is the y-z section. The trajectories of markers are tracked based on the coordinate axis. Tracing can be realized by the axial velocity component ux and the transversal velocity components uy and uz, respectively:
d y d x = u y u x , d z d x = u z u x
Since the particles become disordered at any downstream position in forward tracing method, there is no guarantee that an equal number of traced markers can be found in each cell at the outlet. Using the backward particle tracing (BTP) to construct the distribution matrix, the markers initially fill the recipient cell in the cross section of outlet, and they are traced backward against the flow direction [26,27],
X i = X 0 + i = 1 n ( x i + Δ x x i u u T I u d X )
where Xi = (xi, yi, zi) and X0 = (x0, y0, z0) are the coordinates of the markers at the cross section of inlet and outlet, respectively; I u = [ 1   0   0 ] T ; n = nall − 1, nall is the number of cross section in the mapping method. Due to the integration of axial spatial increment rather than time, the error that is caused by the different time distribution can be eliminated. However, this approach is not effective when the fluid is reflowing.
Figure 2 shows a top view of the grooves on the bottom in a staggered herringbone mixer (SHM). We apply the geometry of SHM and the material properties used in the study of Singh et al. [26]. The length ratio of two arms of every groove is 2:1, and all arms are at 45° to the axial direction; the channel height h = 77 μm, the channel width w = 200 μm, the depth of grooves gd = 17.7 μm, the width of grooves gw = 70.7 μm and the distance between two grooves also equals 70.7 μm; the viscosity and density of the fluid in the micromixers are 0.067 kg·m/s and 1.2 × 103 kg/m3. The average inlet velocity u = 0.2 cm/s. The velocity field is obtained by 144,000 hexahedral elements and 155,031 nodes. Figure 3a–c show the mixing evolutions in a SHM consisting of one groove type, whose mapping matrix is represented by Φ1. In the mapping computations of Figure 3a,b, the tracing cross section is covered with 50 × 50 and 100 × 100 cells, respectively, and each cell is filled with 100 uniformly distributed markers. Figure 3c shows the results in the study of Singh et al. [26] using the above parameters. Figure 3d shows the coarse grain concentration distribution at outlet of SHM with ten grooves. The results in Figure 3 demonstrated the validity of the mapping method that is used in this paper. These parameters of mapping method will be used in the optimization process.

3. Topology Optimization Model of Mixers

When the area of each cell divided by the mapping method is the same, the Equation (4) can be simplified to
J C 0 ( C ) = N N + i = 1 N ( C i C i 0 ) 2
where N is the number of cells divided in cross section. This discrete equation is used as quantitative criterion to measure the mixing performance of the immiscible fluids. In this paper, the key point is how to find a reasonable layout of a micromixer in a designer specified design domain, which minimizes the quantitative criterion and mixing can be obtained using the topology optimization method. Based on the continuity assumption, the fluidic field will be described by the Navier-Stokes equations:
ρ ( u · ) u η Δ u + p = f · u = 0
where u, p, ρ, and η are the velocity, pressure, density, and viscosity of the fluid, respectively; f is the body forces acting on the fluid. In the topology optimization method, an artificial friction force is introduced into the Navier-Stokes equations, which was proposed for the Stokes flow by Borrvall and Petersson [17] and generalized to the Navier-Stokes flow by Gersborg et al. [20] and Olesen et al. [21]. Initially, an artificial porous material is uniformly distributed in the design domain; and then, the artificial porous material forms solid and liquid phases gradually; at last, the high and low impermeability characterize the solid phase and the fluid phase, respectively. The artificial friction force is f = −α(γ)u, and α is the impermeability of the artificial porous material, and γ is the design variable. The design variable γ varies in interval [0,1], where 0 and 1 denote the solid and fluid phases, respectively. The impermeability of porous material α is the interpolation function of design variable γ [17,18,19,20,21,22]:
α ( γ ) = α min + ( α max α min ) q ( 1 γ ) ( q + γ )
where αmax and αmin are the impermeability of the solid phase and the fluid phase, respectively; q is a positive value used to adjust the convexity of the interpolation function; αmin is chosen 0 in the fluidic topology optimization. To obtain the perfect impermeability of solid no-slip boundary, αmax should be infinite, but a finite number has to be chosen to ensure the numerical stability. The layout of the passive micromixer of immiscible fluid can be determined by seeking the distribution of the design variable γ. The topology optimization model based on mapping method is:
min : J C 0 ( C ) = N N + i = 1 N ( C i C i 0 ) 2 s . t . ρ ( u · ) u η Δ u + p = α u    in   Ω · u = 0    in   Ω u = u 0    on   Γ i n [ p I + η ( u + ( u ) T ) ] · n = 0    on   Γ o u t 0 γ 1 X = X 0 + i = 0 n ( x i + Δ x x i u u T I u d X )    on   Γ i n X = X 0    on   Γ o u t C = C 0    on   Γ i n C = Φ ( y i , z i ) C 0    on   Γ o u t Φ i j = M i j M j
where Ω = Ω D Ω C is the computational domain; ΩD is the design domain; ΩC is the channels connected to the inlet, outlet, and design domain.
The constraint optimization model in Equation (15) can be transferred into unconstrained one by Lagrangian multiplier method:
L = J + a ( u , λ u ) Ω + ρ ( ( u · ) u , λ u ) Ω + ( p , λ u ) Ω + ( α u , λ u ) Ω ( g , λ u ) Γ N + ( u u 0 , λ u ) Γ D + ( · u , λ p ) Ω + ( C Φ ( y , z ) C 0 , λ C ) Γ N + ( C C 0 , λ C ) Γ D + ( X X 0 , λ X ) Γ D + ( X i X 0 i = 0 n ( x i + Δ x x i u u T I u d X ) , λ X ) Γ N
where a ( u , λ u ) Ω = Ω η u : λ u d Ω ; ( , ) Ω and ( , ) Γ is the inner product on the computational domain and boundary; λu, λp, λC and λX are the adjoint variable of the velocity field, the pressure field, the coarse grain density on the cross outlet section and the coordinate after particle tracing, respectively.
The variation of the Equation (16) is
δ L = L u δ u + L p δ p + L X δ X + L C δ C + L γ δ γ
According to the Karush-Kuhn-Tucker conditions for the partial differential equation constrained optimization problem, the optimization problem in Equation (15) can be obtained by solving the following adjoint equations [19] (see the Appendix A for the detailed derivation):
{ η Δ λ u ρ ( u · ) λ u + ρ ( u ) · λ u + λ p = α λ u in   Ω · λ u = 0 in   Ω λ u = 0 on   Γ D [ η ( λ u + ( λ u ) T ) λ p I ] · n = ρ ( u · n ) λ u + i = 0 n ( x 0 + Δ x x 0 ( 1 u T I u I u ( u T I u ) 2 I u T ) d X ) λ X on   Γ N λ X = Φ X i λ C on   Γ N λ X = 0 on   Γ D λ C = 0 on   Γ D λ C = 2 N ( C C 0 ) [ N + i = 1 N ( C i C i 0 ) 2 ] 2 on   Γ N
The adjoint sensitivity of the optimization model in Equation (15) is:
D L D γ | Ω = α γ u · λ u
When considering the manufacturability, the micromixer is designed to be a single layer structure, like the SHM [12,13,14,32,33]. In the optimization model, an additional constraint is added to ensure that the design variable in the depth direction have a consistent value. The design variable is defined on the plane x0y (Figure 4), and the adjoint derivative of the optimization model (19) is changed to [22]:
D L D γ = 1 z t z b z b z t α γ u · λ u d z
where zt and zb are the coordinates in z direction of top and bottom surface of the design domain, and zt > zb.
The topology optimization problem is solved by using the gradient-based iterative approach. In the optimization iterations, the Navier-Stokes equations, the backward particles tracing equation, and the corresponding adjoint equations in the weak form are solved by the finite element method using the commercial software COMSOL Multiphysics (version 3.5, COMSOL Inc., Stockholm, Sweden). The method of moving asymptotes (MMA) is used to update the design variable. The optimization iterations are stopped, when the maximal change of the objective value in three consecutive iterations less than 1 × 10−3. The procedure of solving the topology optimization problem for the layout design of passive micromixers is listed in below.
  • Give the initial value of the design variable γ;
  • Solve the Navier-Stokes equations and backward particle tracing equation by the finite element method;
  • Solve the adjoint equation;
  • Compute the adjoint derivative and the corresponding objective and constraint values;
  • Update the design variable by method of moving asymptotes (MMA);
  • Check for convergence; if the stopping conditions are not satisfied, go to 2; and
  • Post-processing

4. Results and Discussion

To demonstrate the validity of the proposed method, passive micromixers of immiscible fluid in a straight microchannel with external driven flow (constant flow rate at inlet) is investigated in the following. The design domain is shown in Figure 4a, and each cubic space with an edge size equals to L is discretized by 20 × 20 × 20 hexahedral elements. The width of micromixer L is 400 μm and the height of design domain H is 240 μm. The viscosity and density of the fluid in the micromixers are 1 × 10−3 Pa·s and 1 × 103 kg/m3. The value of αmax and q in the topology method are chosen 1 × 107 and 0.1. The initial value of the design variable γ is 0.4, which should be between [0,1], as shown in Equation (15). Since the backward particle tracing method is applied, the number of cell in the cross section of outlet is 50 × 50 and 100 markers are traced in each cell. The coarse grain concentrations of two immiscible fluids are set to be dimensionless value 1 and 0 (Figure 4b), respectively. All of the computations are performed on a DELL workstation (DELL Optiplex 7040, Intel Core i7-6700 CPU, 16 gigabyte memory, Round Rock, TX, USA).
Based on the topology optimization model in Section 3, the optimal layouts of the passive micromixers of immiscible fluid with different Reynolds number and length of the design domain are obtained, as shown in Figure 5, Figure 6 and Figure 7. The optimized results strongly depend on the parameters, such as the size of design domain and the Re number. From the comparison between the results in Figure 5 and Figure 6, the obtained layouts of micromixers depend on the selection of Reynolds number. Different values of n are chosen to compare the mixing performance of the micromixers obtained by the topology optimization method (Figure 6 and Figure 7). A larger value of n means a longer mixing length and a more reasonable layout in the micromixers. When the length of the mixer is relatively short, the herringbone type structure is not necessary; when the length of single-periodic becomes longer, the herringbone type structure promotes the mixing effect obviously. By comparing the optimal results, the herringbone type structure is a reasonable topology, which could be used further for detailed shape or size optimization when the mixer has an appropriate length.
Mapping method provides an approximate way to obtain the mixing performance of mixer with spatial periodic layout by multiplying the single cycle mapping matrix in corresponding times. Figure 8 shows the evolution of mixing performance of micromixers with different cycles where the single cycle is the layout in Figure 5. Due to the Equation (4) is only valid for situations where fluid agitation is not obvious, we used Equation (3) to measure the mixing performance. One can see that the mixing effect is enhanced as the number of cycle increases. The effect of numerical diffusion becomes apparent, when the multiplying times of mapping matrix increase. However, a mutual comparison is still possible and valuable.
The CPU time in one optimization iteration step is 500 s for the case that design domain is shown in Figure 4, with n = 2 and Re = 5. Therefore, designing a multi-period structure directly means having a longer computational domain and more CPU time. In contrast, the mapping method can obtain multi-period mixing performance easily by simply multiplying the single cycle matrix in corresponding times (Figure 8). Figure 9 shows the mixing performance of the micromixer that is shown in Figure 5 and SHM shown in Figure 10, which has only one groove and same volume as the micromixer shown in Figure 5. Therefore, there are full of room to adjust the expression of objective in the optimization model (Equation (15)) proposed in this paper. Whereas, the derivation of adjoint sensitivity is unchanged in the case that the chosen objective function is differentiable for the design variable. The whole optimization procedure is valid for immiscible mixer design, as illustrated in this paper.

5. Conclusions

In this paper, a novel method is used to design the layout of the passive micromixers of immiscible fluid has been proposed based on the topology optimization of fluidic flows. The layout of the passive micromixers is determined by solving a topology optimization problem to minimize the mixing measurement. Additionally, the detailed mixing performance is obtained by using mapping method in Lagrangian description without consideration of the diffusion. The mapping method used in this paper is a rough approximated method. However, it can produce results that are in good agreement with the conventional numerical simulation method of up to 10 cycles. The discrepancies may increase with the using of more cycles. The design variable is set to represent the impermeability distribution of the artificial porous medium in the design domain. Based on the adjoint analysis of the topology optimization problem, the design variable is evolved to derive impermeability distribution of the artificial porous medium with low and high levels, which, respectively, correspond to the solid and fluid phases in the micromixer. The manufacturability of the obtained layout is ensured by the manufacturing constraint. Numerical results demonstrated the validity of the proposed method. In addition, this method can be extended to design micromixers with multi-layer structure and active micromixers. These will be investigated in our future work.

Acknowledgments

This research was funded by the National Science Foundation of China under grant Nos. 51675506, 51405465. The authors are grateful to Krister Svanberg for supply of the Matlab codes of MMA.

Author Contributions

Yifan Xu, Yuchen Guo, Yongbo Deng and Zhenyu Liu wrote the topology optimization code; Yongbo Deng, Yifan Xu, Zhenyu Liu and Yuchen Guo derived the theoretical formulas; Yuchen Guo and Zhenyu Liu performed the numerical examples; Yuchen Guo, Zhenyu Liu and Yongbo Deng wrote the paper.

Conflicts of Interest

The authors declare no conflict of interest.

Appendix A

Transfer the constraint model in Equation (15) into unconstrained format by Lagrangian multiplier method:
L = J + a ( u , λ u ) Ω + ρ ( ( u · ) u , λ u ) Ω + ( p , λ u ) Ω + ( α u , λ u ) Ω ( g , λ u ) Γ N + ( u u 0 , λ u ) Γ D + ( · u , λ p ) Ω + ( C Φ ( y , z ) C 0 , λ C ) Γ N + ( C C 0 , λ C ) Γ D + ( X X 0 , λ X ) Γ D + ( X i X 0 i = 0 n ( x i + Δ x x i u u T I u d X ) , λ X ) Γ N
where a ( u , λ u ) Ω = Ω η u : λ u d Ω ; ( , ) Ω and ( , ) Γ is the inner product on the computational domain and boundary; λu, λp, λC, and λX are the adjoint variables of the velocity field, the pressure field, the coarse grain density on the cross outlet section and the coordinate after particle tracing, respectively. λu and λC are vectors as well as λp is scalar. λX = (λx λy λz)T. The variation of the main function Equation (A1) is:
δ L = L u δ u + L p δ p + L X i δ X i + L C δ C + L γ δ γ
According to the Karush-Kuhn-Tucker conditions for the partial differential equation constrained optimization problem, the optimization model in Equation (15) can be solved by solving the following adjoint equations:
L u δ u = 0
L p δ p = 0
L X i δ X i = 0
L X i δ X i = 0
For the Equation (A3):
L u δ u = a ( η δ u , λ u ) Ω + ρ δ ( ( u · ) u , λ u ) Ω + ( α δ u , λ u ) Ω ( · δ u , λ p ) Ω + ( δ u , λ u ) Γ D ( x 0 + Δ x x 0 ( 1 u T I u I u ( u T I u ) 2 I u T ) d x , λ X ) Γ N
where:
a ( η δ u , λ u ) Ω = ( η Δ λ u , δ u ) Ω + ( η ( λ u + ( λ u ) T ) · n , δ u ) Γ N
δ ( ( u · ) u , λ u ) Ω = ( ( u · ) λ u , δ u ) Ω + ( u · λ u , δ u ) Ω + ( ( u · n ) λ u , δ u ) Γ N
( · δ u , λ p ) Ω = ( λ p , δ u ) Ω + ( λ p n , δ u ) Γ N
By inserting Equations (A8)–(A10) into Equation (A7), one can obtain:
L u δ u = Ω [ η Δ λ u ρ ( u · ) λ u + ρ ( u ) · λ u + λ p + α λ u ] · δ u d Ω + Γ D [ λ u ] · δ u d Γ + Γ N [ ( η ( λ u + ( λ u ) T λ p I ) · n + ρ ( u · n ) λ u i = 0 n ( x 0 + Δ x x 0 ( 1 u T I u I u ( u T I u ) 2 I u T ) d X ) λ X ] · δ u d Γ = 0
For the Equation (A4):
L p δ p = Ω [ · λ u ] · δ p d Ω = 0
For the Equation (A5):
L X i δ X i = Γ N [ λ X Φ X i λ C ] · δ X i d Γ + Γ D [ λ X ] · δ X i d Γ = 0
For the Equation (A6):
L C δ C = Γ N [ 2 N ( C C 0 ) [ N + i = 1 N ( C i C i 0 ) 2 ] 2 + λ C ] · δ C d Γ + Γ D [ λ C ] · δ C d Γ = 0
Finally, the adjoint equation can be obtained as:
{ η Δ λ u ρ ( u · ) λ u + ρ ( u ) · λ u + λ p = α λ u in   Ω · λ u = 0 in   Ω λ u = 0 on   Γ D [ η ( λ u + ( λ u ) T ) λ p I ] · n = ρ ( u · n ) λ u + i = 0 n ( x 0 + Δ x x 0 ( 1 u T I u I u ( u T I u ) 2 I u T ) d X ) λ X on   Γ N λ X = Φ X i λ C on   Γ N λ X = 0 on   Γ D λ C = 0 on   Γ D λ c = 2 N ( C C 0 ) [ N + i = 1 N ( C i C i 0 ) 2 ] 2 on   Γ N
According to the Equations (A3)–(A6) into the variation Equation (A2) of main function, one can obtain:
δ L = L γ δ γ = Ω [ α γ u · λ u ] δ γ d Ω
Then the adjoint derivative of the optimization model in Equation (15) is:
D L D γ | Ω = α γ u · λ u

References

  1. Whitesides, G. The lab finally comes to the chip! Lab Chip 2014, 14, 3125–3126. [Google Scholar] [CrossRef] [PubMed]
  2. Manz, A.; Graber, N.; Widmer, H.M. Miniaturized total chemical-analysis systems—A novel concept for chemical sensing. Sens. Actuator B-Chem. 1990, 1, 244–248. [Google Scholar] [CrossRef]
  3. Lee, C.-Y.; Chang, C.-L.; Wang, Y.-N.; Fu, L.-M. Microfluidic mixing: A review. Int. J. Mol. Sci. 2011, 12, 3263–3287. [Google Scholar] [CrossRef] [PubMed]
  4. Suh, Y.K.; Kang, S. A review on mixing in microfluidics. Micromachines 2010, 1, 82–111. [Google Scholar] [CrossRef]
  5. Glasgow, I.; Aubry, N. Enhancement of microfluidic mixing using time pulsing. Lab Chip 2003, 3, 114–120. [Google Scholar] [CrossRef] [PubMed]
  6. Qian, S.; Zhu, J.; Bau, H.H. A stirrer for magnetohydrodynamically controlled minute fluidic networks. Phys. Fluids 2002, 14, 3584–3592. [Google Scholar] [CrossRef]
  7. Chen, L.; Deng, Y.; Zhou, T.; Pan, H.; Liu, Z. A novel elecrtoosmotic micromixer with asymmetric lateral structures and DC electrode arrays. Micromachines 2017, 8, 105. [Google Scholar] [CrossRef]
  8. Ozcelik, A.; Ahmed, D.; Xie, Y.; Nama, N.; Qu, Z.; Nawaz, A.A.; Huang, T.J. An acoustofluidic micromixer via bubble inception and cavitation from microchannel sidewalls. Anal. Chem. 2014, 86, 5083–5088. [Google Scholar] [CrossRef] [PubMed]
  9. Huang, P.H.; Xie, Y.; Ahmed, D.; Rufo, J.; Nama, N.; Chen, Y.; Chan, C.Y.; Huang, T.J. An acoustofluidic micromixer based on oscillating sidewall sharp-edges. Lab Chip 2013, 13, 3847–3852. [Google Scholar] [CrossRef] [PubMed]
  10. Nama, N.; Huang, P.H.; Huang, T.J.; Costanzo, F. Investigation of acoustic streaming patterns around oscillating sharp edges. Lab Chip 2014, 14, 2824–2836. [Google Scholar] [CrossRef] [PubMed]
  11. Nguyen, N.-T.; Wu, Z. Micromixers—A review. J. Micromech. Microeng. 2005, 15, R1–R16. [Google Scholar] [CrossRef]
  12. Stroock, A.D.; Dertinger, S.W.; Ajdari, A.; Mezić, I.; Stone, A.; Whitesides, G.M. Chaotic mixer for microchannels. Science 2002, 295, 647–651. [Google Scholar] [CrossRef] [PubMed]
  13. Liu, Y.; Kim, B.J.; Sung, H.J. Two-fluid mixing in a microchannel. Int. J. Heat Fluid Flow 2004, 25, 986–995. [Google Scholar] [CrossRef]
  14. Camesasca, M.; Kaufman, M.; Manaszloczower, I. Staggered passive micromixers with fractal surface patterning. J. Micromech. Microeng. 2006, 16, 2298–2311. [Google Scholar] [CrossRef]
  15. Lee, S.W.; Kim, D.S.; Lee, S.S.; Kwon, T.H. A split and recombination micromixer fabricated in a PDMS three-dimensional structure. J. Micromech. Microeng. 2006, 16, 1067–1072. [Google Scholar] [CrossRef]
  16. Meijer, H.E.H.; Singh, M.K.; Anderson, P.D. On the performance of static mixers: A quantitative comparison. Prog. Polym. Sci. 2012, 37, 1333–1349. [Google Scholar] [CrossRef]
  17. Borrvall, T.; Petersson, J. Topology optimization of fluids in stokes flow. Int. J. Numer. Methods Fluids 2003, 41, 77–107. [Google Scholar] [CrossRef]
  18. Deng, Y.; Liu, Z.; Wu, Y. Topology optimization of steady and unsteady incompressible Navier-Stokes flows driven by body forces. Struct. Multidiscip. Optim. 2013, 47, 555–570. [Google Scholar] [CrossRef]
  19. Deng, Y.; Liu, Z.; Zhang, P.; Liu, Y.; Wu, Y. Topology optimization of unsteady incompressible Navier-Stokes flows. J. Comput. Phys. 2011, 230, 6688–6708. [Google Scholar] [CrossRef]
  20. Gersborg-Hansen, A.; Sigmund, O.; Haber, R.B. Topology optimization of channel flow problems. Struct. Multidiscip. Optim. 2005, 30, 181–192. [Google Scholar] [CrossRef]
  21. Olesen, L.H.; Okkels, F.; Bruus, H. A high-level programming language implantation of topology optimization applied to steady-state Navier-Stokes flow. Int. J. Numer. Methods Eng. 2006, 65, 975–1001. [Google Scholar] [CrossRef]
  22. Deng, Y.; Liu, Z.; Zhang, P.; Liu, Y.; Gao, Q.; Wu, Y. A flexible layout design method for passive micromixers. Biomed. Microdevices 2012, 14, 929–945. [Google Scholar] [CrossRef] [PubMed]
  23. Zhou, T.; Wang, H.; Shi, L.; Liu, Z.; Joo, S. An enhanced eletroosmotic micromixer with an efficient asymmetric lateral structure. Micromachines 2016, 7, 218. [Google Scholar] [CrossRef]
  24. Ji, Y.; Deng, Y.; Liu, Z.; Zhou, T.; Wu, Y.; Qian, S. Optimal control-based inverse determination of electrode distribution for electroosmotic micromixer. Micromachines 2017, 8, 247. [Google Scholar] [CrossRef]
  25. Andreasen, C.; Gersborg, A.; Sigmund, O. Topology optimization of microfluidic mixers. Int. J. Numer. Meth. Fluids 2009, 61, 498–513. [Google Scholar] [CrossRef]
  26. Singh, M.K.; Kang, T.G.; Meijer, H.E.H.; Anderson, P.D. The mapping method as a toolbox to analyze, design, and optimize micromixers. Microfluid. Nanofluid 2008, 5, 313–325. [Google Scholar] [CrossRef]
  27. Singh, M.K.; Galaktionov, O.S.; Meijer, H.E.H.; Anderson, P.D. A simplified approach to compute distribution matrices for the mapping method. Comput. Chem. Eng. 2009, 33, 1354–1362. [Google Scholar] [CrossRef]
  28. Welander, P. Studies on the general development of motion in a two-dimensional ideal fluid. Tellus 1955, 7, 141–156. [Google Scholar] [CrossRef]
  29. Danckwerts, P.V. The definition and measurement of some characteristics of mixtures. Appl. Sci. Res. 1952, 3, 279–296. [Google Scholar] [CrossRef]
  30. Kruijt, P.G.M.; Galaktionov, O.S.; Anderson, P.D.; Peters, G.W.M.; Meijer, H.E.H. Analyzing mixing in periodic flows by distribution matrices: Mapping method. AIChE J. 2001, 47, 1005–1015. [Google Scholar] [CrossRef]
  31. Kruijt, P.G.M.; Galaktionov, O.S.; Peters, G.W.M.; Meijer, H.E.H. The mapping method for mixing optimizatiom. Part II Transport in a Corotating Twin Screw Extruder. Int. Polym. Proc. 2001, 16, 161–171. [Google Scholar] [CrossRef]
  32. Du, Y.; Zhang, Z.; Yim, C.H.; Lin, M.; Cao, X. A simplified design of the staggered herringbone micromixer. Biomicrofluidics 2010, 4, 024105. [Google Scholar] [CrossRef] [PubMed]
  33. Williams, K.J.; Longmuir, K.J.; Yager, P. A practical guide to the staggered herringbone mixer. Lab Chip 2008, 8, 1121–1129. [Google Scholar] [CrossRef] [PubMed]
Figure 1. The calculation of the mapping coefficient Φij in the mapping matrix is the ratio of the number of markers received by the recipient cell Ωi at X = X0 + ΔX to the initial number of markers in Ωj at X = X0 (in this example Φij = 3 /16) [26,27].
Figure 1. The calculation of the mapping coefficient Φij in the mapping matrix is the ratio of the number of markers received by the recipient cell Ωi at X = X0 + ΔX to the initial number of markers in Ωj at X = X0 (in this example Φij = 3 /16) [26,27].
Micromachines 09 00137 g001
Figure 2. Schematic representation of the grooves in a staggered herringbone mixer (SHM). The mapping matrix Φ1 covers a single groove applying a fully developed velocity field.
Figure 2. Schematic representation of the grooves in a staggered herringbone mixer (SHM). The mapping matrix Φ1 covers a single groove applying a fully developed velocity field.
Micromachines 09 00137 g002
Figure 3. The evolution of coarse grain concentration distribution Ci in a SHM with one groove type: (a) The result of mapping method which the tracing cross section is covered with 50 × 50 cells, and each cell is filled with 100 uniformly distributed markers; (b) The result of mapping method which the tracing cross section is covered with 100 × 100 cells, and each cell is filled with 100 uniformly distributed markers; (c) The results in the study of Singh et al. which the tracing cross section is covered with 200 × 200 cells, and each cell is filled with 256 uniformly distributed markers [26]; (d) The coarse grain concentration distribution at outlet of SHM with ten grooves.
Figure 3. The evolution of coarse grain concentration distribution Ci in a SHM with one groove type: (a) The result of mapping method which the tracing cross section is covered with 50 × 50 cells, and each cell is filled with 100 uniformly distributed markers; (b) The result of mapping method which the tracing cross section is covered with 100 × 100 cells, and each cell is filled with 100 uniformly distributed markers; (c) The results in the study of Singh et al. which the tracing cross section is covered with 200 × 200 cells, and each cell is filled with 256 uniformly distributed markers [26]; (d) The coarse grain concentration distribution at outlet of SHM with ten grooves.
Micromachines 09 00137 g003
Figure 4. (a) Design domain ΩD of the micromixers at the bottom layer of the straight channel. ΩC is the channels connected to the inlet, outlet and design domain ΩD. The length of design domain ΩD is n times of L; (b) the distribution of coarse grain concentration on the outlet.
Figure 4. (a) Design domain ΩD of the micromixers at the bottom layer of the straight channel. ΩC is the channels connected to the inlet, outlet and design domain ΩD. The length of design domain ΩD is n times of L; (b) the distribution of coarse grain concentration on the outlet.
Micromachines 09 00137 g004
Figure 5. (a) Layout of the micromixer for the design domain as shown in Figure 4, where n = 2, Re = 2.5, and projected velocity vector distribution in the cross sections (S1, S2, S3, S4, S5); the mixing measurement J C 0 ( C ) = 0.9141; (b) the distribution of coarse grain concentration on the outlet.
Figure 5. (a) Layout of the micromixer for the design domain as shown in Figure 4, where n = 2, Re = 2.5, and projected velocity vector distribution in the cross sections (S1, S2, S3, S4, S5); the mixing measurement J C 0 ( C ) = 0.9141; (b) the distribution of coarse grain concentration on the outlet.
Micromachines 09 00137 g005
Figure 6. (a) Layout of the micromixer for the design domain as shown in Figure 4, where n = 2, Re = 5, and projected velocity vector distribution in the cross sections (S1, S2, S3, S4, S5); the mixing measurement J C 0 ( C ) = 0.9147; (b) the distribution of coarse grain concentration on the outlet.
Figure 6. (a) Layout of the micromixer for the design domain as shown in Figure 4, where n = 2, Re = 5, and projected velocity vector distribution in the cross sections (S1, S2, S3, S4, S5); the mixing measurement J C 0 ( C ) = 0.9147; (b) the distribution of coarse grain concentration on the outlet.
Micromachines 09 00137 g006
Figure 7. (a) Layout of the micromixer for the design domain as shown in Figure 4, where n = 4, Re = 5, and projected velocity vector distribution in the cross sections (S1, S2, S3, S4, S5, S6, S7, S8, S9); the mixing measurement J C 0 ( C ) = 0.8382; (b) the distribution of coarse grain concentration on the outlet.
Figure 7. (a) Layout of the micromixer for the design domain as shown in Figure 4, where n = 4, Re = 5, and projected velocity vector distribution in the cross sections (S1, S2, S3, S4, S5, S6, S7, S8, S9); the mixing measurement J C 0 ( C ) = 0.8382; (b) the distribution of coarse grain concentration on the outlet.
Micromachines 09 00137 g007aMicromachines 09 00137 g007b
Figure 8. The evolution of mixing performance of micromixers with different cycles which the single cycle is the layout in Figure 5: (ad) are the distribution of coarse grain concentration on the outlet with 1, 2, 5, and 10 cycles, and J C ¯ ( C ) are 0.9882, 0.9779, 0.9404 and 0.8413, respectively.
Figure 8. The evolution of mixing performance of micromixers with different cycles which the single cycle is the layout in Figure 5: (ad) are the distribution of coarse grain concentration on the outlet with 1, 2, 5, and 10 cycles, and J C ¯ ( C ) are 0.9882, 0.9779, 0.9404 and 0.8413, respectively.
Micromachines 09 00137 g008
Figure 9. (a,c) are the distribution of coarse grain concentration on the outlet with 1 and 10 cycles shown in Figure 8a and d, and J C ¯ ( C ) are 0.9882 and 0.8413; (b,d) are the distribution of coarse grain concentration on the outlet of SHM (Figure 10) with 1 and 10 cycles, and J C ¯ ( C ) are 0.9883 and 0.8446.
Figure 9. (a,c) are the distribution of coarse grain concentration on the outlet with 1 and 10 cycles shown in Figure 8a and d, and J C ¯ ( C ) are 0.9882 and 0.8413; (b,d) are the distribution of coarse grain concentration on the outlet of SHM (Figure 10) with 1 and 10 cycles, and J C ¯ ( C ) are 0.9883 and 0.8446.
Micromachines 09 00137 g009
Figure 10. (a) Layout of the SHM in the design domain as shown in Figure 4, where n = 2, Re = 2.5, and the mixing measurement J C ¯ ( C ) = 0.9883; (b) the distribution of coarse grain concentration on the outlet. (the SHM has only one groove and is same volume as the micromixer shown in Figure 5).
Figure 10. (a) Layout of the SHM in the design domain as shown in Figure 4, where n = 2, Re = 2.5, and the mixing measurement J C ¯ ( C ) = 0.9883; (b) the distribution of coarse grain concentration on the outlet. (the SHM has only one groove and is same volume as the micromixer shown in Figure 5).
Micromachines 09 00137 g010
Table 1. Procedure of the iterative approach for solving the topology optimization problem.
Table 1. Procedure of the iterative approach for solving the topology optimization problem.
  • Give the initial value of the design variable γ;
  • Solve the Navier-Stokes equations and backward particle tracing equation by the finite element method;
  • Solve the adjoint equation;
  • Compute the adjoint derivative and the corresponding objective and constraint values;
  • Update the design variable by method of moving asymptotes (MMA);
  • Check for convergence; if the stopping conditions are not satisfied, go to 2; and
  • Post-processing

Share and Cite

MDPI and ACS Style

Guo, Y.; Xu, Y.; Deng, Y.; Liu, Z. Topology Optimization of Passive Micromixers Based on Lagrangian Mapping Method. Micromachines 2018, 9, 137. https://doi.org/10.3390/mi9030137

AMA Style

Guo Y, Xu Y, Deng Y, Liu Z. Topology Optimization of Passive Micromixers Based on Lagrangian Mapping Method. Micromachines. 2018; 9(3):137. https://doi.org/10.3390/mi9030137

Chicago/Turabian Style

Guo, Yuchen, Yifan Xu, Yongbo Deng, and Zhenyu Liu. 2018. "Topology Optimization of Passive Micromixers Based on Lagrangian Mapping Method" Micromachines 9, no. 3: 137. https://doi.org/10.3390/mi9030137

APA Style

Guo, Y., Xu, Y., Deng, Y., & Liu, Z. (2018). Topology Optimization of Passive Micromixers Based on Lagrangian Mapping Method. Micromachines, 9(3), 137. https://doi.org/10.3390/mi9030137

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop