1 Introduction

Microfluidic devices are used in a wide range of applications in biological (Hansen and Quake 2003; Beebe et al. 2002) and chemical (Ehrfeld et al. 2000; Reyes et al. 2002) analyzes. Mixing is of general importance in most microfluidic applications, but is often difficult to achieve since, in microfluidic devices, flows are generally laminar, and turbulence, as commonly used in macro-mixing devices, is absent. Without other means to enhance mixing, mixing by pure diffusion requires long times or long flow lengths. Laminar flow can, however, also produce complicated trajectories of fluid particles, resulting in efficient mixing via chaotic advection (Aref 1984; Ottino 1989). The principles of chaotic advection, which mimic the baker’s transformation: a continuous repetition of stretching and folding, are exploited in different static industrial mixing devices. In microfluidic devices, it is generally difficult to incorporate complex geometries as used in static macro mixers. Consequently, different strategies are applied: (1) to induce transverse flow in pressure driven flows by adopting a three-dimensional static structure in a channel geometry in so-called passive micromixers, and (2) to induce transverse flow by applying external sources for example pressure drop, temperature gradient, acoustic pressure, and magnetic fields in active micromixers.

Bertsch et al. (2001) proposed a miniaturized static micromixer based on a large-scale industrial static mixer geometry. Other researchers proposed “split-recombine” micro mixers (Schönfeld et al. 2004; Park et al. 2004; Jen et al. 2003), which mimic the baker’s transformation and split, stretch, and combine a flow to achieve efficient mixing and to produce uniform striations and lamella. Stroock et al. (2002) experimentally investigated mixing in a three-dimensional staggered herringbone micromixer (henceforth, SHM) and in a three-dimensional slanted groove micromixer (henceforth, SGM). The SHM and the SGM have patterned grooves on the bottom of a rectangular channel to produce helical fluid trajectories inside the channels. Experiments show that the SHM works well in the range of 0 < Re < 100 and that it produces chaotic advection. Kim et al. (2004) proposed an other modification of the SGM using simple grooves combined with barriers on top of the channel known as the barrier-embedded micromixer (BEM). The SHM and the BEM make use of a specific three-dimensional structure to induce a lateral motion of fluids in a periodic manner, giving rise to chaotic mixing. Liu et al. (2000) use a three-dimensional serpentine channel (henceforth, 3D-SC) to achieve chaotic mixing. In this device, mixing depends on inertia that causes the secondary flow. In the Stokes flow regime (Re ≪ 1), where inertia is absent, the device becomes inefficient. Nguyen and Wu (2005) give a comprehensive description of all these chaotic micromixers.

Understanding mixing in these micromixers was also improved by numerical analyzes based on particle tracking (Kang and Kwon 2004; Aubin et al. 2003, 2005; Stroock and McGraw 2004). In spite of the accuracy and the superiority compared to numerical schemes based on the solution of the mass transport equation (which suffer from numerical diffusion if mesh resolution is not fine enough to capture complicated deformations of fluid domains), mixing analysis via particle tracking has several inherent drawbacks. First, it requires the tracking of a huge number of particles to generate high-resolution images at far downstream locations. Second, the time-consuming particle tracking procedure must be repeated at any change in the sequence of repeating units of the mixer. Finally, there is no guarantee that all the space of interest at the desired location will be completely occupied by particles, due to the fact that any ordered array of particles at the inlet become disordered at downstream positions. This leads to the loss of accuracy in quantification of mixing based on particle distributions. In this study, therefore, we apply the mapping method to overcome these inherent disadvantages. Note that the mapping method is only used to analyze distributive mixing, i.e., flows with high Péclet number (Pe≫1) where diffusive mixing contributions are negligible. To extend potential applications of the mapping method, we introduce a new numerical formulation of the method and apply it to analyze and optimize micromixers. In fact, the new approach is far more straight forward as compared to the original mapping approach (Anderson and Meijer 2000; Galaktionov et al. 1997, 2002, 2003; Kruijt et al. 2001a, b), and moreover, it is easily implemented.

We focus on three distinct passive micromixers, the SHM, the BEM, and the 3D-SC, respectively. The three chosen micromixers are specifically considered to show the various aspects of the method as an analysis, design, and optimization tool. For the SHM, we perform in-depth analysis of mixing and show use of the mapping method as an optimization tool. For the BEM, we employ the method to investigate whether different ordering of mixing protocols (functional modules) achieve better mixing; for example, aperiodic sequences can generate better mixing in some cases (Kang et al. 2007). The 3D-SC differs from the two other examples mentioned in the sense that inertia at higher Reynolds numbers induces a back-flow (negative axial velocity) in the channel, while in the Stokes regime in some parts of the channel the axial velocity approaches zero. This specifically requires time-tracking approach to track particles in the flow field. From this analysis we show that the method can be applied as an analysis tool to evaluate mixing in many diverse types of micromixers. We start with describing the basics of the method, and address the detailed procedures of implementation for the three micromixers. For each micromixer, we first introduce the geometrical parameters associated with the mixer, then compute the mapping matrices for each of the characteristic modules of the mixer, and then finally analyze mixing in qualitative, and quantitative way using the method.

2 Mapping method

2.1 Basics of mapping method and original mapping approach

