Volume 19 Number 21 7 November 2019 Pages 3565–3736

Lab on aChip
Devices and applications at the micro- and nanoscale
rsc.li/loc

ISSN 1473-0197

PAPER
J. Wang, N. Zhang et al.
Finding the optimal design of a passive microfluidic mixer
 Lab on a Chip
View Article Online
PAPER View Journal | View Issue
Published on 17 September 2019. Downloaded by Indian Institute of Technology Roorkee on 7/4/2023 1:19:47 PM.

Finding the optimal design of a passive

Cite this: Lab Chip, 2019, 19, 3618
microfluidic mixer†
Junchao Wang, ‡ab Naiyin Zhang, ‡c Jin Chen,a Victor G. J. Rodgers,b
Philip Briskd and William H. Grover *b

The ability to thoroughly mix two fluids is a fundamental need in microfluidics. While a variety of different
microfluidic mixers have been designed by researchers, it remains unknown which (if any) of these mixers
are optimal (that is, which designs provide the most thorough mixing with the smallest possible fluidic
resistance across the mixer). In this work, we automatically designed and rationally optimized a microfluidic
mixer. We accomplished this by first generating a library of thousands of different randomly designed
mixers, then using the non-dominated sorting genetic algorithm II (NSGA-II) to optimize the random chips
in order to achieve Pareto efficiency. Pareto efficiency is a state of allocation of resources (e.g. driving
force) from which it is impossible to reallocate so as to make any one individual criterion better off (e.g.
pressure drop) without making at least one individual criterion (e.g. mixing performance) worse off. After
200 generations of evolution, Pareto efficiency was achieved and the Pareto-optimal front was found. We
examined designs at the Pareto-optimal front and found several design criteria that enhance the mixing
performance of a mixer while minimizing its fluidic resistance; these observations provide new criteria on
how to design optimal microfluidic mixers. Additionally, we compared the designs from NSGA-II with some
Received 9th June 2019, popular microfluidic mixer designs from the literature and found that designs from NSGA-II have lower
Accepted 16th September 2019
fluidic resistance with similar mixing performance. As a proof of concept, we fabricated three mixer designs
from 200 generations of evolution and one conventional popular mixer design and tested the performance
DOI: 10.1039/c9lc00546c
of these four mixers. Using this approach, an optimal design of a passive microfluidic mixer is found and
rsc.li/loc the criteria of designing a passive microfluidic mixer are established.

1 Introduction in contact). Active mixers generally outperform passive

