The Effect of Inertia on the Flow and Mixing Characteristics of a Chaotic Serpentine Mixer

As an extension of our previous study, the flow and mixing characteristics of a serpentine mixer in non-creeping flow conditions are investigated numerically. A periodic velocity field is obtained for each spatially periodic channel with the Reynolds number (Re) ranging from 0.1 to 70 and the channel aspect ratio from 0.25 to one. The flow kinematics is visualized by plotting the manifold of the deforming interface between two fluids. The progress of mixing affected by the Reynolds number and the channel geometry is characterized by a measure of mixing, the intensity of segregation, calculated using the concentration distribution. A mixer with a lower aspect ratio, which is a poor mixer in the creeping flow regime, turns out to be an efficient one above a threshold value of the Reynolds number, Re = 50. This is due to the combined effect of the enhanced rotational motion of fluid particles and back flows near the bends of the channel driven by inertia. As for a mixer with a higher aspect ratio, the intensity of segregation has its maximum around Re = 30, implying that inertia does not always have a positive influence on mixing in this mixer.


Introduction
In microfluidic systems that integrate laboratory functions on a single chip using extremely small liquid volumes in either biology or chemistry, mixing is one of the key tasks in determining the overall success of the system [1][2][3][4]. Under typical operating conditions, flow in the systems is laminar. Thus, one cannot expect turbulent mixing as used in macroscale mixing devices, and mixing primarily is dominated by molecular diffusion. Mixing via molecular diffusion across the interface between two fluids is very slow, which makes efficient mixing in microfluidic devices difficult to achieve [5,6]. The typical diffusion time scale t d is represented by t d ∼ l 2 /D, where l is the characteristic length and D the coefficient of molecular diffusion. To overcome such difficulties, chaotic advection [7][8][9][10], which has been successfully applied to the development of macroscale mixers working in the laminar flow regime, were also adopted in micromixers [11][12][13][14]. It is now well known that, by incorporating chaotic advection, one can achieve enhanced mixing, even in the creeping flow regime (with negligible inertia), which is also the case in microscale flows.
Depending on the way of manipulating working fluids, there are two types of micromixers, passive and active mixers. Passive micromixers make use of a specific three-dimensional structure to induce lateral motions of fluids in a periodic (or aperiodic) manner, resulting in chaotic mixing in a laminar flow [12][13][14][15][16][17]. A passive mixer utilizes no external actuation except for a pressure head or pump, which is used to drive flows into the microfluidic device. In active micromixers, however, active controls over the flow field are introduced to manipulate fluids to be mixed by using moving parts, varying pressure gradients, acoustic waves or electro-magnetic fields [18][19][20][21][22][23][24]. For more information on this topic, we refer to review articles on micromixers [2][3][4].
In this study, we are interested in passive mixing utilizing serpentine channel geometries [15,16,21,25]. Liu et al. [15] proposed a three-dimensional serpentine channel creating two concentric rotational flows moving in opposite directions. Since there is no notable mechanism leading to chaotic mixing, efficient mixing occurs only above a threshold Reynolds number in this channel geometry. As reported by Park and Kwon [26], a bend in the channel generates a rotational motion of fluids, where the rotational motion is directly affected by the shape of the bend. Inspired by the observed flow characteristics of various serpentine-type channels [15,16,21] and the principle of the linked twist maps (LTMs) [1,10,28], we designed a chaotic serpentine mixer (CSM) with two repeating functional units that create the two rotational flows with crossing streamlines. In our previous study [25], we investigated flow and mixing characteristics of the CSM in the creeping flow regime. The "L-shaped" bends in the CSM create helical flows, rotating either clockwise or counterclockwise, depending on its shape. The bypass channels play a role in offsetting the centers of the two rotational flows and splitting their flow path. A CSM with a properly chosen set of geometrical parameters is capable of achieving efficient chaotic mixing in the creeping flow regime with negligible inertia.
Although typical flow phenomena in micromixers occur in the creeping flow regime, sometimes, flow with inertia is needed to get a much higher throughput in a microfluidic system or to achieve better mixing with the aid of inertia [15,25]. Such a non-creeping flow may influence the flow and, therefore, the subsequent mixing performance of a micromixer. Kang et al. [17] showed that the mixing sequence showing the best mixing performance at a certain Reynolds number is not always superior to other mixing sequences at a different Reynolds number regime. Therefore, to ensure the usefulness of the CSM in a wide range of flow conditions, it is necessary to examine the effect of inertia on the flow and mixing characteristics of the CSM.
As an extension to our previous work [25], in this study, the effect of inertia on the fluid flow and the mixing performance of the CSM will be examined using a Reynolds number ranging from 0.1 to 70. As for the geometrical parameters, we will restrict our interest to the aspect ratio of the L-shaped channel, since it is known as the most critical geometrical parameter having influence on the flow characteristics and subsequent mixing performance [25,27]. The detailed flow characteristics and the evolution of the inter-fluid interface in flow are demonstrated with an emphasis on the effect of inertia and the aspect ratio of the channel. The progress of mixing in a mixer geometry is characterized using a measured of mixing, called the intensity of segregation.