Chaotic mixing of viscous liquids in laminar flows is usually based on the situation where the baker’s transformation is applied a number of times on a specified volume of material. Spencer and Wiley (1951) suggested that the distribution of material in such flows can be handled quite well by the use of matrix methods. The mapping method describes the transport of a conservative quantity without taking into account diffusion from one state to another by means of a mapping matrix, describing the transport of fluid from an initial cross section to a final one (for spatially periodic flows) or from an initial time to a final time (for time-periodic flows). A micromixing device composed of repetitive sequences (repeating units), the advantage of using mapping method becomes apparent, since it requires just a one-time computation of the deformation induced by the flow during a fixed flow in the time Δt (for time-periodic flows) or fixed flow in the length Δl (for spatially-periodic flows). The effect of the flow for any number or combination of cycles (t = n Δt or l = n Δl) can then be evaluated by a repeated multiplication (n − times) of the distribution matrix with a prescribed concentration distribution (vector) at the inlet. Since these multiplications take only few CPU seconds, this brings a huge benefit over conventional particle tracking techniques in which the tracking is cumbersomely repeated from the first to the last period to analyze mixing. In addition, in the mapping method, a change in the sequence of repetitive units in a mixer requires just a different order of distribution matrices to be multiplied, whereas the conventional approaches require the full re-computation of the particle tracking. The method also allows to compute quantitative measures of mixing, such as the volume- or area-averaged, and flux-weighted discrete intensity of segregation, and the scale of segregation (Danckwerts 1952; Galaktionov et al. 2003).

Numerically, the original mapping method (Anderson and Meijer 2000;Galaktionov et al. 1997, 2002, 2003; Kruijt et al. 2001a, b) exploits the above idea as follows: a distribution matrix \({\varvec{\varphi}}\) is formed to store information about the distribution of fluid from one cross-section to the next due to a specified flow. To obtain the coefficients of the distribution matrix, the initial cross section of the flow domain is subdivided into a large number of discrete cells (N) of identical size. During flow, the material from a donor cell is transferred to different recipient cells. The fraction of material that is transferred from the donor cell to a recipient cell gives the distribution coefficient of the donor cell with respect to the recipient cell. Thus, in total N cells form a distribution matrix of the order N × N. The discrete coefficient φ ij equals the fraction of deformed sub-domain Ω j at z = z 0 + Δz that is found in the original sub-domain Ω i at z = z 0:

$$ \varphi_{ij} = \frac{\int_{\Upomega_j|_{z=z_0+\Updelta z} \bigcap \Upomega_i|_{z=z_0}} \hbox{d}A}{\int_{\Upomega_j|_{z=z_0}} \hbox{d}A}. $$
(1)

Tracking all interfaces of all N cells during a flow over a distance Δz can be done as we have demonstrated for different flows (Anderson and Meijer 2000; Galaktionov et al. 1997, 2002, 2003; Kruijt et al. 2001a, b), but it is cumbersome to track interfaces experiencing complicated deformation patterns. Therefore, here we present an alternative approach that is much simpler to implement.

2.2 A new formulation of the mapping method

A schematic representation of how the mapping coefficients are calculated in new formulation of the method is shown in Fig. 1. To approximate the coefficients of the mapping matrix (or distribution matrix), K markers inside all cells are tracked. The markers are uniformly distributed in the cells. Then, to determine the final distribution of markers, they are advected during the flow from z = z 0 to z = z 0 + Δz. If the number of markers in the donor cell Ω j is M j at z = z 0 and the number of markers found after tracking in the recipient cell Ω i is M ij at z = z 0 + Δz, then the mapping coefficient Φ ij is calculated as:

$$ \Upphi_{ij} = \frac{M_{ij}}{M_j}. $$
(2)

In other words, the coefficient Φ ij is the measure of the fraction of total flux of the cell Ω j donated to the cell Ω i . If the number of markers tracked is large enough then Φ ij approaches φ ij .

Fig. 1
figure 1

Illustration of the computation the mapping coefficient Φ ij in the mapping matrix \({\varvec{\Upphi}}.\) The cell Ω j at z = z 0 is covered with a number of markers that are tracked during flow in Δz (to arrive at the final cross section z = z 0 + Δz). The ratio of the number of markers received by the recipient cell Ω i to the initial number of markers in Ω j is determined (in this example Φ ij is 6/25)

Full size image

The elegance of this mapping method is that if one wants to analyze mixing-related scalar quantities, like the concentration vector \({\user2{C}} \in {{\mathbb{R}}}^{N\times 1}\) (N is the number of cells) defined on initial cells, then the concentration evolution C 1 after the deformation can be obtained by simply multiplying the mapping matrix \({\varvec{\Upphi}}\) with the initial concentration vector C 0:

$$ {\user2{C}}^{1} = {\varvec{\Upphi}}{\user2{C}}^0. $$
(3)

Note that C represents coarse-grained description of volume fraction (dimensionless concentration) of a marker fluid in a mixture of two marker fluids with identical material properties, and its component C i describes the concentration (volume fraction) locally averaged in the cell Ω i . For repetitive mixing, the same operation is repeated multiple times on the same mass and, hence, the concentration evolution after n steps is given by \({\user2{C}}^n ={\varvec{\Upphi}}^n {\user2{C}}^0.\) For sufficiently large n, the matrix \({\varvec{\Upphi}}^n\) will not be sparse and it becomes so large that it can even not be stored anymore. This is due to the fact that after performing the operation n times, material from one cell is advected to a large part of the whole cross-section, especially in the case of chaotic advection. Instead of studying \({\varvec{\Upphi}}^n,\) the evolution of the concentration after n steps C n is computed in sequence as follows:

$$ {\user2{C}}^{i+1} = {\varvec{\Upphi}}{\user2{C}}^i,\quad \hbox{hence} \quad {\user2{C}}^n= (\underbrace{{\varvec{\Upphi}}({\varvec{\Upphi}}(\ldots ({\varvec{\Upphi}}}_{n\;{\rm times}}{\user2{C}}^0)))). $$
(4)

Thus, the mapping matrix \({\varvec{\Upphi}}\) is determined only once and is utilized a number of times to study the evolution of concentration in the flow field. The computation of mapping matrices is expensive, and may take several CPU hours, but, once calculated, the necessary matrix–vector multiplication only takes a few CPU seconds to process the results. The mapping matrix calculations are easily parallelized (Galaktionov et al. 1997).

