Optimized Simulation and Validation of Particle Advection in Asymmetric Staggered Herringbone Type Micromixers

Abstract

This paper presents and compares two different strategies in the numerical simulation of passive microfluidic mixers based on chaotic advection. In addition to flow velocity field calculations, concentration distributions of molecules and trajectories of microscale particles were determined and compared to evaluate the performance of the applied modeling approaches in the proposed geometries. A staggered herringbone type micromixer (SHM) was selected and studied in order to demonstrate finite element modeling issues. The selected microstructures were fabricated by a soft lithography technique, utilizing multilayer SU-8 epoxy-based photoresist as a molding replica for polydimethylsiloxane (PDMS) casting. The mixing processes in the microfluidic systems were characterized by applying molecular and particle (cell) solutions and adequate microscopic visualization techniques. We proved that modeling of the molecular concentration field is more costly, in regards to computational time, than the particle trajectory based method. However, both approaches showed adequate qualitative agreement with the experimental results.

1. Introduction

The precisely controlled manipulation of fluids in microanalytical systems or microchemical reactors is a key issue in terms of the use of these devices. In addition to fluidic transport, integrated functional microfluidic elements (pumps, mixers, separators, etc.) are also essential building blocks of sample preparation systems. The reliable modeling of the microscale fluidic processes in these systems is of critical importance to their economical design and development. Accordingly, the finite-element modeling of these microstructures is of current interest, with respect to their use in analytical devices [1,2]. The reliable modeling of biological fluids, with a precise estimation of particle motion, is a significant challenge due to the non-Newtonian behavior of the media and the special geometric and physical properties of the cells [3,4]. Although fast and cost effective, the prediction of the functional performance of preliminary analysis of the increasingly complex microfluidic systems is problematic considering that a coarse mesh resolution can deteriorate the resulting solutions. Our aim was to compare different modeling strategies to suggest a truly economical way to simulate demonstrative microfluidic systems, and to find a compromise between detailed solutions and fast analyses.

One of the basic functions of microfluidic systems integrable into microanalytical devices such as microchemical reactors is the dilution and complete mixing of the analyte with an adequate buffer solution or different reactants to ensure homogeneous concentration distribution over the critical active area of the architecture [5,6]. However, the mixing possibilities are highly limited at the microscale since turbulent flow cannot build up due to the dominant viscous forces. Evolving transversal (or chaotic) advection [7] could be a promising and effective mixing method in case of microfluidic applications considering stable and laminar flow in low Reynolds regime [2,8]. We have chosen chaotic advection as a complex process to demonstrate the performance of the different modeling approaches. Micromixers could be the key components in sample preparation units of proposed lab-on-a-chip systems providing efficient mixing of sample and analyte fluids or reagents in spite of the laminar flow present at the microscale.

1.1. Mixing at the Microscale

Mixing can be explained as a gradient driven transport process increasing entropy and decreasing global inhomogeneity of particle concentration, temperature and phases. Various effective mixing strategies were demonstrated at the macroscale such as molecular diffusion, turbulent diffusion, advection and Taylor-dispersion. In contrast to macroscale fluid dynamics where turbulent flows are easily achieved, microscale flow processes cannot generate turbulence due to the dominant viscosity of the system. Thus, the efficient mixing possibilities are limited. This is because the small characteristic dimensions of the microfluidic systems liquid flows are laminar and characterized by low Reynolds number. In such a laminar fluidic environment, the species can be transported by molecular diffusion between the component streams, where a dynamically diffusive interface is created with predictable geometry. Advection is defined as a transport process generated by the fluid flow causing a nonlinear or even chaotic distribution of molecules. This phenomenon can be realized using peculiar geometries. The Taylor-dispersion can be explained as a special advection generated by the sheared flow with mixing efficiency two orders of magnitude greater than strategies based on molecular diffusion [9].

1.2. Chaotic Advection