mixers, but integrating an external force or field in the chip
Mixing is one of the fundamental functions in microfluidic adds unwanted complexity and cost. Passive mixers are
chips. For the past decade, a wide variety of different simpler and more economical, but increasing the area and
microfluidic mixers have been designed.1 Microfluidic mixers time of contact between the two fluids has undesirable
are usually categorized as either “active” (an external energy consequences: increasing contact area by lengthening the
force or an external physical field is present to accelerate channel containing the two fluids adds unwanted additional
mixing phenomenon) or “passive” (mixing is accomplished fluidic resistance to the channel, and increasing contact time
only by diffusion and is dependent only on the area of contact by slowing the flow rate decreases the overall throughput of
between the two fluids and the amount of time the fluids are the microfluidic chip.2 Thus, there is an unmet need for
mixer designs that combine high mixing performance with
low fluidic resistance and high flow rates.
a
Key Laboratory of RF Circuits and Systems, Ministry of Education, and Zhejiang Several studies have been conducted on the optimization
Provincial Laboratory of Integrated Circuit Design, Hangzhou Dianzi University,
of standard microfluidic mixer designs. Li et al. optimized a
China. E-mail: junchao@hdu.edu.cn; Tel: +86 571 8691 9078
b
Department of Bioengineering, University of California Riverside, Riverside, CA,
chaotic microfluidic mixer using lattice Boltzmann method.3
USA. E-mail: wgrover@engr.ucr.edu; Tel: +1 951 827 4311 Hertzog et al. used an optimized microfluidic mixer to study
c
College of Life Information Science and Instrument Engineering, Hangzhou the protein folding kinetics.4 Two continuous studies from
Dianzi University, China Wang et al. focused on the optimization of the layout of
d
Department of Computer Science and Engineering, University of California
obstacles for enhanced mixing in microchannels for different
Riverside, Riverside, CA, USA
† Electronic supplementary information (ESI) available. See DOI: 10.1039/
applications using a fluid dynamics software.5,6 Hossain
c9lc00546c et al. conducted research of optimizing a modified Tesla
‡ These authors contributed equally to this work. structure based on topology optimization.7 Finally, Cortes-

3618 | Lab Chip, 2019, 19, 3618–3627 This journal is © The Royal Society of Chemistry 2019
 View Article Online

Lab on a Chip Paper

Quiroz et al. optimized a grooved microfluidic mixer using a give us new useful design criteria to use when manually
multi-objective optimization approach.8 In these optimization designing mixers?
processes, the design criteria of the post-optimized mixer In this work, we set out to answer the question, is it
designs remained unchanged compared to their original possible to find the most optimized mixer within certain
designs, which results in a limited improvement of the conditions? Specifically, are we able to explore the
mixing performance. For instance, Hossain et al. optimized performance boundary of how good a microfluidic mixer can
Published on 17 September 2019. Downloaded by Indian Institute of Technology Roorkee on 7/4/2023 1:19:47 PM.

the mixing performance of the Tesla structure mixer, but the possibly be within a certain limit on fluidic resistance?
basic design of the mixer remained unchanged.7 Here, we developed an approach to automatically design
Occasionally, researchers develop new microfluidic mixer and optimize passive microfluidic mixers for specific
designs that have advantages over existing designs. For conditions. We accomplished this in two steps. First, we
example, Fu et al. designed a rapid vortex microfluidic mixer generated a library of more than six thousand different
utilizing double-heart chambers that achieved a 92% mixing random mixer designs and simulated the performance of
ratio at Reynolds numbers as low as Re = 1 (ref. 9) and Wang each of them. We have previously used this technique to
et al. used triangular posts in a conventional Y-shaped generate designs of functional microfluidic chips that can
microfluidic mixer to increase the mixing performance.10 But deliver solutes of any desired concentrations.11 Second, we
are these mixer designs really optimal, or are there better used the non-dominated sorting genetic algorithm II (NSGA-
designs waiting to be discovered? With an infinite variety of II)12 to optimize multiple design parameters of our
possible designs, and only a relatively small number of microfluidic mixer at the same time. NSGA-II is one of the
researchers exploring this design space, progress toward multi-objective evolutionary algorithms (MOEAs), which
better mixers is frustratingly slow. helped our mixer designs to increase their mixing
This situation inspired us to ask, is it possible to design a performance, achieve Pareto efficiency, and find the Pareto-
microfluidic mixer from scratch by computer algorithm optimal front. Pareto efficiency is a state of allocation of
without needing microfluidics expertise at design phase? If resources (e.g. fluid driving force) from which it is
so, is it possible that these automatically-designed mixers will impossible to reallocate so as to make any one individual

Fig. 1 (A) Schematic of a simulated microfluidic mixer unit. A simulated unit has two inlets and two outlets. Between inlets and outlets is a 500
μm × 500 μm design domain. In the design domain, each mixer has ten cylindrical posts with random sizes and locations. Different cylinder posts
were allowed to overlap to create additional structures like walls. (B) The predicted fluid velocity field of a typical mixer unit. This velocity field is
used for simulating the solute concentration profile in the mixer. (C) The predicted pressure profile of the mixer unit. This pressure profile is used
to characterize the fluidic resistance of the mixer. (D) The predicted solute concentration profile of the mixer unit. This concentration profile is
used to determine the mixing performance of this mixer unit.

This journal is © The Royal Society of Chemistry 2019 Lab Chip, 2019, 19, 3618–3627 | 3619
 View Article Online

Paper Lab on a Chip

criterion better off (e.g. pressure drop) without making one generated designs as well as NSGA-II designs. The specific
or more individual criterion (e.g. mixing performance) worse code we write to generate 6069 mixer designs is available in
off. The Pareto-optimal front is the set of all Pareto efficient ESI.†
allocations, which is conventionally visualized as a boundary
in a graph of performance. To investigate the Pareto 2.2 Simulating mixer performance
efficiency of our system, random mixer designs and NSGA-II
Published on 17 September 2019. Downloaded by Indian Institute of Technology Roorkee on 7/4/2023 1:19:47 PM.

All simulations were performed using the finite element

mixer designs were visualized in the same graph (Fig. 4A) of analysis software COMSOL Multiphysics (COMSOL Inc.,
mixing performance vs. fluidic resistance. Burlington, MA). We used the software's MATLAB API to
Our results showed that after 200 generations of evolution, automate the simulation process. The laminar flow physics
the mixer designs converged near the true Pareto-optimal module and transport of dilute species physics module as
front; this allowed us to explore the fundamental performance well as two stationary solvers were used in COMSOL
limits of a microfluidic mixer. A user can select a design from Multiphysics. In the laminar flow physics module in
these optimized mixers and be confident that the design is COMSOL Multiphysics, each inlet was assigned an inlet
optimal for a given fluidic resistance. Additionally, we boundary condition of 1 mm s−1 normal inflow velocity, and
identified certain design trends in the optimized mixers, each outlet was assigned an outlet boundary condition of 0
manually designed several mixers that incorporate these Pa pressure. The remaining boundaries were walls (no-slip
trends, and compared the performance of our manually- boundary condition), and the material filling the channels
designed mixers to that of our automatically-designed optimal was water under incompressible flow. In the transport of
mixers. In each case, our automatically-designed and dilute species physics module, inlet 1 is assigned an inflow
optimized mixers equaled or exceeded the mixing concentration of 1 mmol L−1 and inlet 2 is assigned an inflow
performance of conventional designed mixers. Finally, to concentration of 0 mmol L−1. The two outlets were assigned
confirm that the mixers designed by our algorithm function as as outflows. The solute diffusion coefficient of fluorescein
predicted, we chose three optimum mixer designs at (4.25 × 10−10 m2 s−1) was used in simulation in order to
corresponding minimized resistance conditions from the represent the mixing behavior of small molecules.13
Pareto-optimal front and one chip from conventional designed Fig. 1B and C show the calculated velocity field and pressure
mixers and fabricated corresponding polydimethylsiloxane field of one mixer unit design, and Fig. 1D shows the
(PDMS) microfluidic chips for experimental verification. concentration mixing field of the same design. The
corresponding script for simulating the performance of the
2 Materials and methods mixer designs is available in ESI.†
2.1 Generating initial random mixer designs
We created our first generation of passive microfluidic mixers 2.3 Evolving mixer designs with NSGA-II
by generating mixer designs at random.11 Of course, there is The genetic algorithm NSGA-II12 was used to evolve
an essentially limitless variety of possible mixer designs, so optimized versions of our passive random mixers. A flow
we applied certain constraints to our random designs. Fig. 1A chart representation of our custom NSGA-II implementation
shows the basic design template of our random mixers. Each is shown in Fig. 2. The fitness function for fluidic resistance
mixer has two inlets, two outlets, and a 500 μm × 500 μm (SP) is
design domain where the random mixing structures are
located. In the design domain are ten cylindrical posts with SP = P2 − P1 (1)
random sizes and locations. Ten cylindrical structures were
chosen as a balance between computational resources and where P2 is the pressure at the outlets and P1 is the pressure
achieving as many different mixing features as possible at the inlets. This means that the smaller the pressure drop
within a 500 μm × 500 μm design domain. Each cylindrical across the mixer, the better the performance of the mixture.
structure acted as a building unit of the mixing features in The fitness function for mixing performance is defined as the
each design, and ten cylindrical structures were good enough mixing score (SC),
to represent the diversity of mixing features. For example,
two or more cylindrical posts can overlap, which enables the SC = (1 − C1) + (C2 − 0) (2)
mixer designs to also include non-circular features (like
walls, inverted-L or S-shape). In addition, cylindrical where C1 is the average concentration of outlet 1; 1 − C1
structures naturally create smoother fluid streamlines than calculates the average concentration difference between inlet
square or triangular structures, which is helpful for reducing 1 and outlet 1; C2 is the average concentration of outlet 2; C2
the overall fluidic resistance of the mixer. In total, 6096 − 0 calculates the average concentration difference between
different mixer designs were generated and stored in a inlet 2 and outlet 2. This indicates that the closer the
database. Finally, in addition to randomly-generated designs, concentrations of the fluids in outlet 1 and outlet 2, the
we also manually designed five mixer units based on our better performance of the mixer. Eqn (2) is not our only
experience so as to compare them with the randomly choice for quantifying the mixing performance. For example,

3620 | Lab Chip, 2019, 19, 3618–3627 This journal is © The Royal Society of Chemistry 2019
 View Article Online

Lab on a Chip Paper

Out of all 6069 random mixer designs, one design was

chosen as generation 0 (G0, the parent of first generation in
our evolutionary algorithm). G0 was chosen for two reasons.
First, G0 is one of the top performers, which is located at the
edge of all random designs in the pressure drop vs. mixing
score map (Fig. 4A), a location that might already be close to
Published on 17 September 2019. Downloaded by Indian Institute of Technology Roorkee on 7/4/2023 1:19:47 PM.

the Pareto-optimal front. Second, after investigating the post

layout of G0, we found all ten cylindrical posts were located
around the center region of the design domain, which might
let our algorithm to have a higher probability to explore as
many post layouts as possible in a limited number of evolution
runs. The non-dominant sorting operation, selection, crossover,
and mutation operators were then conducted so as to generate
the new population of designs. After that, numerical
simulations were performed using the same simulation
parameters as the randomly generated designs above. In total,
200 generations were calculated to find the Pareto-optimal
front. The corresponding guidelines for implementing NSGA-II
for specific applications are available in ESI.†

2.4 The robustness of NSGA-II

NSGA-II is one of the most popular multi-objective
evolutionary algorithms (MOEAs), and the algorithm has
been applied in many different fields.14–17 To verify the
performance and robustness of MOEAs, computer scientists
have developed guidelines and carefully selected a number of
test problems.18,19 The objective functions of these test
problems are complicated mathematical functions, whose
graphs could be convex, nonconvex, disconnected or even
Fig. 2 A flow chart depicting our custom NSGA-II process for nonuniformly spaced. Computer scientists then used MOEAs
optimizing mixers and finding the Pareto-optimal front. The overall
to find the minimum solutions which were satisfying all the
goal was to minimize the pressure drop (fluidic resistance) of mixer
designs while increasing the mixing performance. Numerical simulation objective functions at the same time. After testing with nine
was conducted by COMSOL Multiphysics and MATLAB. Typical genetic test problem sets, NSGA-II demonstrated strong robustness
operators (selection, crossover, mutation) were conducted after a among different MOEAs.12
non-dominant sorting operator. After that, the population of next
generation mixer designs were generated and repeated in the loop
until the optimization criterion was satisfied. 2.5 Functional chip design, fabrication and experiments
Design generation 0 (G0), generation 60 (G60), generation 120
(G120) and conventional design C (shown in Fig. 4) were
((1 − C1) + (C2 − 0))2 or C1/C2 are acceptable fitness functions chosen to be fabricated by conventional soft-lithography.20
as well. However, eqn (2) has several advantages in this study. Since each mixer unit (Fig. 1) predicted by our method has
For example, eqn (2) results in a range of scores that is limited mixing performance in the 500 μm × 500 μm design
normalized from 0 to 1. For instance, a perfectly mixed domain, we put 11 identical mixing units in a chain to
solution would have C1 = C2 = 0.5 mmol L−1, SC would equal amplify the mixing performance. As shown in Fig. 3A and C,
1; a perfectly un-mixed solution would have C1 = 1 mmol L−1, each fabricated microfluidic mixer has 11 mixing units and
C2 = 0 mmol L−1, SC would equal 0. Additionally, since eqn each mixing unit duplicates the structure of the design
(1) had been defined as a first-order equation, it was natural domain from G0, G60, G120, and conventional design C.
to define the fitness function for mixing performance as a Treating all 11 copies of each mixer identically is not ideal.
first-order equation as well because both optimization criteria However, our decision to do so represents a trade-off between
are equally important in our optimization process. Using two practical issues. On one hand, since the amount of
fitness function like ((1 − C1) + (C2 − 0))2 for mixing mixing provided by a single mixer is rather small and would
performance would probably yield similar results. However, be difficult to accurately quantify experimentally, we needed a
plots of the Pareto-optimal front would be distorted by the way to amplify the net mixing we observed in our
non-uniform spacing of mixing scores. For these reasons, experiments. On the other hand, we recognize that the second
eqn (2) was used to quantify the mixing performance in this mixer in a series receives a slightly different concentration
work. profile than the first mixer, the third mixer receives a different

This journal is © The Royal Society of Chemistry 2019 Lab Chip, 2019, 19, 3618–3627 | 3621
 View Article Online

Paper Lab on a Chip

Published on 17 September 2019. Downloaded by Indian Institute of Technology Roorkee on 7/4/2023 1:19:47 PM.

Fig. 3 (A) The SU-8 mold of mixer designs G0, G60, G120, and conventional design C. (B) In each case inlet 1 was injected with FD&C Blue #1 and
inlet 2 was injected with FD&C Red #3. Silastic laboratory tubings were used to connect the syringe pump and sample reservoirs with the
microfluidic mixers. Post-mixed fluids were collected in these reservoirs for quantification of the mixing performance. (C) Photographs of the G60
microfluidic mixer. As with all designs, this G60 microfluidic mixer has two inlets and two outlets. The main channel consists of 11 identical mixing
units, each of which corresponds to the structure of pre-simulated G60 design. (D) Photographs of the first mixing unit and the last mixing unit of
G60 microfluidic mixer in operation. Two dyes (red color from the top inlet, and orange color from the bottom inlet) flow into the chip. The
different layouts created by the 10 posts in the design domain affected the mixing performance and the resistance of each mixer design, causing
the two fluids to be significantly mixed after passing through 11 mixing units and exiting through the two outlets.

profile than the second, and so on. Accurately predicting the pressures at each inlet will give similar results, since the two
unique behavior of the N unit mixer would require knowledge sets of inlet and outlet channels in our mixer chips have the
of the unique behavior of all N-1 unit mixers upstream. This same length and resistance to flow). The volumetric flow rate
interdependence would add enormous computational for each syringe was set to 3 μL min−1. The inlets of chips
complexity to our task: finding an optimal series of different were connected to the syringes by Silastic laboratory tubing
mixer designs would require simulating and optimizing (Dow Corning, Michigan, US) and the outlets were connected
thousands of versions of all 11 designs connected together, a to the sample reservoirs by tubing as well. Fig. 3D shows the
computational task that is far outside of our capabilities. first mixing unit and the last mixing unit of G60 microfluidic
As a trade-off between experimental verification and mixer in operation. For visualization, inlet 1 contains a dye
computational feasibility, we experimentally tested different that appears orange on our imaging system (FD&C Red #3)
mixer designs using 11 identical copies of each mixer. and inlet 2 contains a dye that appears red (FD&C Blue #1).
The main features (cylindrical posts) of the mixer were The chips were imaged using an optical microscope
designed by our algorithms and exported into DXF files.21 (Olympus BX51, Tokyo, Japan). For quantification of the
Based on these DXF files, additional features including inlets mixing performance, inlet 1 was injected with water and inlet
and outlets were designed manually in AutoCAD (Autodesk, 2 was injected with FD&C Red #3. To make sure that the
San Rafael, CA) and then written to a transparent mask. system reached steady state and all the air bubbles went
Negative photoresist (SU-8 25, Microchem, MA) was spin- away, samples were collected from the chip after 5 minutes
coated on a 4 inch polished silicon wafer to fabricate the SU- of flow and visual confirmation for the nonexistence of
8 mold as shown in Fig. 3A. The channel width was 200 μm bubbles in the channels was conducted with the microscope.
and the channels depth was 50 μm, which was consistent After that, samples from both outlets were collected into
with our simulation models. After that, a volumetric ratio of tubes. A standard curve was plotted based on standard dye
10 : 1 mixture of PDMS (Sylgard 184, Dow Corning, MI) and concentration. Both the standard curve and samples were
curing agent were poured onto the SU-8 mold. After analyzed using a UV-VIS-NIR spectrophotometer at 530 nm
degassing and curing, the PDMS replica was peeled off from (Shimadzu UV3600, Kyoto, Japan). Finally, the concentration
the master and punched on top for inlet and outlet. Finally, a of each sample was calculated using the standard curve.
plasma cleaner was used to change the surface properties of
the PDMS replica and the glass slides in order to create a 3 Results and discussion
PDMS–glass bond. 3.1 Finding the Pareto-optimal front
As shown in Fig. 3B, a dual channel syringe pump was
used to provide the driving force for the microfluidic mixers Fig. 4A plots pressure drop versus mixing score for each of
(although we expect that pressure-driven flow using identical the randomly-generated designs (small blue circles), NSGA-

3622 | Lab Chip, 2019, 19, 3618–3627 This journal is © The Royal Society of Chemistry 2019
 View Article Online

Lab on a Chip Paper

Published on 17 September 2019. Downloaded by Indian Institute of Technology Roorkee on 7/4/2023 1:19:47 PM.

Fig. 4 (A) Pressure drop (SP) vs. mixing score (SC) for random microfluidic mixer designs (blue dots) and NSGA-II-evolved mixer designs (red stars).
The random designs are distributed in the bottom-left corner while the NSGA-II designs are at the boundary of all random designs. By connecting
all the NSGA-II designs, we can draw a Pareto-optimal front (dashed line). Conventional designs (B–F) (yellow stars) lie above the Pareto-optimal
front, which means that their mixing performance is not as good as NSGA-II designs with a certain pressure drop. (G0–G200) The concentration
profiles of NSGA-II designs of generation 0, 25, 50, 75, 100, 125, 150, 175 and 200. Additional concentration profiles, pressure profiles, and velocity
fields of 0–200 generations are available in ESI.†

II-evolved optimal designs (red stars), and conventional Conventional design C has a pressure drop of 5.36 Pa and a
designs (yellow stars). The NSGA-II designs distribute at the mixing score of 0.59 mmol L−1. Conventional design D has a
boundary of the randomly generated designs. This means pressure drop of 1.29 Pa and a mixing score of 0.38 mmol
that NSGA-II successfully found the Pareto-optimal front. L−1. Adding their mixing performance to the plot in
To achieve a similar mixing score, NSGA-II designs always Fig. 4A (yellow stars) shows that designs B, C, and D all lie
need less pressure drop or generate less resistance in a above the evolved designs (red stars). This tells us that the
mixer unit. In other words, within a certain pressure drop common mixer designs for microfluidics still have potential
condition, the NSGA-II designs will always have better to be optimized.
mixing performance than the random-design mixers. Since Fig. 4 (G0–G200) are the concentration profiles of NSGA-II
G0 was randomly designed, its pressure drop still had designs in generations 0, 25, 50, 75, 100, 125, 150, 175 and
space to be minimized. That is why we observed a small 200. As the generation number increases, the mixing
decrease in pressure drop during the initial 20 generations. performance improves and the mixer geometry converges
After that, as the mixing score increased, the pressure drop into an S-shaped line of cylinders. The S-shape suggests that
increased as well. NSGA-II selects S-shaped designs as elite designs and retains
Fig. 4B–D are three common microfluidic mixer designs the S-shaped feature into the next generations. The S-shape
being constrained to our design domain using cylinder posts could increase the mixing contact area as well as minimizing
to map the geometry. Conventional design B has a pressure the fluidic resistance. The small gaps between each post also
drop of 0.98 Pa and a mixing score of 0.35 mmol L−1. appear to be crucial to the performance of the mixer. From

This journal is © The Royal Society of Chemistry 2019 Lab Chip, 2019, 19, 3618–3627 | 3623
 View Article Online

Paper Lab on a Chip

the concentration and pressure profiles of each generation

(see ESI†), we know that each small gap allows fluid with no
chance to mix (solute concentration around 0 mmol L−1) to
go through the S-shape and reduce the overall fluidic
resistance of the mixer. We are unaware of any similar
designs that have been created by conventional manual
Published on 17 September 2019. Downloaded by Indian Institute of Technology Roorkee on 7/4/2023 1:19:47 PM.

design methods. Finally, Fig. 4E and F are manually-designed

mixers that are inspired by NSGA-II designs. Conventional
design E has a pressure drop of 2.87 Pa and a mixing score is
0.59 mmol L−1. Conventional design F has a pressure drop of
0.81 Pa and a score of 0.40 mmol L−1. Their performance
(gold stars marked E and F on Fig. 4A) is close to the Pareto-
optimal front but they do not have small gaps in dark blue
area (solute concentration around 0 mmol L−1) to reduce the
fluidic resistance.

3.2 Experimental verification

To demonstrate the functionalities of the automatically-
designed microfluidic mixers, three evolved designs (G0, G60,
and G120 in Fig. 4) and one of the conventional designs (C in
Fig. 4) were chosen to be fabricated and tested. The
concentration profiles of each last (11th) mixing unit of G0,
G60, G120, and conventional design C in operation are shown
in Fig. 5A. As generation number increased, the red fluid
occupied more region in the design domain, invading the area
of the orange fluid, which enhanced the mixing phenomenon.
The small gaps created by posts in the red fluid region
contributed to minimizing the fluidic resistance. The mixing
in conventional design C simply relied on diffusion between
two fluid without any structure to minimize the fluidic
resistance. As shown in Fig. 5B, the mixing scores of G0, G60,
G120 and conventional design C were 0.459 mmol L−1, 0.498
mmol L−1, 0.613 mmol L−1 and 0.762 mmol L−1, respectively;
and the predicted pressure drop of the 11 mixing units by
Fig. 5 (A) The concentration profiles of each last (11th) mixing unit of
COMSOL Multiphysics were 5.5 Pa, 5.39 Pa, 8.36 Pa and 58.96 G0, G60, G120, and conventional design C in operation. (B) Mixing
Pa, respectively. After evolving for 120 generations, the mixing scores and predicted pressure drop of four fabricated mixers. The left
score improved by 33.6% while the cost (the pressure drop y-axis indicates the mixing score and the right y-axis indicates the
generated by 11 mixing units) increased by 52%. In contrast, predicted pressure drop of total 11 mixing units of each mixer. Three
measurements were conducted for each point of the mixing scores;
though conventional design C had a better mixing score than
error bars indicate ±1 standard deviation.
G120 (0.762 mmol L−1 vs. 0.613 mmol L−1), the pressure drop
cost was tremendous (58.96 Pa vs. 8.36 Pa).
In order to mathematically quantify the performance of
the mixers using the mixing score and the pressure drop, eqn is clear that although conventional design C has a higher
(3) is defined as follows, mixing score than G0, G60 or G120, its Mcost is more than 5
times higher than the Mcost of G0, G60 and G120 in average.
SP The Mcost first dropped 10% from G0 to G60 and increased
M cost ¼ (3)
SC 25% from G60 to G120. Since G0 was randomly designed and
only located close to the Pareto-optimal front, the 10% drop
where SP is the pressure drop that is defined in eqn (1) and indicated that G0 to G60 was approaching the ideal Pareto-
calculated by COMSOL Multiphysics; SC is the mixing score optimal front, in which the increase of SC was along with the
defined in eqn (2) and measured by experiments; Mcost decrease of SP in a certain range. The 25% increase indicated
indicates the mixing cost of the mixers by calculating the
fraction of the pressure drop and the mixing score. The Table 1 Mixing cost of four tested mixers
physical meaning of Mcost is how much pressure (the driving
force) we need in order to achieve 1.0 mmol L−1 mixing score. G0 G60 G120 Conventional design C
−1
The Mcost of four tested mixers is summarized in Table 1. It Mcost (Pa mmol L) 12.0 10.9 13.7 76.9

3624 | Lab Chip, 2019, 19, 3618–3627 This journal is © The Royal Society of Chemistry 2019
 View Article Online

Lab on a Chip Paper

that the mixer designs of G60 to G120 had lain on the Pareto- (around 0.68 mmol L−1), the design from the first run has a
optimal front, in which the increase of SC had to be along lower pressure drop. This indicates that the Pareto-optimal
with the increase of SP as well and the increment of SP is front found from the first run is closer to the ideal Pareto-
larger than the increment of SC. Overall, the Mcost gives us a optimal front. Although the geometries resulting from the two
quantitative way to calculate how much a specific mixer evolution runs are different, they do share two important
design can be optimized. similarities. First, they both created a narrow gap near the left
Published on 17 September 2019. Downloaded by Indian Institute of Technology Roorkee on 7/4/2023 1:19:47 PM.

edge with a large cylindrical post. Second, they used the rest
of the cylinder posts to generate a wall containing small gaps
3.3 Rational design inspired by two separate evolution runs in the dark blue area (around 0 mmol L−1) so as to minimize
Since we only have 200 populations in each generation while the pressure drop. So, why did the design from the first
the size and position of cylinder posts in the design domain evolution run have a lower pressure drop? From the
are infinite, the Pareto-optimal front we found is close to the concentration profiles, we can see that in the first run design,
ideal Pareto-optimal front theoretically. It is interesting to fluid had a longer contact time and contact area before
investigate how separate evolution runs will affect the design entering the critical gap (generated by the largest post).
of mixers. Fig. 6 shows the comparison of NGSA-II designs Additionally, it seems that the second-evolution designs only
between two separate evolution runs. Fig. 6A is the first run used eight posts to create a wall instead of ten. Two upper
(the results are the same as those in Fig. 4). Fig. 6B is the cylinder posts (pointed by gray arrows) seem to have no
second run, and in this run we found that the geometry function but increase the fluidic resistance of this design.
converged into a Y-shape instead of an S-shape. While the The main reason two separate runs falls into two different
mixing scores of these two separate evolution runs are similar local minimums is due to the limited populations in each
generation. In order to get highly identical results between
separated runs, we could include more populations in each
generation. However, predicting the velocity fields and
concentration profiles of more populations would be
computationally expensive. It took us several hours to
simulate 200 different velocity fields and concentration
profiles for only one generation even the simulations were
processed in a workstation with a Intel 10-core Xeon Silver
CPU and 64 GB RAM. Fortunately, the comparison between
two separated runs gave us a perfect example – if we want to
reduce the overall pressure drop, we could design some gaps
in the region where fluid get no chance to mix without
hurting the overall mixing performance.

4 Conclusions
We demonstrated how to optimize a functional microfluidic
mixer for two parameters, pressure drop and mixing score,
using NSGA-II. We accomplished this by using MATLAB and
COMSOL Multiphysics as our simulation platform and
implementing NSGA-II in MATLAB. We found the pressure
drop versus mixing score Pareto-optimal front. After that, we
compared the designs at the Pareto-optimal front with
conventional designs and random designs. Our simulations
indicate that designs from NSGA-II have lower pressure drops
than designs by conventional methods or random designs
while achieving a similar mixing performance. Based on the
NSGA-II designs, we have a better understanding about how to
design a microfluidic mixer rationally: a mixer should have a
constriction to increase contact area and contact time between
Fig. 6 (A) The NSGA-II design selected at the end of the first run of
evolution. (B) The NSGA-II design selected at the end of a second run
the fluids, as well as some features that are not for mixing but
of evolution. In the second run, the geometry converged into a rather for reducing the overall resistance of the mixer.
Y-shape. To achieve a similar mixing score as the S-shaped design
from run 1, the Y-shaped design from run 2 will have a higher fluidic 4.1 Limitations
resistance. Two gray arrows indicate the inefficient use of two posts by
the NSGA-II algorithm, which seemed to only increase the resistance The optimum mixer designs generated by our algorithm have
instead of improving the mixing performance. certain constraints. Since the boundary conditions of our

This journal is © The Royal Society of Chemistry 2019 Lab Chip, 2019, 19, 3618–3627 | 3625
 View Article Online

Paper Lab on a Chip

optimizing system were set to a fixed value, the optimum the target cells.31 In this case, coupled with our previous
mixer design is only investigated and experimentally verified work (MOPSA, microfluidics-optimized particle simulation
when Reynold number is around 0.4. Our current work algorithm32), NSGA-II could be used to optimize an inertial
utilizes a constant volumetric flow rate to drive fluid flow. microfluidic chip so as to increase the separation
Since the two sets of inlet and outlet channels in our mixer performance while minimizing the damage to cells from
chips have the same length and resistance to flow, we expect shear stress.
Published on 17 September 2019. Downloaded by Indian Institute of Technology Roorkee on 7/4/2023 1:19:47 PM.

that our results would also hold true for pressure-driven flow
using identical pressures at each inlet. It is possible that the Conflicts of interest
optimal mixing features (the layout of posts) will change
according to different boundary conditions. For instance, a There are no conflicts of interest to declare.
certain optimum mixer design may still have potential to be
optimized if we apply it on inertial microfluidics, centrifugal Acknowledgements
microfluidics or capillary microfluidics. As for other methods This work was supported by National Natural Science
for driving fluid flow (like inertial, centrifugal, or capillary Foundation of China (No. 61827806 and 61871161).
forces), as long as we can accurately model the physical
phenomena involved in those methods, our micromixer References
generating algorithm can optimize those microfluidics
designs based on the predicted results. Accurately modeling 1 C.-Y. Lee, C.-L. Chang, Y.-N. Wang and L.-M. Fu, Int. J. Mol.
inertial, centrifugal, and capillary forces in microfluidics is Sci., 2011, 12, 3263–3287.
itself a challenging and active research topic, one that is 2 C.-Y. Lee, W.-T. Wang, C.-C. Liu and L.-M. Fu, Chem. Eng. J.,
beyond the scope of this work. But as improved modeling 2016, 288, 146–160.
techniques are developed in the future, these techniques can 3 C. Li and T. Chen, Sens. Actuators, B, 2005, 106,
be integrated with our micromixer generating algorithm to 871–877.
develop optimal mixer designs for systems that utilize those 4 D. E. Hertzog, B. Ivorra, B. Mohammadi, O. Bakajin and
forces. J. G. Santiago, Anal. Chem., 2006, 78, 4299–4306.
In addition, our simulations were constrained in a 2-D 5 H. Wang, P. Iovenitti, E. Harvey and S. Masood,
plane, and the height of the channel was not defined. In J. Micromech. Microeng., 2003, 13, 801.
theory, a 3-D simulation would be a more accurate way to 6 H. Wang, P. Iovenitti, E. Harvey and S. Masood, Smart
capture the mixing and diffusion phenomenon. For instance, Mater. Struct., 2002, 11, 662.
Stroock et al. designed a chaotic micromixer with a 7 S. Hossain, M. A. Ansari, A. Husain and K.-Y. Kim, Chem.
herringbone design in the bottom of the channels that could Eng. J., 2010, 158, 305–314.
introduce mixing of different lamina in 3-D.22 Cha et al. 8 C. A. Cortes-Quiroz, A. Azarbadegan, M. Zangeneh and A.
presented a PDMS-based micromixer with 3-D structures for Goto, Chem. Eng. J., 2010, 160, 852–864.
rapid mixing performance.23 However, the computational 9 L.-M. Fu, W.-C. Fang, H.-H. Hou, Y.-N. Wang and T.-F.
expense of implementing such a large number of simulations Hong, Chem. Eng. J., 2014, 249, 246–251.
in 3-D was prohibitive for us. Therefore, the performance of 10 L. Wang, S. Ma, X. Wang, H. Bi and X. Han, Asia-Pac. J.
3-dimensional micromixers could not be compared with our Chem. Eng., 2014, 9, 877–885.
optimum mixer designs. 11 J. Wang, P. Brisk and W. H. Grover, Lab Chip, 2016, 16,
4212–4219.
12 K. Deb, A. Pratap, S. Agarwal and T. Meyarivan, IEEE Trans.
4.2 Future directions Evol. Comput., 2002, 6, 182–197.
Our approach is not limited to optimize only mixer 13 Z. Petrášek and P. Schwille, Biophys. J., 2008, 94, 1437–1448.
performance—it should be able to optimize additional 14 B. Ghasemishabankareh, M. Ozlen and X. Li, International
parameters as well. For instance, the overall chip size is also Conference on Evolutionary Multi-Criterion Optimization,
a key aspect of a microfluidic chip. Instead of constraining a 2019, pp. 541–552.
mixer into a fixed design domain, we could try to minimize 15 M. Hesami, R. Naderi, M. Tohidfar and M. Yoosefzadeh-
the size of the design domain as well. Additionally, Najafabadi, Front. Plant Sci., 2019, 10, 869.
microfluidic mixers are just one of many components in 16 M. H. Shojaeefard, S. E. Hosseini and J. Zare, Struct.
microfluidic chips.24,25 We are confident that our approach Multidiscip. Optim., 2019, 1–17.
can be applied to other applications of microfluidics. For 17 S. Wang, D. Zhao, J. Yuan, H. Li and Y. Gao, Electr. Power
example, cell sorting is a major application in Syst. Res., 2019, 175, 105893.
microfluidics.26 In various sorting technologies, inertial 18 K. Deb, Evol. Comput., 1999, 7, 205–230.
microfluidics shows the potential to efficiently separate 19 E. Zitzler, K. Deb and L. Thiele, Evol. Comput., 2000, 8,
different cells based on the sizes of cells.27–30 However, 173–195.
inertial microfluidic devices are usually operated at a high 20 Y. Xia and G. M. Whitesides, Angew. Chem., Int. Ed.,
Reynolds number, at which shear stress could be harmful to 1998, 37, 550–575.

3626 | Lab Chip, 2019, 19, 3618–3627 This journal is © The Royal Society of Chemistry 2019
 View Article Online

Lab on a Chip Paper

21 Autodesk, DXF Reference, 2011, http://images.autodesk.com/ 27 D. Di Carlo, Lab Chip, 2009, 9, 3038–3046.
adsk/files/autocad_2012_pdf_dxf-reference_enu.pdf. 28 S. S. Kuntaegowdanahalli, A. A. S. Bhagat, G. Kumar and I.
22 A. D. Stroock, S. K. Dertinger, A. Ajdari, I. Mezić, H. A. Stone Papautsky, Lab Chip, 2009, 9, 2973–2980.
and G. M. Whitesides, Science, 2002, 295, 647–651. 29 S. C. Hur, N. K. Henderson-MacLennan, E. R. McCabe and
23 J. Cha, J. Kim, S.-K. Ryu, J. Park, Y. Jeong, S. Park, S. Park, D. Di Carlo, Lab Chip, 2011, 11, 912–920.
H. C. Kim and K. Chun, J. Micromech. Microeng., 2006, 16, 30 A. A. S. Bhagat, H. W. Hou, L. D. Li, C. T. Lim and J. Han,
Published on 17 September 2019. Downloaded by Indian Institute of Technology Roorkee on 7/4/2023 1:19:47 PM.

1778. Lab Chip, 2011, 11, 1870–1878.

24 G. M. Whitesides, Nature, 2006, 442, 368. 31 M. G. Lee, J. H. Shin, C. Y. Bae, S. Choi and J.-K. Park, Anal.
25 T. M. Squires and S. R. Quake, Rev. Mod. Phys., 2005, 77, 977. Chem., 2013, 85, 6213–6218.
26 C. W. Shields IV, C. D. Reyes and G. P. López, Lab Chip, 32 J. Wang, V. G. Rodgers, P. Brisk and W. H. Grover,
2015, 15, 1230–1249. Biomicrofluidics, 2017, 11, 034121.

This journal is © The Royal Society of Chemistry 2019 Lab Chip, 2019, 19, 3618–3627 | 3627