2.3 Measure of mixing

To quantify mixing and to compare the performance of mixers, we employ the intensity of segregation as a measure of mixing defined as the second-moment variance of concentration distribution (Danckwerts 1952):

$$ I = \frac{\sigma_c^2}{{\bar{c}}(1-{\bar{c}})}, $$
(5)

where σ 2 c is the variance in the concentration over entire domain Ω defined as:

$$ \sigma_c^2 = {\left\langle \left(c(x)-{\bar{c}}\right)^2 \right\rangle}_\Upomega, $$
(6)

where c(x) denotes the volume fraction (dimensionless concentration) of a fluid in a mixture of two fluids at a point x, and \({\bar{c}}\) the average volume fraction of the fluid in the whole domain. When no diffusion is present, c(x) will either be 1 or 0. Therefore, I will always be equal to 1, independent of the distribution. To avoid this situation, the coarse grain concentration C i (Welander 1955) on a finite cell Ω i is defined:

$$ {C}_i = {\left\langle c(x) \right\rangle}_{\Upomega_i}. $$
(7)

In the coarse grain description of concentration, C i can take values between and including, 0 and 1. From Eq. (5), we define the flux-weighted discrete intensity of segregation (I d ) for distributive mixing without diffusion. The number of cells (sub-domains) used to compute the discrete intensity of segregation is chosen equal to the number of cells used to compute the mapping matrices. If the cells are uniform in size, and thus the flux-weighted discrete intensity of segregation can be simplified as follows:

$$ I_d = \frac{1}{{\bar{C}}(1-{\bar{C}})}\frac{1} {F}\sum^N_{i=1}(C_i-{\bar{C}})^2 f_i, $$
(8)

where the average concentration \({\bar{C}}\) is

$$ {\bar{C}} = \frac{1}{F}\sum\limits_{i=1}^NC_i f_i,\quad F=\sum\limits_{i=1}^Nf_i. $$
(9)

The term f i is the volumetric flux through cell number i and F is the total flux through the mixer. The intensity of segregation I d is a measure of the deviation of the local concentration from the ideal situation (perfectly mixed case), which represents a homogeneous state of the mixture. In a perfectly mixed system, I d = 0, while in a completely segregated system, I d = 1. As found by Galaktionov et al. (2002, 2003) the flux-weighted definition (see Eq. (8)) of the intensity of segregation is much better suited for analyzing continuous mixers than area- or volume-averaged definitions of the intensity of segregation. This is due the fact that the real influence of an unmixed island on the value of I d is proportional to the flux, carrying this island. Note that the number of cells covering a cross-sectional plane in the mapping method must be fixed for comparing various mixers and numerous layouts of a mixer. This is due to the fact that the coarse graining mixing measure I d is dependent on the cell size of the mapping. The cell size tells about the minimum striation thickness between two mixing fluids which can be resolved, and below this size the fluids are assumed to be completely mixed. In our study we have used 200 × 200 grid to cover a cross section of interest to analyze various layouts.

3 Application to the staggered herringbone micromixer

First, we apply the mapping method to analyze and optimize the design of SHMs. The SHM geometry is subdivided into four functional mixing modules, and for each module a mapping matrix is computed. Combining these four matrices in different ways enables us to investigate various designs of the SHM.

3.1 Mixer geometry

Figure 2 shows, schematically, a two-dimensional view of the grooves on the bottom of the rectangular channel in a SHM. Their presence induces transverse flows inside the channel. Every groove is composed of two parts, a long arm and a short arm (length ratio of arms considered here is 2:1), both at 45° to the axial direction; two series of grooves comprise a single cycle. We apply the geometry of the SHM used in the study of Kang and Kwon (2004) with a channel height h = 77 μm and a channel width w = 200 μm, the depth of the grooves g d = 17.7 μm, their width g w = 70.7 μm and the offset distance between two consecutive grooves in a half cycle also equals 70.7 μm. The original SHM consists of six grooves per half cycle. Five different SHM geometries are investigated with groove depths g d of 10, 20, 30, 40, and 50 μm. Fluent V6.1.22 is used to obtain the periodic velocity field. Structured hexahedral meshes are constructed by using Gambit V2.1.2. The material properties of a mixture of glycerol (80%) and water (20%) with a density ρ = 1,200 kg/m3 and viscosity μ = 0.067 kg m/s, are prescribed for all the simulations. The fluid is assumed to be Newtonian and incompressible. To obtain a fully developed velocity field for one cycle, periodic boundary conditions are prescribed with the mass flow rate corresponding to an average inlet velocity u of 0.2 cm/s. The Reynolds number Re (=ρu D h /μ, where D h is hydraulic diameter) is found to be below 10−2 for all simulations.

Fig. 2
figure 2

Schematic representation of the grooves in the bottom of a rectangular channel of a SHM. a Definition of a cycle. The mapping matrices \({\varvec{\Upphi}}_{\rm L}\;\hbox{and}\;{\varvec{\Upphi}}_{\rm R}\) cover a single groove applying a fully developed velocity field. b The mapping matrices \({\varvec{\Upphi}}_{t{\rm LR}}\;\hbox{and}\;{\varvec{\Upphi}}_{t{\rm RL}}\) cover the two transition regions, where flow is not developed due to entrance and exit effects caused by changing symmetry of two types of grooves (the geometrical features used are: h = 77 μm, w = 200 μm and g d  = 17.7 μm, g w = 70.7 μm, and length of one periodic unit (cycle) = 1,992 μm)

Full size image

3.2 Defining the mapping matrices of the SHM