The transversal or chaotic advection refers to a transport process where advection, i.e., particle transport is generated by the fluid flow. Even in case of a simple laminar flow, velocity distribution can lead to a chaotic motion of particles without the need for turbulence. In chaotic advection, the flow velocity components are constant over space and time in a steady flow whereas in a turbulent system these are considered to be random. In a three-dimensional steady state advection flow, the representative streamlines of the system or the particle trajectories in a particle tracing setup can cross each other, thus fulfilling the requirements of minimum complexity for chaotic behavior [10]. Advection can become chaotic in a two-dimensional unsteady flow (unsteady meaning time-dependent perturbations), and in three-dimensional steady or unsteady flows [11]. The importance of chaotic advection, therefore, lies in the possibility to enhance mixing in the laminar flow regimes by orders of magnitude compared to molecular diffusion, thus enabling the use of smaller microfluidic devices with slightly more complex channel geometries (such as staggered herringbone type double layered channels) to greatly outperform straight channels in mixing and sample preparation tasks known to be of critical importance at the microscale.

2. Experimental Section

2.1. Applied Geometries

To demonstrate and compare different modeling strategies for mixing processes, various staggered herringbone microstructures were designed with different parameters. The microchannel cross-section dimension was set to 100 μm (height) × 200 μm (width) with the groove depth of 100 μm (Figure 1). The grooves are patterned at a 45° angle to the x axis asymmetrically to divide the channel into one third and two thirds as suggested by Stroock et al. [7]. One mixing unit was defined as two sets of consecutive herringbone grooves with two different cusp positions as presented in Figure 1. The grooves comprising the same half unit are intended to build transversal advection from the same side of the channel via the generation of counter-rotating vortices. The total length of the mixer channel was fixed at 6000 μm. For detailed modeling and experimental verification of the advective processes and mixing performances, different microfluidic geometries were applied using the following parameters: width of the grooves varied between 30 and 40 μm, the number of grooves per half unit varied between 4 and 6 and the number of mixing units varied between 6 and 4 as summarized in Table 1.

Figure 1. Parameters of the Staggered Herringbone Mixer (SHM): (a) width (W) and height (H) of the channel, depth (h) and width of the herringbones (a), and the schematic structure of the overall microfluidic system (b) including inlets and outlet and one mixing unit.

Figure 1. Parameters of the Staggered Herringbone Mixer (SHM): (a) width (W) and height (H) of the channel, depth (h) and width of the herringbones (a), and the schematic structure of the overall microfluidic system (b) including inlets and outlet and one mixing unit.

Table 1. Structural parameters of the applied micromixers. Width of the herringbone grooves, number of herringbones per half unit, and the number of mixing units were considered in our study.

Herringbone groove width	Number of herringbones per half unit	Number of mixing units	Name
30 μm	4	6	30/4/6
30 μm	6	4	30/6/4
35 μm	4	5	35/4/5
35 μm	5	4	35/5/4
40 μm	4	5	40/4/5
40 μm	5	4	40/5/4

Table 1. Structural parameters of the applied micromixers. Width of the herringbone grooves, number of herringbones per half unit, and the number of mixing units were considered in our study.

Herringbone groove width	Number of herringbones per half unit	Number of mixing units	Name
30 μm	4	6	30/4/6
30 μm	6	4	30/6/4
35 μm	4	5	35/4/5
35 μm	5	4	35/5/4
40 μm	4	5	40/4/5
40 μm	5	4	40/5/4

2.2. Modeling

This paper focusses on two different approaches regarding the modeling of microscale mixing. The concentration-based method calculates the diffusion of a reagent and estimates the concentration field as a result. The particle tracing-based method calculates particle trajectories with information about the path of each particle along the channel. Both aspects are based on a pre-calculated stationary velocity field as a result of solving the Navier-Stokes equation [2] numerically with finite element method. Models were built in COMSOL Multiphysics version 4.4 [12]. Laminar inflow boundary condition was applied on the inlet with 0.002 m/s average velocity in case of each numerical study of this series and zero pressure was set at the outlet. On the channel walls, no slip boundary conditions were defined for the laminar flow model as summarized in Table 2. The properties of room temperature water (density: 1000 kg/m³, kinematic viscosity: 10⁻⁶ m²/s) were applied as material parameters. The maximal 0.0112 Reynolds number was calculated by the simulator solver considering 0.002 m/s average inlet velocity. The average Reynolds number was estimated to be 0.00244.