Mixer Geometry
As is well known, repeated stretching and folding of fluid elements is the key mechanism leading to chaotic mixing. According to the theory of the linked twist maps (LTMs) [10], a necessary condition for chaos is the crossing of streamlines, which should occur at different times in a two-dimensional time-periodic flow. In the case of three-dimensional spatially periodic flows, two successive flow portraits generated by two mixing protocols should have intersecting streamlines, when projected onto the same plane. Figure 1 depicts the periodic unit for the chaotic serpentine mixer (CSM) that consists of two functional units, repeating periodically and constituting one spatial period of the mixer. Since the mixer has a simple two-layer structure, the channel can be fabricated by a mass-production technique, such as micro-injection molding or hot embossing. The two functional units are designed such that they split fluid elements and create two asymmetric counter-rotating flows with crossing streamlines in an overall sense. Therefore, the mixer can be thought of as a suitable fit within the linked twist maps (LTMs) framework. The streamline portraits are shown at the lower-right corner of the figure. The arrows in the streamlines indicate the direction of the rotational flows viewed from the exit of the channel. In the schematic drawing of the streamlines, dashed and solid lines indicate streamlines during the first and second half periods in the periodic unit, respectively. Here, W 1 is the width of the "L-shaped" channel, W 2 the width of the bypass channel and H the height of one layer of the channel. The width of the inlet W is W = W 1 + W 2 and the length of one period is L = 8W 1 . The origin of the coordinate system is located at the center of the inlet surface. In the figure, Ω, Γ w , Γ i and Γ o denote the entire fluid domain, the wall boundary, the inlet boundary and the outlet boundary, respectively.
In the previous study, it was found that the two geometrical parameters α and β are the key design variables with influences on flow kinematics and resulting mixing, where α = H/W 1 and β = W 2 /(W 1 + W 2 ). The first variable α is the aspect ratio of the "L-shaped" channel and β is the ratio of the width of the bypass channel to the whole width of the rectangular channel. Figure 1. A channel geometry of the chaotic serpentine mixer that consists of "L-shaped" structures and bypass channels, leading to two counter-rotating flow portraits. Overall streamline portraits are shown at the lower-right corner of the figure. The arrows represent the direction of the rotational flows viewed from the exit of the channel. Here, W 1 is the width of the "L-shaped" channel, W 2 the width of the bypass channel and H the height of one layer of the channel. The length of one period is L = 8W 1 .