To implement the mapping approach for optimization, the SHM is subdivided into four independent functional mixing modules for which the mapping matrices are denoted by \({\varvec{\Upphi}}_{\rm L}, {\varvec{\Upphi}}_{t{\rm LR}}, {\varvec{\Upphi}}_{\rm R}\;\hbox{and}\; {\varvec{\Upphi}}_{t{\rm RL}}\) (see Fig. 2). The mapping matrix for a single groove in the first half cycle is denoted by \({\varvec{\Upphi}}_{\rm L}\) (index L for left). Similarly, mapping matrix in the second half cycle is denoted by \({\varvec{\Upphi}}_{\rm R}\) (index R for right). Both are computed based on a fully developed velocity field. The grid used to obtain the periodic velocity field for one cycle consists of 567,840 hexahedral cells and 609,070 nodes. When the grooves of the SHM switch symmetry, a transition region occurs: first, halfway a cycle, and second, between cycles (see Fig. 2b). In the transition regions, the velocity field is not developed and comprises all exit and entrance effects when two types of grooves in the SHM switch symmetry. Transition regions can play an important role in mixing (Galaktionov et al. 2002, 2003), and hence, they must be included in the analyzes. The mapping matrix for the first transition is denoted by \({\varvec{\Upphi}}_{t{\rm LR}}\) (index small t for transition), and the one for the second transition is denoted by \({\varvec{\Upphi}}_{t{\rm RL}}.\) Any SHM can be mapped by a combination of these four matrices \({\varvec{\Upphi}}_{\rm L}, {\varvec{\Upphi}}_{t{\rm LR}}, {\varvec{\Upphi}}_{\rm R}\; \hbox{and}\;{\varvec{\Upphi}}_{t{\rm RL}}.\) Since the two sets of grooves in two half cycles and the two transitions are mirror images, there is no need to calculate \({\varvec{\Upphi}}_{\rm R}\;\hbox{and}\;{\varvec{\Upphi}}_{t{\rm RL}}\) via mapping, since they can be computed via mirroring \({\varvec{\Upphi}}_{\rm L}\;\hbox{and}\; {\varvec{\Upphi}}_{t{\rm LR}},\) respectively.

In the mapping computations, the cross section of interest is covered with a 200 × 200 grid, and each cell is filled with 256 uniformly distributed (in a 16 × 16 pattern) passive markers (compare with Fig. 1). For a fixed number of cells covering a cross-sectional plane in the mapping method, the mixing index I d should converge with the increase in the number of particles per cell (NPPC). In case of the SHM, we show the dependence of flux-weighted intensity of segregation on the NPPC (see Fig. 3). From this plot, it is clear that the NPPC above a critical value (in this case 64) provides converged mixing measure. In our computations, we used 256 number of particles per cell. The critical value of the NPPC may be effected by nature of flow, whether it is chaotic or regular, and how a mixer is subdivided into various mixing modules to compute separate individual matrices representative of the mixing modules.

Fig. 3
figure 3

The dependence of mixing index I d on the NPPC (up to 10th mixing cycle) for the SHM (the geometrical features used are: h = 77 μm, w = 200 μm and g d  = 17.7 μm, g w  = 70.7  μm, and length of one periodic unit (cycle) = 1,992 μm)

Full size image

The trajectory of markers is tracked by using the axial co-ordinate, rather than time. This can be realized by dividing the transversal velocity components u x and u y by the axial velocity component u z :

$$ \frac{\hbox{d}x}{\hbox{d}z} = \frac{u_{x}}{u_{z}},\quad \frac{\hbox{d}y} {\hbox{d}z} = \frac{u_{y}}{u_{z}}. $$
(10)

This axial integration approach is useful because integration is done with respect to the spatial increment along the axial direction rather than time, eliminating the effects of different residence time distributions. Note that this approach is only valid for the systems where back-flows are not present. It turned out to be advantageous to use the backward (reverse) particle tracking (BPT) to track the tracers to obtain the mapping coefficients. In other words, tracers originally filling the recipient cell are tracked backward against the flow direction. Equation (10) is integrated by the fourth order Runge–Kutta Bulrish Store scheme with the adaptive step size selection of Press et al. (1992). To find the velocity at any arbitrary point, interpolation using the basis function is applied (Galaktionov et al. 1997, 2002, 2003).

3.3 Combining mapping matrices to achieve various designs

The computed mapping matrices \({\varvec{\Upphi}}_{\rm L}, {\varvec{\Upphi}}_{t{\rm LR}}, {\varvec{\Upphi}}_{\rm R}\;\hbox{and}\;{\varvec{\Upphi}}_{t{\rm RL}}\) are combined to obtain the concentration distribution for a SHM with a desired number of grooves per half cycle for a number of cycles. This provides a simple and a computationally inexpensive way to evaluate different designs. To analyze the concentration evolution in a SHM with 10 grooves per half cycle, the first cycle, C 1, is obtained with:

$$ {\user2{C}}^1 = \underbrace{({\varvec{\Upphi}}_{\rm L}({\varvec{\Upphi}}_{\rm L}(\ldots}_ {8\; {\rm times}}({\varvec{\Upphi}}_{t{\rm LR}}\underbrace{({\varvec{\Upphi}}_{\rm R} ({\varvec{\Upphi}}_{\rm R}(\ldots({\varvec{\Upphi}}_{\rm R}}_{8\;{\rm times}}{\user2{C}}^0)))))))), $$
(11)

where C 0 is the prescribed initial concentration distribution. We neglect the inlet and outlet effects of the total system. This simplification is obviously less severe for longer SHM geometries (e.g., 20 cycles). Note that the transition regions \({\varvec{\Upphi}}_{t{\rm LR}}\;\hbox{and}\;{\varvec{\Upphi}}_{t{\rm RL}}\) contain two \({\varvec{\Upphi}}_{\rm L}\) and two \({\varvec{\Upphi}}_{\rm R}\) grooves to compensate for local entrance and exit effects at these transitions. These contributions must be taken into account and hence the equation for C 1 contains only 8 (rather than 10) grooves per half cycle. To increase readability we introduce a simpler notation:

$$ {\varvec{\Upphi}}_{n{\rm L}}= (\underbrace{{\varvec{\Upphi}}_{\rm L}({\varvec{\Upphi}}_{\rm L}(\ldots}_{n\;{\rm times}}, $$
(12)

where \({\varvec{\Upphi}}_{n{\rm L}}\) represents the matrix-vector multiplication of matrix \({\varvec{\Upphi}}_{\rm L}\) in sequence for n-times on a given initial concentration distribution. Note that here notation \({\varvec{\Upphi}}_{n{\rm L}}\) is used instead of \({\varvec{\Upphi}}^{n{\rm L}}\) to describe the matrix-vector multiplication in a sequence (recall Eq. (4)). Similarly, \({\varvec{\Upphi}}_{n{\rm R}}\) is defined. In calculating the concentration evolution for the second cycle C 2, three contributions from transition regions are met: two \({\varvec{\Upphi}}_{t{\rm LR}}\) and one \({\varvec{\Upphi}}_{t{\rm RL}}.\) Hence, the number of grooves per half cycle in intermediate parts (6) is calculated by taking into account the number of grooves (2) that are part of these (2) transition regions. Consequently, the concentration evolution for cycle 2, C 2, and that of cycle 3, C 3, can be obtained as follows:

$$ {\user2{C}}^2 = ({\varvec{\Upphi}}_{\rm 8L}\underbrace{({\varvec{\Upphi}}_{t{\rm LR}}({\varvec{\Upphi}} _{\rm 6R}({\varvec{\Upphi}}_{t{\rm RL}}({\varvec{\Upphi}}_{\rm 6L}}_{\rm repeating\; unit}({\varvec{\Upphi}}_{t{\rm LR}}({\varvec{\Upphi}}_{\rm 8R}{\user2{C}}^0))))))), $$
(13)

and

$$ {{\user2{C}}^3 = ({\varvec{\Upphi}}_{\rm 8L}\underbrace{({\varvec{\Upphi}}_{t{\rm LR}}({\varvec{\Upphi}}_ {\rm 6R}({\varvec{\Upphi}}_{t{\rm RL}}({\varvec{\Upphi}}_{\rm 6L}}_{\rm repeating\; unit}} \underbrace{({\varvec{\Upphi}}_{t{\rm LR}}({\varvec{\Upphi}}_{\rm 6R}({\varvec{\Upphi}}_ {t{\rm RL}}({\varvec{\Upphi}}_{\rm 6L}}_{\rm repeating\; unit}({\varvec{\Upphi}}_{t{\rm LR}}({\varvec{\Upphi}}_{\rm 8R}{\user2{C}}^0)))))))))). $$
(14)

One easily recognizes the repeating unit, which we denote by \({\varvec{\Upphi}}_{\rm RU}\) (index RU for repeating unit). To calculate concentration evolutions of other cycles, this notation can be used. For example, C 3 can be re-written, and C 4 can be expressed using the repeating unit as follows:

$$ {\user2{C}}^3 = ({\varvec{\Upphi}}_{\rm 8L}(2\times{\varvec{\Upphi}}_{\rm RU}({\varvec{\Upphi}}_{t{\rm LR}} ({\varvec{\Upphi}}_{\rm 8R}{\user2{C}}^0)))), $$
(15)
$$ {\user2{C}}^4 = ({\varvec{\Upphi}}_{\rm 8L}(3\times{\varvec{\Upphi}}_{\rm RU}({\varvec{\Upphi}}_{t{\rm LR}}( {\varvec{\Upphi}}_{\rm 8R}{\user2{C}}^0)))), $$
(16)

where, \(2\times{\varvec{\Upphi}}_{\rm RU}\;\hbox{and}\; 3\times{\varvec{\Upphi}}_{\rm RU}\) represents two- and three-times repetition of the repeat unit \({\varvec{\Upphi}}_{\rm RU}.\) Hence, from the above, it is clear that one can calculate concentration evolutions in the SHM with 10 numbers of grooves per half cycle for any number of cycles. The concentration evolution for other designs of the micromixer with a different number of grooves per half cycle can be obtained in a similar way (basically only changing the numbers 8 and 6 in above equations) and mixing quality can be characterized quantitatively using the intensity of segregation. This is our basis of optimization the number of grooves later. From the above it is clear that the minimum number of grooves per half cycle that can be analyzed equals 4, since that equals the number of grooves involved in the entrance and exit regions of the transition regions. To clarify the above steps, here we show mixing evolutions in a simple SHM consisting of only one type of groove whose mapping matrix is represented by \({\varvec{\Upphi}}_{\rm L}\) (as shown in Fig. 2a). The concentration distribution C 1 after 1 groove is obtained by multiplying the matrix \({\varvec{\Upphi}}_{\rm L}\) with the concentration vector C 0, next C 2 is found multiplying \({\varvec{\Upphi}}_{\rm L}\) with C 1, etc. In this way, the mixing evolutions obtained are shown in Fig. 4, illustrating the basic cross-sectional deformation induced by the presence of the asymmetric grooves.

Fig. 4
figure 4

Evolution of the concentration distribution C i in a SHM with only one groove type, shown after the ith groove and computed by repeated multiplication of matrix \({\varvec{\Upphi}}_{\rm L}\) with evolved concentration vectors, \({\user2{C}}^{i+1}= {\varvec{\Upphi}}_{\rm L}{\user2{C}}^i\) (the geometrical features used are: h = 77 μm, w = 200 μm and g d  = 17.7 μm, g w  = 70.7 μm, and length of one periodic unit (cycle)= 1,992 μm)

Full size image

Now the procedure has become clear, we show next the mixing analyzes results in detail for various SHM designs.

3.4 Validation of the mapping method

To validate the method, we compare the mixing results obtained by mapping with those of numerical mixing results of Kang and Kwon (2004) and experimental mixing patterns of Stroock et al. (2002). To diminish numerical diffusion, and to simplify comparison with results reported in Kang and Kwon (2004), in this sub-section we compute only one mapping matrix \({\varvec{\Upphi}}\) representative of one whole cycle (repeating unit). The concentration distributions after each cycle are given by \({\user2{C}}^1= {\varvec{\Upphi}}{\user2{C}}^0, {\user2{C}}^2= {\varvec{\Upphi}}{\user2{C}}^1,\) etc. Figure 5 depicts the evolution of mixing for five cycles, comparing the results from mapping with the numerical results of Kang and Kwon (2004) and the experimental results of Stroock et al. (2002). Excellent agreement is found.

Fig. 5
figure 5

Evolution of the concentration distribution in a SHM during 5 total cycles a mapping results, b experimental results (confocal micrographs) from Stroock et al. (2002), and c numerical results from Kang and Kwon (2004). The geometrical features used are: h = 77 μm, w = 200 μm and g d  = 17.7 μm, g w  = 70.7 μm, and length of one periodic unit (cycle)= 1,992 μm

Full size image

3.5 Effect of groove depth and groove number on mixing

It is known that mixing in the SHM is sensitive to the groove depth and even an increase or decrease by 10% can affect mixing (Aubin et al. 2005; Bennet and Wiggins 2003). Here, five groove depths, g d = 10, 20, 30, 40, and 50 μm, are chosen and the Reynolds number is fixed to be 0.004 for all the five cases. To obtain concentration evolutions, for each groove depth a mapping matrix is computed, representative of the repeating unit. Figure 6 shows the results of mixing and reveals that deeper grooves provide better mixing, especially at the initial stage of mixing. To quantify mixing, we compute the flux-weighed intensity of segregation vs. the pressure drop up to 20 mixing cycles. The pressure drops per cycle for g d = 10, 20, 30, 40, and 50 μm are found to be 675, 660, 654, 651, and 649 Pa (N/m2), respectively, revealing that as groove depth increases the pressure drop per cycle slightly decreases. Figure 7 shows the logarithm of intensity of segregation vs. pressure drop for the five SHMs and, in accordance with the literature (Aubin et al. 2005; Bennet and Wiggins 2003), the fastest decrease in intensity of segregation (for a given pressure drop) is found for the 50 μm deep groove, but the 20, 30, and 40 μm depths are quite equivalent to 50 μm, especially at longer mixer lengths (or pressure drops). At g d  = 10 μm, mixing is the poorest even at longer mixer lengths, indicating that a minimum groove depth is necessary to induce chaotic advection. We conclude that the groove depth should exceed a critical value of roughly 20% of the channel depth to induce chaotic mixing. Once this critical value is exceeded, the effect of groove depth on mixing becomes insignificant, although, in the initial stages of mixing, slightly better mixing for the deeper grooves is achieved.

Fig. 6
figure 6

Effect of groove depth on the evolution of the concentration distribution in a SHM for designs with grooves depths of a 10, b 20, c 30, d 40, and e 50 μm, respectively (Channel depth is 77 μm; channel width 200 μm and length of one periodic unit (cycle) 1,992 μm)

Full size image
Fig. 7
figure 7

The intensity of segregation versus pressure drop (up to 20th mixing cycle) for designs with groove depths of 10, 20, 30, 40, and 50 μm

Full size image

Apart from groove depth, the number of grooves per half cycle is a design parameter that can be optimized. Periodic alternation of groove patterns in the SHM after each half cycle results in crossing streamlines and chaotic mixing. In general there is an optimum interval such that the total length stretch is maximum for a fixed length (spatially periodic flows) or time (time-periodic flows) of mixing. For example, in time-periodic 2-D cavity flows, Ottino and co-workers (Ottino 1989) found an optimal time period which maximizes mixing. In the SHM, optimization concerns the amount of stretching during the flow in half a cycle.

We used the geometrical parameters as mentioned in Sect. 3.1 to compute the mapping matrices for four mixing modules (see Fig. 2) as described in Sect. 3.2. Using these matrices, the concentration evolution and the corresponding intensity of segregation in SHMs with different number of grooves per half cycle is computed by the technique described in Sects. 3.2 and 3.3, respectively. To calculate the pressure drop per cycle, we first calculate the pressure drop across a single left groove (ΔP L ), using the fully developed velocity field as used for the calculation of \({\varvec{\Upphi}}_{\rm L},\) and the pressure drop across the left transition region (ΔP tLR), as used for the computation of \({\varvec{\Upphi}}_{t{\rm LR}}.\) Next, ΔP R is found equal to ΔP L, and ΔP tRL equal to ΔP tLR. By adding the individual contributions, the total pressure drop for the whole series of micromixers with different number of grooves per half cycle is obtained. Next, to find the optimum design, we plot the logarithm of the discrete intensity of segregation versus pressure drop in Fig. 8, which shows the results for SHMs with 6, 8, 10, and 14 numbers of grooves. From this plot, the smallest value of the intensity of segregation at a given pressure drop (a vertical line in Fig. 8) is found for 10 grooves per half cycle. Alternatively, a given mixing quality of e.g., log10 I d = −2.5 (a horizontal line in Fig. 8) is obtained with pressure drops close to 10.4, 9.5, 8.3, and 10.7 KN/m2 for SHMs with 6, 8, 10, and 14 grooves per half cycle, respectively, yielding the same optimum of 10 grooves (lowest energy used). This analysis reveals that for a fixed value of the transverse to axial velocity (which is decided by the geometrical parameters of a given SHM), there exists a minimum number of grooves per half cycle where mixing is optimum, and for the analyzed geometry this minimum proves to be 10.

Fig. 8
figure 8

The intensity of segregation versus pressure drop (up to 20th mixing cycle) for SHM designs with 6, 8, 10, and 14 grooves per half cycle. Clearly the optimum number of grooves is 10 given the lowest pressure drop (energy used) for the same mixing quality

Full size image

From the above analyzes, we conclude that various designs of the SHM can be analyzed using the mapping method in an efficient way and that the optimum design can be found.

4 Application to the barrier-embedded micromixer

Various combinations of two mixing protocols (functional modules) of the BEM provide numerous competitive designs. In this section, we demonstrate that the mapping method can be used as an efficient tool to analyze various layouts of these micromixer.

4.1 Mixer geometry

Figure 9a shows one periodic unit of the BEM with a barrier on the top-mid surface of a rectangular channel. The BEM can be thought of composed of two repeating units, say protocol P1 and P2, as shown in Fig. 9b. The first repeating unit P1 is a simple rectangular channel with six slanted grooves on the bottom surface inducing an overall rotational flow and the second repeating unit P2 has the same channel geometry as P1 except for a barrier on top-mid surface inducing two co-rotating flows. We choose the same geometrical features as used in the SHM (g d = 17.7 μm, g w = 70.7 μm, h = 77 μm, and w = 200 μm). The barrier height is 2/3 h, its thickness is 50 μm, and its length corresponds to the length of six grooves. To solve the periodic velocity field for P1 and P2, the same boundary conditions, and material properties as used for the SHM mentioned in Sect. 3.1 are applied. The grid used for P1 consists of 469,800 hexahedral cells and 498,318 nodes and for P2 it consists of 421,400 hexahedral cells and 454,518 nodes.

Fig. 9
figure 9

BEM. a A typical periodic unit of a barrier embedded mixer (BEM) with 12 grooves on the bottom surface and the barrier on the top. The barrier length corresponds to that of six grooves, b schematics of mixing protocols of BEM, P1 and P2, seen from the top and front. The first protocol, P1, is a rectangular channel with slanted grooves, while the second protocol, P2, consists of a barrier located exactly at mid of the top surface. The gray and black areas represent grooves and a barrier, respectively

Full size image

4.2 Defining mapping matrices for the BEM

Different ordering of the mixing protocols of the BEM, P1 and P2, results in different designs of the micromixer. Two mapping matrices \({\varvec{\Upphi}}_{1}\;\hbox{and}\;{\varvec{\Upphi}}_{2}\) are computed, representative for the mixing protocols P1 and P2. Next, any specific design with a periodic or aperiodic sequence can be evaluated using sequential multiplication of respective matrices (as specified in a given sequence of the two protocols) with a prescribed concentration vector. The number of tracers per cell and the number of mapping cells used to compute mapping matrices are same as for the SHM (see section 3.2 and schematically Fig. 1).

4.3 Mixing analysis in the BEM

Only three representative designs are considered as illustrative examples. The first design is periodic sequence composed of only P1, the second one is the periodic alternation of P1 and P2, and the third is an aperiodic (random) sequence of P1 and P2, as used by Liu et al. (1994) in cavity flows:

$$ P11: {{\mathbf{1\;1\;1\;1\;1\;1}}} \cdots {{\mathbf{1\;1\;1\;1\;1\;1}}}, $$
(17)
$$ P12: {{\mathbf{1\;2\;1\;2\;1\;2}}} \cdots {{\mathbf{1\;2\;1\;2\;1\;2}}}, $$
(18)
$$ AP12: \underbrace{{\mathbf{12}}} \underbrace{{\mathbf{21}}} \underbrace{{\mathbf{2112}}} \underbrace{{\mathbf{21121221}}} \underbrace{{\mathbf{2112122112212112}}}, $$
(19)

where a boldface number 1 represents the protocol P1, 2 the protocol P2, (see Fig. 9b). All the sequences consists of 30 mixing protocols. The sequence P11 is the SGM, P12 is the BEM, and AP12 is the variation of the BEM, which is composed of P1 and P2 in a recursive way.

Figure 10 shows the mixing evolution at down-channel positions z = 4, 10, 20, and 30L, where L is length of a mixing protocol (=1,078 μm). Figure 11 shows the quantitative characterization of mixing using the flux-weighted intensity of segregation. In the case of P11, the flow is only able to rotate fluid in monotonic way around the elliptic point without significant increase in mixing, while mixing in the protocol P12 is almost chaotic except for several tiny unmixed islands dispersed in whole cross-section (see Fig. 10). The aperiodic sequence AP12 also reveals chaotic mixing in most of the domain, far better than P11 from the viewpoint of mixing, but not better than P12. Further optimization of the designs, however, still seems to be possible, but it is beyond the scope of this paper.

Fig. 10
figure 10

Evolution of mixing patterns at several down-channel positions, z = 4, 10, 20, and 30 L, for the two periodic sequences and one aperiodic sequence. The initial concentration at z = 0 is shown on the top of a. a P11, b P12, and c AP12 (the geometrical features used are: h = 77 μm, w = 200 μm and g d  = 17.7 μm, g w  = 70.7 μm, and length of one periodic repeat unit = 1,078 μm)

Full size image
Fig. 11
figure 11

The intensity of segregation to quantify mixing in periodic and aperiodic sequences P11, P12, and AP12 composed of P1 and P2 protocols. Here, L is length of the one periodic repeat unit (=1,078 μm)

Full size image

5 Application to the three-dimensional serpentine channel

Finally, we consider the 3D-SC to show that the method can also work as an analysis tool in devices producing complicated flows. In this mixing device, due to the presence of back-flows, it is necessary to employ time-tracking rather than axial tracking, as adopted for the SHM and the BEM, to compute particle distributions at down-channel positions.

5.1 Mixer geometry

Figure 12 shows one periodic unit of the three-dimensional serpentine channel used in the experimental study of Liu et al. (2000). The basic building block is a “C-shaped” section. The geometrical features adopted are as follows: the inlet and outlet cross-sections are all 300  μm wide and 150 μm high, the length of C-shaped section is 900 μm and in total the length of the channel is 1,200 μm. Four Reynolds numbers, Re = 0.01, 10, 50, and 70, are considered to analyze dependence of mixing on inertia in the channel. To obtain the fully developed velocity field for one cycle of the channel at each Reynolds number, periodic boundary conditions with mass flow rate are prescribed to solve the Navier-Stokes equations. The grid used to obtain the periodic velocity field for the channel consists of 500,000 hexahedral cells and 529176 nodes.

Fig. 12
figure 12

Schematic of periodic unit of three-dimensional serpentine channel consisting of C-shaped building block. The inlet and outlet cross-sections are all 300 μm wide and 150 μm deep. The length of C-shaped section is 900 μm and in total the length of the channel is 1,200 μm

Full size image

5.2 Defining mapping matrix for the 3D-SC

Since the fluid flow at higher Reynolds numbers induces back-flows (u z  < 0), while in Stokes flow regime the axial velocity approaches to zero when fluid moves in perpendicular to the axial direction (z-direction), the axial integration Eq. (10) to track tracer positions fails. Therefore, to compute the coefficients of the mapping matrix \({\varvec{\Upphi}}\) (see Eq. (2)), we apply time integration to find the position of the tracers:

$$ \frac{\hbox{d}x}{\hbox{d}t} = u_{x}, \quad \frac{\hbox{d}y}{\hbox{d}t} = u_{y},\quad \frac{\hbox{d}z}{\hbox{d}t} = u_{z}. $$
(20)

Time-tracking approach is computationally more expensive as compared to axial tracking. This is due to the two facts. First, instead of solving two equations as in axial tracking (see Eq. (10)), one is required to solve three equations. Second, the large variations in axial velocity in the channel brings different residence time distributions of tracers according to their cross-sectional positions and, hence, tracers reach to final position at different times. A tracer close to walls requires much more time steps than one in the center. Therefore, to track a huge number of tracers using time integration (see Eq. (20)) to the end period of the mixer is cumbersome. However, using the mapping method is a better option since it requires time integrations to be performed only once for the representative repeat unit.

5.3 Mixing analysis in the 3D-SC

We study the progress of mixing for four Reynolds numbers, Re =  0.01, 10, 50, and 70. Liu et al. (2000) demonstrated experimentally that chaotic mixing can be achieved in this device, once inertia is significant.

Figure 13 shows the mixing evolutions along the down-channel positions after 1, 2, 3, 4, and 10 mixing cycles for Re = 0.01, 10, 50, and 70. As the Reynolds number increases, stretching and folding of interfaces becomes more vigorous and it is evident that the flow at Re = 70 is capable of producing the best mixing, while in Stokes flow regime (Re = 0.01) the flow is totally incapable to mix the two fluids. The onset of chaos in the whole cross-section requires a minimum cross-over Reynolds number. Figure 14 shows the quantitative comparison of mixing, using the flux-weighted intensity of segregation for the four Reynolds numbers. Higher Reynolds numbers provide better mixing, and between Re = 10 and Re = 50 we detect a change from regular to chaotic mixing. The mixing rate at the highest Reynolds number is best, but of course requires more energy.

Fig. 13
figure 13

Effect of increasing Reynolds number on the evolution of mixing patterns in the serpentine channel at several down-channel positions after 1, 2, 3, 4, and 10 cycles of mixing for a Re = 0.01, b Re = 10, c Re = 50, and d Re = 70. The initial concentration at z = 0 is shown on the top of a, and the length of one cycle is 1,200 μm

Full size image
Fig. 14
figure 14

The intensity of segregation plot up to the 10 cycles of mixing in the serpentine channel for four Reynolds numbers Re = 0.01, 10, 50, and 70

Full size image

This analysis shows that the mapping method is capable of analyzing mixing in quite complicated types of flow.

6 Conclusions

In this paper, we developed a new approach to compute the distribution matrices in a mapping method, and showed the importance of the method as an efficient toolbox to analyze, design, and optimize micromixers. The concentration distribution evolution can be computed and quantified via the use of the discrete flux-weighted intensity of segregation. To show the capabilities of the method, we investigated mixing in three well known micromixers: the SHM, the BEM, and the 3D-SC.

To map the SHM, the geometry was subdivided into four functional modules, for each of which a mapping matrix was computed. Different combinations of matrices result in concentration evolutions and the corresponding mixing measures for different lengths and lay-outs. This was the basis for the optimization of an important design parameter like the groove number per half cycle. In the SHM investigated (with channel height 77 μm and channel width of 200 μm), a groove depth of 10 μm is clearly insufficient to induce chaotic advection, while for the other four depths investigated, g d  = 20, 30, 40 and 50 μm, the deepest one gives the best mixing performance. However, at a sufficiently long distance from the inlet, differences among these four groove depths disappear. The optimum number of grooves per half cycle proves to be 10. Various designs of the BEM can be realized by combining different sequences of the two mixing protocols. The mixing analysis by the multiplication of respective matrices with the specified concentration at the inlet proves to be a very efficient way to predict the best possible design. This example, and the SHM analyzes, showed that the mapping method can work as an engineering design tool to find out an optimal design for numerous micromixers.

As for the 3D-SC, the flow characteristics are different, as compared with the above two examples, due to the presence of back-flows. This requires equation of motions to be integrated with respect to time, which is time consuming. However, mapping method requires this computationally expensive tracking step to perform only once for a repetitive unit. The outcome of mixing analysis in the channel using this approach indicates that the mixing quality is highly dependent on inertia, and as inertia increases mixing improves. In the analyzed range of Reynolds numbers, only flows at Re = 50 and 70 induce chaotic mixing, while flow at Re = 10 induces regular mixing, and at Re = 0.01 produces no mixing at all.