Table 2. Summary of boundary conditions for the model.

Boundary	Model	Boundary condition	Value
Inlet	CFD	Laminar inflow with average flow velocity	0.002 m/s
Inlet	Trajectory	Particle inlet	6000 particles with uniform density
Inflow1	Concentration	Concentration	100 mol/m3
Inflow2	Concentration	Concentration	0 mol/m3
Channel wall	CFD	No slip	-
Channel wall	Trajectory	Bounce (for the trajectory model)	-
Channel wall	Concentration	No flux	-

Table 2. Summary of boundary conditions for the model.

Boundary	Model	Boundary condition	Value
Inlet	CFD	Laminar inflow with average flow velocity	0.002 m/s
Inlet	Trajectory	Particle inlet	6000 particles with uniform density
Inflow1	Concentration	Concentration	100 mol/m3
Inflow2	Concentration	Concentration	0 mol/m3
Channel wall	CFD	No slip	-
Channel wall	Trajectory	Bounce (for the trajectory model)	-
Channel wall	Concentration	No flux	-

In case of the concentration-based modeling, the concentration field of the diluted molecular size particles can be described by stationary scalar transport. In the finite element model, a step function initial concentration distribution was defined at the inlet of the mixer channel according to the influx of the particle solution. The concentration values at the mesh nodes were set to 100 mol/m³ on the left side and zero on the right side of the inlet section. At the channel walls and the outlet, no flux and outflow boundary conditions were defined, respectively. Diffusion coefficient of the model molecule was set to 10⁻¹⁰ m²/s. The mesh resolution was refined to avoid numerical diffusion effects [13], which could be a key issue in modeling molecular diffusion processes and typically induces false mixing effects in the channel at low mesh resolutions (see Table 3 for more details).

The particle tracing based model calculates, follows and depicts the individual particle trajectories according to the hydrodynamic drag force described by Stokes’ law. The lift force was not considered due to the low Reynolds number and the small difference between the particle and water density. For recording the individual particle trajectories a total number of 6000 spherical particles were released from the left half of the inlet plain of the mixer channel after the steady-state velocity field had been obtained. Properties of the model particles were set to be in correspondence with the principal properties of red blood cells (density: 1100 kg/m³, particle diameter: 6 μm [14,15]). Spherical particle geometry was used as an approximation of the cell geometry and this approach was in accordance with the experimental methods, as well. Particle trajectories follow the fluid streamlines established by prior hydrodynamic simulations, thus reducing the need for higher mesh resolution and higher computational demand as compared to the concentration-based model.

Tetrahedral mesh was applied for both the convergence study and the Finite Element Modeling of the microfluidic mixing processes in the different geometries described in Table 1, as well. Mesh element size for the convergence study is listed in Table 3. For the resulting concentration calculations, we used approximately 9,000,000 elements and for the particle trajectory calculations we used approximately 1,000,000 mesh elements.

Table 3. Mesh parameters for the convergence study on the 40/5/4 geometry.

Mesh name	Number of elements	Average element size (μm)	Standard deviation of element size (μm)
Fine	11,828,167	2.3752	1.8527
Medium	4,922,845	3.7606	2.4241
Low	1,788,920	5.4645	3.2021
Coarse	121,404	13.8776	8.5472

Table 3. Mesh parameters for the convergence study on the 40/5/4 geometry.

Mesh name	Number of elements	Average element size (μm)	Standard deviation of element size (μm)
Fine	11,828,167	2.3752	1.8527
Medium	4,922,845	3.7606	2.4241
Low	1,788,920	5.4645	3.2021
Coarse	121,404	13.8776	8.5472

2.3. Fabrication