Modeling and Numerical Methods
Since the CSM is composed of a periodically repeating mixing protocol, it is assumed that the velocity field of a repeating unit is also spatially periodic. As seen from its application, in previous studies [17,25,29,30], this proves to be a good assumption when investigating the progress of mixing in a micromixer composed of repeating mixing protocols, each with their own flow characteristics. Thus, the computational domain for the flow and mixing analysis is defined as one periodic unit of the mixer, as shown in Figure 1. In the mixing analysis, the velocity field in a periodic unit is repeatedly reused in order to analyze mixing in a micromixer to the period up to which one wants to analyze.
The periodic velocity field is obtained by solving a set of governing equations for the flow problem with proper boundary conditions and constraints. For the simplicity of analysis, the two fluids to be mixed are assumed to be incompressible and Newtonian with the same density and viscosity. The steady Naiver-Stokes equations and the continuity equation for an incompressible fluid are: where ρ is the density, u the velocity, p the pressure and µ the viscosity. The boundary conditions and constraint equations are given by: Equation (3) is the wall boundary condition and Equation (4) is the constraint equation for the flow rate at the inlet. Here, Q is the imposed flow rate at the inlet boundary Γ i and n the outward unit normal vector at the inlet boundary. The periodicity of the velocity is warranted by the constraints for the velocity, Equation (5), in which u i and u o represent the velocities at the inlet and the outlet nodal points (facing each other), respectively. To complete the flow problem, a reference pressurep, sayp = 0, is specified at one point at the outlet. The relative influence of viscous and inertial effects in fluid flow is measured by the Reynolds number Re, defined by: where U is the average velocity at the inlet surface, D h is the hydraulic diameter of the inlet, defined by , and ν is the kinematic viscosity of the fluid, defined by ν = µ/ρ. The computational domain is discretized into 27-node hexahedral elements to solve the flow problem. The typical element size is defined as having at least 20 elements along the channel thickness direction. The Galerkin finite element method and a Lagrange multiplier method are employed to solve the flow problem governed by Equations (1)-(5). We use tri-quadratic interpolation for the velocity and tri-linear interpolation for the pressure, the so-called Q 2 /Q 1 element. The constraints for the flow rate and the periodic velocity, Equations (4) and (5), are enforced using Lagrange multipliers in a similar way as described in a previous work [25]. The resulting matrix equation is solved using a parallel sparse direct solver, PARDISO [31], on a symmetric multiprocessing (SMP) machine with two octa-core processors.
For a given velocity solution of a mixer, an interface between two fluids at the inlet is tracked to gain insight into flow characteristics and to visualize the progress of mixing in a mixer at a given Reynolds number. In particle tracking, particles are assumed to be passive; there are no interactions, such as fluid-particle and particle-particle interactions, so that particles are just convected along the given velocity field of a micromixer. During the particle tracking procedure, one can identify the initial locations of the particles by their intrinsic color assigned in the beginning. The particle positions at the outlet can be used as inlet positions for the following period. The periodic velocity field is used repeatedly to track particles to the period up to which one wants to analyze. Given the velocity field in a micromixer, the position vector x of each tracer is updated by solving an ordinary differential equation: where t is the time and x 0 the initial position of a tracer. A fourth-order Runge-Kutta method is employed to numerically integrate Equation (7) with respect to time.
In the mixing analysis, we solve the steady advection-diffusion equation for the concentration c, given by: where u is the velocity obtained in the flow problem, N the diffusive flux due to the concentration gradient, defined by N = −D∇c, and D the diffusivity of the fluid. The boundary conditions for the equation are summarized as follows: where c j−1 is the concentration at the inlet of the j − th period and n · N is the flux normal to the boundaries, Γ w and Γ o . Assuming two fluids are introduced through a T-type inlet channel in the beginning, the concentration at inlet nodal points is either zero or one, depending on the species of the fluid; in this study, c 0 = 0 in the left half and c 0 = 1 in the right half. Once the concentration c j at the outlet of the j − th period (as part of the solution of the mixing problem) is obtained, c j is used as an inlet condition of the mixing problem for the next period j + 1. The relative ratio of the rate of convective transport to the rate of molecular diffusion driven by the concentration gradient is represented by the Péclet number (P e), given by: To discretize Equation (8), the streamline-upwind/Petrov-Galerkin (SUPG) method [32] is employed to prevent possible numerical oscillations in the concentration field, and tri-linear interpolation is used for the concentration c. The resulting matrix equation is solved using the same solution technique employed to solve the flow problem.

