Enhanced Microfluidic Mixing via a Tricritical Spiral Vortex Instability

Experimental measurements and numerical simulations are made on fluid flow through cross-slot devices with a range of aspect (depth:width) ratios, 0.4<alpha<3.87. For low Reynolds numbers Re, the flow is symmetric and a sharp boundary exists between fluid streams entering the cross-slot from opposite directions. Above an alpha-dependent critical value Re_c, the flow undergoes a symmetry-breaking bifurcation (though remains steady and laminar) and a spiral vortex structure develops about the central axis of the outflow channel. An order parameter characterizing the instability grows according to a sixth-order Landau potential, and shows a progression from second order to first order transitions as alpha increases. A tricritical point occurs for alpha ~ 0.55. The spiral vortex acts as a mixing region in the flow field and this phenomenon can be used to drive enhanced mixing in microfluidic devices.

The ability of fluids to mix is greatly enhanced by turbulence, which occurs at large values of the Reynolds number Re ≡ U L/ν, where U and L are characteristic velocity and length scales respectively and ν is the kinematic viscosity of the fluid. Small length scales tend to suppress Re, making it difficult to develop turbulent mixing in microfluidic devices. Achieving efficient mixing via diffusion dominated processes presents a major technological challenge to the expanding field of lab-on-a-chip development.
The cross-slot device is emerging as a promising platform for enhancing mixing of fluids in micro-scale geometries [1,2]. The planar elongational flow field generated by the cross-slot geometry ( Fig. 1(a)) has found applications in many research areas including for studies of macromolecular dynamics [3][4][5], extensional rheometry and elastic instabilities of viscoelastic fluids [6][7][8][9], hydrodynamic trapping [10], and imposing controlled deformations to complex biological structures [11][12][13][14]. It has been known since the early 1990's that such intersecting flows are prone to instability beyond a modest critical Reynolds number Re c ∼ O(10-100). In very deep crossslot flow channels the observed instability takes the form of a stack of three-dimensional vortical structures that appear in the central cross-over region [15,16], while in the related 4-roll mill apparatus a flow field incorporating a single helical region has been observed [17]. In this work, with cross-slots of modest aspect ratios typical of microfluidic devices (depth to width ratios between 0.4 and 3.87), we have found that the instability occurs above an aspect ratio-dependent value of Re c and results in a simple straight vortex that extends downstream along the outlet channels, see Fig. 1(b). This flow instability has been shown to promote mechanical scission of polymer Confocal microscopy is performed in z-planes, which are scanned through the full depth of the device and used to reconstruct an image in the x = 0 plane (shaded (green) region). (b) Threedimensional (3D) rendering of a vortex structure observed for the flow of water at Re = 75.8 in a cross-slot with α = 1. The image is generated from z-plane images spaced at δz = 5µm and has been cropped around the central vortex. The volume shown corresponds to the fluorescently-dyed fluid stream. Also see Movie M1 for an animated version of Fig. 1(b) [19].
chains [18], and has more recently been demonstrated to enhance mixing between the two incoming fluid streams in micro-sized devices [1,2]. Numerical simulations have shown that the mixing quality is higher and Re c is lower in the cross-slots than in the more well-known T-shaped mixing device with equivalent channel dimensions [2].
In this Letter, we report the results of detailed experimental and numerical studies of the spiral vortex flow instability in cross-slots with a range of aspect ratios and over a wide range of Re. In contrast to previous studies [2,15,16], we identify appropriate order parameters that characterize the instability as a function of Re in each case. In particular, we present a systematic analysis in terms of bifurcation theory analogous to the Landau theory of phase transitions. The observed phenomena are well described by a Landau-type sixth-order polynomial potential [20][21][22], with parameters that show the transition develops from second order to first order as α increases, passing through a tricritical point for α ≈ 0.55. These data can be fully described by scaling theory, and the universal scaling function and behavior are measured near the tricritical point. Improved understanding and characterization of stability conditions for flows through intersecting geometries is vital for the optimization of many laboratory microfluidic experiments and also practical lab-on-a-chip designs, including for the specific goal of enhancing the mixing of fluids in channels with small dimensions operating at low Re.
The cross-slot device ( Fig. 1(a)) consists of two bisecting rectangular channels of width w and depth d, with the aspect ratio defined by α = d/w. In the experiments four devices are utilized, all with d = 1.2 mm but with w varied in order to provide α = 0.49, 1.00, 1.85, and 3.87. An inlet length of at least 12.5w ensures a fully-developed flow before the fluid reaches the central region of each device. Newtonian fluid (water) is pumped into two opposing channels (along the y-direction) and exits through the two opposing outlet channels (along the x direction). One of the incoming fluid streams is fluorescently-dyed with rhodamine B (concentration = 10 µM, Sigma) and a laser-scanning confocal microscope (Zeiss LSM 780) is employed to examine the interface between fluid streams where they meet in the central cross-over region. Imaging in closely-spaced planes through the depth of the flow cell (−d/2 ≤ z ≤ d/2) allows accurate reconstruction of an image in the x = 0 plane, see (green) shaded area in Fig. 1(a). Experiments are performed at 24 • C over a range of Re = ρwU/µ, with the fluid density ρ = 997.1 kg m − and the dynamic viscosity µ = 9.1 × 10 −4 Pa s. The average flow velocity U within each channel of the cross-slot is controlled using a precision dual syringe pump (Harvard PHD Ultra). While for low Re the interface between fluid streams is sharp and vertical over the y = 0 plane, beyond a fairly moderate critical value Re c the flow bifurcates and breaks symmetry (but remains steady and laminar) and intricate spiral vortex structures develop, see Fig. 1
The numerical method solves the equations of motion ρDu/Dt = −∇p+µ∇ 2 u and mass conservation ∇· u = 0, for laminar flow of a Newtonian incompressible fluid by using a fully-implicit, second-order finite volume method [7,23,24]. For the inertial flows considered here, relevant features of the method are the treatment of the partial time derivative ∂u/∂t using the three time-level pressure correction algorithm [25] and the treatment of the convective term u· ∇u, which is discretized using the 3rd order high-resolution CUBISTA scheme [26]. In this way we maintain 2nd-order accuracy in space and time.
Non-uniform orthogonal meshes are deployed on the 3D cross-slot geometry of Fig. 1(a), with the central cube (for α = 1) having 25 3 uniform control volumes (CV's) on the base mesh, and 35 3 or 71 3 CV's on more refined meshes. The CV size expands slowly towards the two inlets, along the channel aligned with the y-axis, and towards the two outlets, along the x-axis (see Fig. S1 [19]). Theoretical [27] 2D fully-developed velocity (v th (x, z)) profiles are applied at the inlets, while at the outlets zero axial gradients are assumed for velocity components (∂/∂x = 0), pressure is linearly extrapolated and the total flow rate is forced to satisfy overall mass conservation (Q out = 2Q in , where Q in = U wd for each inlet section). Note that we do not force the flow rate entering via each inlet to divide equally between each outlet, but only that the flow rate through each outlet is equal, as in the experiments. The numerical simulations explore devices with aspect ratios set equal to the four experimental devices as well as additional values of aspect ratio near to the tricritical value. Great care is exercised to guarantee that the numerical results are well converged and do not significantly with the computational mesh. In Fig. S1 [19], we illustrate the various meshes employed and compare the results obtained for a cross-slot with α = 1, showing good accuracy.
In Fig. 2(a-d) we present reconstructed confocal images of the x = 0 plane (spanning −w/2 ≤ y ≤ w/2, −d/2 ≤ z ≤ d/2) for the cross-slot device with α = 1, which show the evolution of the vortical structure as Re increases. Qualitatively similar behavior is observed in all four geometries (see Fig. S2 and Movies M2-M5 [19]): at low Re the interface between fluid streams is sharp and vertical as expected, however as Re increases above a critical value the spiral vortex forms abruptly about the x-axis. With further increases in Re, the spatial extent of the spiral expands, as does the number of turns of the arms. The value of the critical Reynolds number Re c,exp for the appearance of the spiral in the experiments depends on the aspect ratio, but is broadly consistent with values reported previously in related experiments [2,[15][16][17] (Re c,exp ≈ 100, 40, 24 and 26 for α = 0.49, 1.00, 1.85 and 3.87, respectively). In each of the four devices the central vortex forms with a favored orientation about the x-axis (anticlockwise for α = 0.49 and α = 1.00 and clockwise for α = 1.85 and α = 3.87, see Fig. S2). The orientation is presumably biased by some minor geometrical imperfections in the devices and causes the bifurcation to follow a favored branch. In only a few instances in the cross-slots with α = 1.85 and α = 3.87 at higher Re did we observe the vortex to form in the unfavored sense (see Fig. S3 [19] for examples). The vortex structures observed in the x = 0 plane can be well-fitted by elliptical logarithmic spirals (see Fig. S4 [19]).
Numerical simulations result in remarkably good agreement with the experiment, as shown by the streamline plots in Fig. 2(e-h) that can be directly compared with the experimental results in Fig. 2(a-d). Note that in the numerical simulations at α = 1.00 there is hysteresis in the transition and both symmetric and asymmetric steady solutions can be obtained for Re = 42.8 ( Fig. 2(f)), depending on whether Re is quasistatically increased or decreased from below or above the transition, respectively. Both possible solutions are presented in Fig. S5 [19]. Additional numerical plots of velocity vector fields at various Re spanning the transition for the case of α = 1.00 are provided in Fig. S6 [19]. Our experimental method is not amenable to a genuinely quasistatic variation of Re since the flow is necessarily interrupted between image acquisition at each Re in order to refill the syringes with fluid. Therefore our experiment can not resolve the hysteresis in the transition.
We treat our experimental data following the approach suggested by Stroock et al. [28] in order to assess the degree of mixing between the incoming fluid streams, M exp . This parameter is computed using the standard deviation of the pixel intensity evaluated over images such as those shown in Fig. 2(a-d) and Fig. S2 [19]: where I is the grayscale pixel value (normalized between 0 and 1, corresponding to the minimum and maximum grayscale intensities in the images) and indicates an average over all pixels in the field of view. For completely segregated fluid streams, there is a binary distribution of pixel intensities, the standard deviation = 0.5 and hence M exp = 0; this is essentially the case illustrated in Fig. 2(a). Complete mixing between the fluid streams would mean a uniform pixel intensity of I = 0.5 over the entire field of view, therefore a standard deviation of zero and hence M exp = 1. In Fig. 2(i), we plot the mixing parameter, M exp , as a function of the Reynolds number for all four experimental aspect ratios. Fig. 2(i) shows a general reduction in Re c,exp and an increasingly large and rapid increase in M exp as α increases. We also note that the transition is imperfect; there is a small but distinct curvature in M exp near the onset, which is consistent with the presence of their being a favored branch for the transition (clockwise or anticlockwise spirals). We treat the favored branch as the positive branch regardless of the spiral orientation.
In the case of the numerical simulations, we find the growth of the instability is best described by the following parameter: where v max | x=y=0 is the maximum value of the ycomponent of the velocity measured along the z-axis at each Re (for a symmetric flow M num ≡ 0). A pictorial description of M num is given in Fig. S6 [19]. Fig. 2(j) shows a plot of M num as a function of the control parameter ε = (Re−Re c )/Re c for various aspect ratio cross-slot devices. The data in Fig. 2(j) is fitted using a Landautype model with "free energy" F , given by: where the order parameter ψ ≡ M num , in the numerical case. Even terms are included in Eq. (3) since for a perfect system F should be independent of the sign of ψ (i.e. the handedness of the spiral). The lowest order asymmetric "field" term (hψ) can account for imperfections that bias the handedness. In the numerical simulations the vortex can form in the clockwise or anticlockwise orientation with similar probablity, hence h = 0 (for all aspect ratios). For given values of the parameters h, g, and k, the value of ψ(ε) corresponds to the extrema of F , ∂F/∂ψ = 0, giving: The ratio of the coefficients g/k > 0 for α < ∼ 0.55 (corresponding to a forward bifurcation or a second-order transition) and decreases monotonically with increasing α, turning negative for α > ∼ 0.55 at which point the bifurcation becomes backwards, or first order and an increasing large hysteresis loop grows as α increases. The numerical data obtained for α = 0.55 are well-fitted by Eq. (4) with g = 0, (cyan) stars in Fig. 2(j), and thus corresponds to a tricritical transition. Fitting of the numerical data with Eq. (4) provides the values of Re c shown in Fig. S7 [19]. For α > ∼ 0.55, due to the hysteresis loop we can find two values of the critical Reynolds number: Re c denotes the value found for quasistatic increases in Re from below, while Re * c denotes the value found for quasistatic decreases in Re from above the onset. The numerically predicted hysteresis loops for α > ∼ 0.55 can not be observed experimentally, however we find that our experimental critical Reynolds numbers (Re c,exp ) agree very well (within approx. ±5%) with the numerical values of Re * c (Fig. S7). In Fig. 3(a) we fit our experimental mixing parameter with Eq. (4), setting the order parameter ψ ≡ M exp and fixing Re c to the value determined numerically for each aspect ratio device. We introduce a small coefficient 0.002 < h < 0.003 into Eq. (4) that accounts for geometrical imperfections in the experimental devices. In general we find that our experimental order parameter is well described by the same sextic Landau potential that describes the numerical data.
The theory for tricritical points [29] predicts a universal scaling form for the canonical order parameter φ ≡ ψ k/|g| in terms of the control parameter r ≡ εk/g 2 , showing data collapse in the form: In Fig. 3(b), we have plotted both the experimental and numerical results in scaled form and show the comparison with the data collapse prediction. The agreement is excellent and confirms our identification of the mixing transition as a tricritical one.
In summary, we have demonstrated that the spiral instability observed for Newtonian flow in the cross-slot device driven far from equilibrium above a critical Reynolds number, can be well-described by a Landau model analogous to that used near equilibrium tricritical points. We note a certain similarity of our results with those presented for flow in low-aspect ratio Taylor-Couette devices above a critical angular velocity (e.g. [21,22,30]), however the instability in the cross-slot geometry can be harnessed to promote mixing of Newtonian fluids at the quite modest Reynolds numbers accessible in microfluidic devices. Such improved mixing could benefit the efficiency of lab-on-a-chip applications such as chemical synthesis, antibody antigen binding reactions and bioassay. It is likely that bifurcation phenomena reported in related flow geometries (such as in the T-shaped micromixer, for example [2,31]) could also be characterized by the tricritical point formalism presented here.