Polydimethylsiloxane (PDMS) [16] soft lithography is a conventionally applied fabrication method for polymer microfluidic structures. Microfabrication into PDMS needs moulding and polymerising on the photoresist replica, then removal by peel-off from the moulding form as schematically demonstrated in Figure 2. Micropatterned photoresist structures such as SU-8 are ideally suitable for such moulding forms. For our proposed mixer structure, a multi-layered SU-8 epoxy based negative photoresist was used [17,18,19] as a master replica for PDMS molding. Lin et al. reported substrate fabrication for 3D herringbone mixer by femtolaser direct writing [20] as an alternative fabrication method. An improved 3D SU-8 multilayer fabrication process was developed to be applicable for formation of microfluidic systems with high aspect ratio sidewalls and advanced functional parts. SU-8 layers were patterned by subsequent spin-coatings (Brewer Science Cee 200CBX spin-coater [21]), lithographic exposures (Süss MicroTech MA6 mask aligner [22]) and a final development step, as illustrated in Figure 2a. The layers of the SU-8 replica contained different functional components of the micromixer. The main microfluidic channel and reservoirs were fabricated in the first layer while the inlets and outlets as well as the herringbone structures were created in the second layer. PDMS prepolymer was poured onto the developed replica and polymerized in two days in room conditions. At the final step of the fabrication process, PDMS was sealed to glass by low temperature bonding after oxygen plasma treatment applying 200 W plasma power, 100 kPa chamber pressure and 1400–1900 sccm oxygen flow (Terra Universal Plasma Preen Cleaner/Etcher [23]). Hence, stable covalent bonds were developed between the surfaces.

These materials have excellent structural, mechanical, optical and technological properties which make them suitable for application in microtechnology [24]. PDMS is a silicon-based organic polymer: (H₃C)₃SiO[Si(CH₃)₂O]_nSi(CH₃)₃ which has become a versatile material to realize microfluidic test structures due to its transparency, flexibility, reliable geometry transfer, price, and biocompatibility. This material provides an excellent solution for fast prototyping of simple microfluidic test structures, and can also be convenient for more complex demonstrator applications, although its applicability for large scale production is not yet proven considering its long term chemical and biological stability.

Figure 2. (a) SU-8 multilayer lithography steps and (b) the 3D SU-8 molding replica imaged by Scanning Electron Microscopy (SEM).

Figure 2. (a) SU-8 multilayer lithography steps and (b) the 3D SU-8 molding replica imaged by Scanning Electron Microscopy (SEM).

2.4. Measurement

To validate the model results, both modeling strategies were tested in the microfabricated fluidic systems. The concentration-based mixing process was characterized by simultaneous injection of diluted yellow and blue food color dye (molecular weight is approximately 500 g/mol) [25], and the mixing of the two solutions was observed. We used an upright microscope (Zeiss AxioScope A1 [26]) with bright field color imaging applying adequate flow parameters corresponding to the modeled values (4 × 10⁻¹¹ m³/s inlet flow rate).

To validate the particle-based modeling method, yeast cells were diluted in room temperature, with the water being of a similar volume as set in the simulation. The size distribution of the applied yeast cells was verified by Diatron Abacus Junior 30 hematology analyzer [27] and compared to the size histogram of human red blood cells (RBC) as presented in Figure 3. The RBC concentration was around 4.2–5.4 × 10¹² L⁻¹ (woman) and the applied yeast cell concentration was 7.12 × 10¹⁰ L⁻¹ to avoid clogging and multiple scattering during the measurements. The hematology analyzer uses 25 μL sample volume for one test. The spherical shape used in the model calculations is a good approximation for the yeast cells. To enhance and visualize cell trajectories and the resulting mixing states in the microchannel, dark field imaging was used. This imaging method facilitates the recording of the light scattered from the cells crossing the light beam, and we estimated the local cell concentration from the lateral distribution of the scattered light intensities, i.e. from the local brightness levels of the image. A relative light intensity was calculated by transforming the recorded intensities to follow the following conservation law:

$$\frac{{\displaystyle {\int }_0^{w}I\left(x_0,y\right)}dy}{w}=cont.=\frac{{\displaystyle {\int }_0^{w}I\left(0,y\right)}dy}{w}=\frac{{\displaystyle {\int }_0^{w/2}1}dy+{\displaystyle {\int }_{w/2}^{w}0}dy}{w}=\frac{1}{2}$$

(1)