Flow Characteristics
Firstly, a Poincaré section is plotted to examine the dynamical systems of a half periodic unit composed of an L-shaped channel and a bypass channel. As a representative example, Figure 2 depicts the Poincaré section obtained using the periodic velocity field at Re = 0.1 for the half periodic unit of the channel geometry shown in Figure 1. Here, the aspect ratio of the L-shaped rectangular channel is α = 1, and the ratio of the width of the bypass channel to the whole width is β = 1/3.
From the Poincaré section, one can clearly observe an overall rotational motion of the fluid particles. The flow field has transverse components in addition to the down-channel component, resulting in the rotational motion. Due to this rotational motion, a big KAM (Kolmogorov-Arnold-Moser) boundary is present on the right side of the cross-section. If only the half periodic unit is repeated to form a micromixer, one cannot expect efficient mixing in the mixer, due to the presence of unmixed islands. As shown in Figure 1, however, two consecutive functional units exhibiting mirror symmetry (in this case, with respect to the yz-plane) are able to generate two asymmetric rotational flows. If the two units form a periodic unit, then the two flows form crossing streamlines, which causes the flows to possibly induce chaotic advection. We found that, in our previous study, globally chaotic mixing is achieved only at a properly chosen set of geometrical parameters in the creeping flow regime [25].
The most important feature of a micromixer is its flow characteristics, which ultimately determine the mixing performance. To investigate the flow characteristics of the CSM in the non-creeping flow regime, we solved the flow problem for the mixer by changing the value of α from 0.25 to one and the Reynolds number from 0.1 to 70. In addition, we fix the value of β to be 1/3, making the aspect ratio α the only geometrical variable in numerical simulations. In these figures, the velocity vectors are projected onto a plane at z = 0.25L, where L is the length of the channel and z represents the coordinate along the down-channel direction (see Figure 1 for the definition of the coordinate system). By comparing Figure 3 with Figure 4, we observe that the aspect ratio of the L-shaped channel has a significant influence on the cross-sectional velocity components. In the case of α = 0.25, two rotating flows with opposite directions of rotation are created as the Reynolds number increases, whereas in the case of α = 1, only one rotating flow is observed at Reynolds numbers larger than 10. Thus, the structure of the flow field varies with the Reynolds number and the aspect ratio of the channel. In both cases, the higher the Reynolds number, the stronger the rotational flow is in the cross-section.

Interface Tracking
The evolution of the interface between two fluids is visualized through the process of particle tracking. Figures 5 and 6 show the manifolds of deforming interfaces in the CSM with α = 0.25 and one at four Reynolds numbers, Re = 10, 30, 50 and 70. The figures illustrate the split, the rotation and the recombination of the interface with the fluid moving from right to left. Initially, the interface consists of 5000 passive tracers with assigned colors dependent on their vertical locations at the inlet. The interfacial line is tracked down to the outlet to show the deformation pattern, which is dependent on the Reynolds number. Figure 5 depicts the evolution of the interface in a single period of the CSM with α = 0.25 at the four Reynolds numbers. At the lowest Reynolds number (Re = 10), the flow portraits in the channel reveals a clockwise rotational motion, followed by a counterclockwise rotation (viewed from the exit). At the first bend, as shown in Figure 5, the vertical interface rotates during the half period. During the next half period, the interface experiences both a split and rotation, thereby increasing the number of interfaces at the exit. One can clearly see the topological change (due to split and rotation) of the interface through the colors being assigned to the interface. Due to repeating of this split and rotation, the flow characteristics are thought to be the mechanism inducing mixing. However, as the Reynolds number increases, the flow becomes more complicated, due to the enhanced rotational flows observed in the cross-sectional velocity plots. In addition to the enhanced rotational flow, one can observe stronger back flows, followed by rotational flows near the bends of the channel, that make the interface even more complex. At Re = 70, as depicted in Figure 5d, the interface at the outlet after one periodic unit is more complicated than that at Re = 10, due to the combined effect of back flows and rotational flows driven by inertia. One can expect more enhanced mixing would occur within a few periods at Re = 70, which needs to be confirmed in the mixing analysis that will be discussed shortly. Figure 6 shows the evolution of the interface for the CSM at α = 1. In this case, we observed a single rotational flow along with its strength increasing along with an increase in the Reynolds number. At higher Reynolds number flows, a twisted interface is formed due to the rotational flow. Similar back flows, as observed in the flow with α = 0.25, are also found near the bends. From the manifolds of the deforming interface, it is not clear whether or not inertia has a positive influence on the mixing as the Reynolds number increases. Thus, a detailed mixing analysis is required to elucidate the mixing performance of the mixer with the Reynolds number.