where x₀ and y represents the local coordinates in the channel at a given cross-section and w denotes the width of the channel. This calibration ensures that the integral of the relative scattered light intensity is constant and equal to the initial value (0.5 at the inlet).

Figure 3. Size distribution of yeast cells and red blood cells. The two cell types have a similar range of cell diameters making yeast cells a well suited model of human cells in our experiments. The significant peaks in the volume distribution of the RBCs and the yeast cells are around 100 fL, and the spherical equivalent diameter belonging to this value is 5.76 μm. In the size distribution of the yeast cell, there is another peak at around 30 fL with an equivalent diameter of 4.16 μm. The total number of counted yeast and human red blood cells was 1.78 × 10⁶ and 1.25 × 10⁸, respectively.

Figure 3. Size distribution of yeast cells and red blood cells. The two cell types have a similar range of cell diameters making yeast cells a well suited model of human cells in our experiments. The significant peaks in the volume distribution of the RBCs and the yeast cells are around 100 fL, and the spherical equivalent diameter belonging to this value is 5.76 μm. In the size distribution of the yeast cell, there is another peak at around 30 fL with an equivalent diameter of 4.16 μm. The total number of counted yeast and human red blood cells was 1.78 × 10⁶ and 1.25 × 10⁸, respectively.

3. Results and Discussion

3.1. Modeling Results

At first, the laminar flow field was modeled as a basis for estimation of both the diffusion-based transport of diluted molecular species and the advection-based particle trajectories. Detailed analysis and visualization of the velocity field in Figure 4 highlights the rotating streamlines due to anisotropic fluid resistance generated by the secondary channel grooves.

The applied mesh resolution is a key issue in concentration-based modeling of micromixers as proven by our work. The numerical diffusion demonstrated a significant role during modeling the mixing process with lower mesh resolutions (mainly above 5 μm mesh element size). Figure 5 clearly shows the effect of numerical diffusion in the case of the herringbone micromixer after one mixing unit. In this case, the numerical diffusion may induce significant error in the estimation of the resulting concentration field as demonstrated by the vanishing phase stratification in the lowest resolution numerical solution (Figure 5d).

Mesh convergence study was implemented for the trajectory model as well and the results were summarized in Figure 6. As a demonstrative model parameter, the mixing efficiencies were calculated in case of different mesh resolutions in the 40/5/4 geometry after one mixing unit considering the ratio of the cells on the right side and the left side of the channel. A slight and saturating improvement in the modeling results could be experienced, but the effect was significantly smaller compared to the diffusion model.

To observe the evolution of the mixing states four planes of interest were defined locally: one at the inlet, one after the first and one after the second mixing unit and one at the outlet (Figure 7). Concentration fields in these planes (Figure 7a–d) show the mixing effect of the rotation considering the additional diffusion that occurs on the increased interface between the two fluids due to transversal advection.

Figure 4. Velocity field in the cross section of the channel and grooves. Coloring and arrows denote the amplitude (blue: left direction, yellow to red: right direction) and the direction of the y component of the local flow velocity vector, respectively. Rotating effect of the herringbone shaped grooves is clearly observable due to the flow direction modification effect of the grooves.

Figure 4. Velocity field in the cross section of the channel and grooves. Coloring and arrows denote the amplitude (blue: left direction, yellow to red: right direction) and the direction of the y component of the local flow velocity vector, respectively. Rotating effect of the herringbone shaped grooves is clearly observable due to the flow direction modification effect of the grooves.

Figure 5. Effect of mesh element size on the quality and reliability of the numerical solution in case of the concentration based model. Concentration fields after the first mixing unit with fine (a), medium (b), low (c) and coarse (d) mesh resolutions. The numerical diffusion caused by the poor mesh resolution is well observable on the more extensive green colored areas on (b) and (d) compared to the respective areas on (a) and (c).

Figure 5. Effect of mesh element size on the quality and reliability of the numerical solution in case of the concentration based model. Concentration fields after the first mixing unit with fine (a), medium (b), low (c) and coarse (d) mesh resolutions. The numerical diffusion caused by the poor mesh resolution is well observable on the more extensive green colored areas on (b) and (d) compared to the respective areas on (a) and (c).

Figure 6. Mesh convergence study for the trajectory model, by calculating the ratio of the number of transferred and remained particles. The difference between values is calculated from to the actual and the previous mesh resolution results (depicted with grey bars). Poincaré maps at the outlet for the extremely coarse and the fine mesh show minor differences in the particle trajectories compared to the diffusion model.

Figure 6. Mesh convergence study for the trajectory model, by calculating the ratio of the number of transferred and remained particles. The difference between values is calculated from to the actual and the previous mesh resolution results (depicted with grey bars). Poincaré maps at the outlet for the extremely coarse and the fine mesh show minor differences in the particle trajectories compared to the diffusion model.

Figure 7. Mixing evolution of two different solutions modeled by concentration based method (a–d) and particles modeled by trajectory based method (e–h) along the microchannel. The mixing distributions are qualitatively similar which is observable on the similar distributions of shapes on (a–d) compared to (e–h), respectively.

Figure 7. Mixing evolution of two different solutions modeled by concentration based method (a–d) and particles modeled by trajectory based method (e–h) along the microchannel. The mixing distributions are qualitatively similar which is observable on the similar distributions of shapes on (a–d) compared to (e–h), respectively.

For modeling cell (or particle) distribution in the studied micromixer the particle tracing module of COMSOL Multiphysics was applied. Poincaré maps (Figure 7e–h) show the actual location of individual particles at the selected planes. The different points correspond to the plain sections of trajectories of the monitored particles. The color of the particles denotes their initial location along the y axis (blue: left, red: right) of the channel to support visualization of the cell displacements relative to each other. As revealed in Figure 7h, stratified particle “shells” have evolved on the right side of the channel due to the transversal advection effect. On the left side of the channel, the particles also show a layered structure with the buffer solution.

To determine the effects of the selected geometric parameters of the proposed herringbone mixers (see parameters in Table 1), the mixing efficiencies were characterized and compared to each other in the case of trajectory-based models. By creating the outline of the Poincaré sections at the outlet plane of the channel, the change of its area estimates the number of transferred particles from the initial channel-side to the other. We defined a mixing efficiency to be better if it corresponds to a larger transferred particle shell area relative to the cross section of the channel. Efficiency values were calculated by considering the ratio of the number of particles on the left (initial) side of the channel and the number of particles on the right side of the channel. Figure 8 and Table 4 show that mixing effect improves with increased width of the grooves as well as by applying more mixing units rather than more grooves per unit. The effect of the depth of the grooves has been characterized previously by Fürjes et al. [2] and it was found that saturating efficiency can be achieved when the depth of the grooves are similar in range to the channel depth.

Figure 8. Outlines of the mixed particles at the outlet plane of the mixer channel. The color of the lines denote the width of the grooves, line type indicates the number of grooves. Shrinkage of the area with large number of particles at their initial side (left) and expansion of the area rich in particles on the right side is well observable, which results in the narrowing of the unmixed region in the middle of the channel.

Figure 8. Outlines of the mixed particles at the outlet plane of the mixer channel. The color of the lines denote the width of the grooves, line type indicates the number of grooves. Shrinkage of the area with large number of particles at their initial side (left) and expansion of the area rich in particles on the right side is well observable, which results in the narrowing of the unmixed region in the middle of the channel.

Table 4. Mixing efficiencies for the proposed geometries.

Geometry	40/5/4	40/4/5	35/5/4	35/4/5	30/4/6	30/6/4
Efficiency	0.4526	0.4513	0.3655	0.2998	0.2105	0.1984

Table 4. Mixing efficiencies for the proposed geometries.

Geometry	40/5/4	40/4/5	35/5/4	35/4/5	30/4/6	30/6/4
Efficiency	0.4526	0.4513	0.3655	0.2998	0.2105	0.1984

3.2. Experimental Validation

To validate the applied numerical simulation methods, micromixer structures were fabricated and characterized with bright field microscopy for the concentration-based and dark field microscopy for the particle-based approach. Figure 9 shows a qualitative comparison of the modeled and the measured concentration fields near the channel outlet after three mixing units. Blue and yellow regions indicate high concentration levels of the unmixed food color dies while green color denotes the mixed regions. The deducted relative concentration curves recorded along the line section of the mixer channel by ImageJ [28] represent the mixing efficiency of the herringbone mixer.