Mixing Analysis
Even though the manifold of the deforming interface in a specific micromixer can qualitatively help us visualize the progress of mixing, one cannot estimate the rate of mixing and the mixing state at the end of the micromixer, which are of great importance in a practical application. From the viewpoint of designing a micromixer, one should know how fast mixing would be and how long the length of the micromixer needs to be in order to achieve complete mixing. In this regard, we carried out mixing analyses for the CSM using the obtained velocity solutions for the Reynolds numbers ranging from 0.1 to 70. At a certain cross-section of a micromixer, the distribution of the concentration c can reflect the status of mixing, both qualitatively and quantitatively.
Concentration distributions, which are obtained as solutions of the advection-diffusion equation, Euqation (8), for the CSM with α = 0.25 and α = 1, are depicted in Figures 7 and 8    In the case of the mixer with a lower aspect ratio (α = 0.25), we observed an increased rotational motion of the fluids from the velocity vector plots (see Figure 3) and an increased interfacial length in particle tracking, as the Reynolds number increases. Figure 7 illustrates enhanced mixing with increasing inertia, which is caused by the combined effect of the rotational motion and back flows in the L-shaped channel. Comparing mixing at Re = 1 with that of Re = 70, a significant enhancement in mixing is observed. Thus, we can conclude that inertia has a positive influence on the mixing performance of the mixer.
When using a higher aspect ratio (α = 1), mixing at the intermediate Reynolds number (Re = 30) is not as effective as that of the other two Reynolds numbers. When inertia is negligible, as reported in [25], mixing in the CSM with α = 1 relies on the split-and-recombine (SAR) principle [33]. The formation of laminated patterns are disturbed by the rotational flow at Re = 30, leading to less effective mixing as compared with the mixing at Re = 1. However, flow beyond Re = 30 results in better mixing than that at the lower Reynolds number flows. This is due to the effective stirring caused by the rotational motion and back flows in the L-shaped channel, as observed in the mixing analysis of the CSM with α = 0.25 and Re = 70. Now, the remaining task is to quantify the progress of mixing using the concentration distributions. The intensity of segregation [34] is employed to characterize the performance of mixing that is dependent on the aspect ratio (α) of the channel and the Reynolds number (Re). The intensity of segregation I is defined as the second-moment variance of concentration distribution at the outlet surface of a periodic unit: In Euqation (12), σ 2 c is the variance in the concentration over the rectangular cross-sectional area Ω c , defined by: where A c is the area of the cross-section, c(x) the concentration at a position x andc the average concentration defined by:c The denominator of Euqation (12) is the variance of a completely segregated system, where c is either zero or one. In equal volume mixing,c = 0.5, which is the case in this study. The value of I always varies from unity to zero with the progress of mixing. Thus, the intensity of segregation I is one when segregation is complete (unmixed state), while I = 0 when the concentration is uniform (c =c at all points, ideally mixed state). The intensity of segregation (I) is a measure of the deviation of the local concentration from the perfectly mixed situation, which represents a homogeneous state of the mixture [17,30]. The evolution of the intensity of segregation I along the normalized down-channel position z/L is plotted in Figure 9 to quantitatively compare the relative performance of the CSM at five Reynolds numbers, Re = 1, 10, 30, 50 and 70. We also plot the intensity of segregation at a down-channel position z = 5L as a function of the Reynolds number (see Figure 10). In the CSM with α = 0.25, the rate of mixing does not show any significant change up to a Reynolds number of Re = 50 (see Figure 9a). A noticeable increase in the mixing performance is observed in a flow around Re = 70, confirming the qualitative estimation of the mixing performance in the low aspect ratio channel, as mentioned above. While, in the CSM with α = 1, the rate of mixing is the fastest at Re = 70 and the slowest at Re = 30. Figure 10 depicts the intensity of segregation at z = 5L as a function of the Reynolds number at the three values of α. From Figure 10, we conclude that, in the CSM with higher aspect ratio, inertia decelerates the progress of mixing in a lower Reynolds number regime, such as Re ≤ 30, but inertia accelerates the mixing above a threshold value of the Reynolds number, Re = 30. For readers who are interested in the pressure drop with the increase of the aspect ratio and the Reynolds number, we also plot the pressure drop along the down-channel direction in Figure 11. In this figure, the ordinate is the dimensionless pressure p * , defined by p * = p/∆p 1 , where ∆p 1 is the pressure drop in a single periodic unit of each channel geometry at a reference Reynolds number Re = 1. As reported in [35], the performance of a mixer can be evaluated using two different criteria, energy consumption measured by the pressure drop and its compactness measured by the length of the mixer. In microfluidic applications where length becomes a major issue, mixing in a short length at the expense of higher pressure drop is an option if the compactness of a microfluidic device is of great importance. The reduced mixing length, which is the minimum length needed to achieve a desired mixing state (in this study, measured by the intensity of segregation) may compensate for the higher pressure drop in a higher Reynolds number flow.

Conclusions
We numerically investigated the flow and mixing characteristics of the chaotic serpentine mixer (CSM) in non-creep flow conditions. Fluid flow and mixing in the CSM were analyzed for the varying Reynolds number, ranging from 0.1 to 70. As for the geometrical parameters, we focused on the effect of the aspect ratio of the channel since; in the previous study, it was found to be the most critical factor to have an influence on the mixing. The periodic velocity field was obtained for each spatially periodic channel geometry and used in interface tracking and a mixing analysis of a specific mixer.
The Poincaré section at one representative channel geometry that is composed of an L-shaped channel and a bypass channel shows that the flow in the channel induces an overall rotational motion of fluid particles. With the increasing Reynolds number, the strength of the rotation motion increased regardless of the aspect ratio for the channel. Two vortical flows were observed in a channel with α = 0.25, while a single vortical flow was observed in a channel with α = 1, as the Reynolds number increases. The progress of mixing was visualized by plotting the manifold of the deforming interface between two fluids and the cross-sectional distribution of the concentration field obtained as a solution of the advection-diffusion equation. A measure of mixing, called the intensity of segregation, was used to characterize the progress of mixing. As the Reynolds number increases, mixing in the CSM with α = 0.25 shows a noticeable increase above a threshold value of the Reynolds number (Re = 50), as compared with that in the creeping flow regime. This is due to the combined effect of the enhanced rotational motion of fluid particles and back flows near the bends of the channel driven by inertia. In the case of a mixer with α = 1, the intensity of segregation had its maximum around Re = 30, implying that inertia does not always have a positive influence on mixing in this mixer. The CSM with α = 1 always showed better mixing than the mixer with α = 0.25 within the range of the Reynolds number investigated in this study. Finally, it is worthwhile to mention that the CSM with α = 0.25, showing poor mixing in the creeping flow regime, can be an efficient mixer around Re = 70. A microchannel with a lower aspect ratio is favored in mass production using hot embossing or micro-injection molding, since the microchannel can be easily ejected from a mold insert.