When mixing the yeast cell solution in our device, the layered structure of the flow is well observable in accordance with the particle trajectory model results (Figure 10). The cell depleted and cell enriched regions (dark and bright regions in the channel) correspond to the modeled particle advections that evolved in the micromixer structure.

Figure 9. Modeled and measured concentration fields in case of diffusion based approach. (a) Result of finite element concentration modeling. (b) Mixing of food dies in the microchannel monitored by bright field microscopy imaging. (c) The recorded histograms represent the calculated and measured relative concentrations along a line section perpendicular to the channel direction.

Figure 9. Modeled and measured concentration fields in case of diffusion based approach. (a) Result of finite element concentration modeling. (b) Mixing of food dies in the microchannel monitored by bright field microscopy imaging. (c) The recorded histograms represent the calculated and measured relative concentrations along a line section perpendicular to the channel direction.

Figure 10. Relative scattered light intensity distribution corresponding to the local yeast cell concentrations recorded by dark field microscopy. (a) The cell-depleted region is clearly observable and coinciding with the food color dye concentration field presented in Figure 9. (b) The histogram extracted from the intensity field along a line section perpendicular to the channel direction is in good agreement with the modeled Poincaré map evolved at the corresponding plane.

Figure 10. Relative scattered light intensity distribution corresponding to the local yeast cell concentrations recorded by dark field microscopy. (a) The cell-depleted region is clearly observable and coinciding with the food color dye concentration field presented in Figure 9. (b) The histogram extracted from the intensity field along a line section perpendicular to the channel direction is in good agreement with the modeled Poincaré map evolved at the corresponding plane.

Both modeling methods demonstrated the impressive characteristics of the herringbone micromixer. The multilayer streamline structure developed by the rotating effects and the evolving transversal advection due to the special microchannel geometry were clearly visible in both cases. The enhanced mixing by molecular diffusion on the additional dynamic surfaces was well captured, and concentration fields showed a good qualitative agreement between simulated and observed results, proving the reliability of these modeling approaches. The particle trajectory model provides accelerated qualitative information on the functional behaviour of the proposed microfluidic systems. However, for a precise description and prediction of concentration dependent processes (e.g., in microreactors), the application of high resolution concentration-based models is required. The results obtained with COMSOL Multiphysics provide a solid basis for further development and testing of more complex microfluidic systems with cost-efficient numerical modeling.

4. Conclusions

Numerical modeling of complex microfluidic systems is a cost-efficient way of developing and testing novel structures. However, models have to be verified and set up properly to obtain reliable results. In this study, we compared two significantly different modeling approaches and verified the results experimentally to characterize the applied strategies regarding their mesh resolution sensitivity, reliability, performance and computational demand. For this purpose, the staggered herringbone micromixer was selected due to its high mixing efficiency and complex 3D microstructure. The mixing process was numerically modeled applying both concentration-based and particle trajectory-based approaches.

The concentration-based method showed a high sensitivity to mesh resolution having a significant effect on false numerical diffusion in model results. Therefore, a high resolution mesh with high computational cost was needed to obtain accurate results. Application of the particle tracing approach resulted in a much lower computational cost (approximately 20 times lower) and makes an extensive geometric parameter optimization more feasible. Both concentration and particle trajectory-based models were evaluated by representative measurements, and adequate accordance in the resulted mixing states was experienced. Comparing the estimated mixing efficiencies of the various structures, we verified that the performance of the mixer can be improved by increasing the width of the grooves and by increasing the number of herringbones per mixing unit.

Cell and particle manipulation is a key issue in the development of Lab-on-a-chip applications (e.g., CTC (circulating tumor cell) separation and counting [29]). Therefore, the quality of model simulations is critical in the effective development of these new devices. Our test structure—the staggered herringbone micromixer—was integrated into a microfluidic sample preparation system developed for polymer-based photonic biosensor applications in the P3SENS project [30].