Abstract
This paper presented a direct numerical simulation (DNS) study on the elasticity-induced irregular flow, passive mixing, and scalar evolution in the curvilinear microchannel. The mixing enhancement was achieved at vanishingly low-Reynolds-number chaotic flow raised by elastic instabilities. Along with the mixing process, the passive scalar transportation carried by the flow was greatly affected by the flow structure and the underlying interaction between microstructures of viscoelastic fluid and flow structure itself. The simulations are conducted for a wide range of viscoelasticity. As the elastic effect exceeds the critical value, the flow tends to a chaotic state, while the evolution of scalar gets strong and fast, showing excellent agreement with experimental results. For the temporal changing of scalar gradients, they vary rapidly in the form of isosurfaces, with the shape of “rolls” in the bulk and evolving into “threads” near the wall. That indicates that the flow fields should be related to the deformation of viscoelastic micromolecules. The probability distribution function analysis between micromolecular deformation and flow field deformation shows that the main direction of molecular stretching is perpendicular to the main direction of flow field deformation. It implies they are weakly correlated, due to the confinement of channel wall.
1. Introduction
Chaotic flow occurs for the viscoelastic fluids at very low Reynolds number (Re). Since the inertial effect is trivial, such kind of flow state is called purely elastic instability or elastic turbulence which have been identified in a vast range of macroscopic and microscopic flow geometries with curvilinear streamlines, such as Taylor-Couette device, parallel-disk, cone-and-disk, and serpentine microchannels [1–3]. The viscoelastic fluids are typically the solutions of polymers or surfactants, and they have analogous flexible microstructures which deform with the flow. Considering the nature of elastic effect of viscoelastic fluid in different flows, these instabilities are driven by large positive primary normal stress difference N1 and suppressed by negative second normal stress difference N2. On the balance of fluid nature and the flow geometries, the elegant Pakdel-McKinley condition unifies the critical threshold for the onset of elastic instabilities [4]. Combination between fluid elasticity and curvatures of device is concluded to be the premise to excite the purely elastic instabilities or elastic turbulence. Without the confinement of curved walls, instabilities can hardly happen for very-low-Reynolds-number flow if no external perturbation is imposed. It is worth noting that the slow three-dimensional recirculation flow of shear-banded fluids in the straight microchannel is due to the interfacially driven undulation caused by the jump in second normal stress across the interface [5, 6]. And if the disturbance is not sufficiently large, the disturbance itself will damp by viscous forces [7]. As the elastic instability grows, the purely elasticity-driven turbulence hereby emerges. Also, recently there are numerical and experimental evidence of that the nonlinear elastic instability and even elastic turbulence could be generated in the straight microchannels [7, 8]. Note that their irregular flow state should be sustained by the outside perturbation, for instance, like numerically introducing a sinusoidal force term into the momentum equation and also constitutive equation or placing a variable number of obstacles in the center microchannel. In contrast, those perturbations are not necessary for curved geometries since they could maintain the nonlinearity by themselves. For the excitation and persistence of elastic turbulence in curved geometries, their functioning mechanism is suggested to be similar to “dynamo effect” [9]. The perturbation imposed in the flow induces the additional wavy deformation on the curved streamlines which associate the flexible long-chain molecular structures. And in turn, such deformations in the long-chain molecules raise local variation of stresses to the flow field as feedback. Thus, perturbation reveals itself in a positive way with the viscoelastic fluid flow, and the flow accelerates the perturbation via the flexible molecular motion. As such nonlinearity of elastic effect increases, the flow is proven to be chaotic, which resembles the developed turbulence condition in the statistic characteristics [3, 4]. However, it behaves stochastic in time, but smooth in space, which differs from the hydrodynamic turbulence.
As aforementioned, the elastic instabilities and turbulence stem from the peculiar fluid structures and their motion, which lead to the fast energy transportation. Therefore, as to the passive particles carried by those elasticity-induced irregular flows, they could be transported much forcefully than that in Newtonian fluids cases. Newtonian fluids are known to completely remain laminar flow at very-low Reynolds. So the elastic turbulence is promising for the mixing process and reagents reaction in the microfluidic applications, besides of the electrokinetic effect, electrophoresis, acoustic stream, magnetic driving and other active methods [10–12] to promote mixing. It could be widely employed for most of microfluidic manipulation circumstances, as a passive way of mixing. In this regard, it is probably the optimal option for microfluidic field where needs to mix fluids by the flow itself.
Chaotic mixing by viscoelastic fluid flow in the microfluidic devices has been intensively investigated, with a wide variety of layouts of microchannels [4, 11, 13]. And by the curved microchannel, efficient mixing enhancement has been achieved due to the chaotic flow. Whilst the mixing effect by elastic turbulence has been studied a lot in aspects of mixing length, mixing time, and scalar statistic moment, studies of scalar evolution and motion with time during the mixing process have been much scarcer, especially in the confined microchannels. Ottino and Wiggins proposed that the key factor to mixing enhancement lies in the effectiveness of the magnitude of cross-sectional flow, that is, to create as many contact interfaces as possible between two fluids in the cross-section of the flow direction [14]. The stretching and folding of the interfaces indicate that the flow structure forms Swiss-roll-like shape, which could be accomplished with the flexibility and elasticity of viscoelastic fluids in the motion. Along with the rolling process, the flow-carried scalars are thus transferred by this kind of “forced” chaotic convection, and their evolution is in a ramp-cliff variation temporally and spatially [15, 16]. As for a developed inertial turbulence, the scalar evolution has been studied in a three-dimensional straight channel flow, and it is found that the scalar dissipation shows a thin-layer structure, which is strongly related to the orientation of scalar gradient [17]. Considering the disparities of how the turbulence flow is driven, the scalar evolution should be accordingly behaving unusual in an elasticity-induced irregular flow. And that should be bound up with the motion and deformation of internal microstructures of working fluids.
For scalar transportation in the viscoelastic fluids flow in the curvilinear microchannel, so far there is no information of scalar gradient evolution and its relationship to fluid nature have been reported, to our knowledge. The present authors are of the opinion that a detailed perspective of scalar transportation for a three-dimensional chaotic viscoelastic flow has not been reached yet. The previous work by coauthors [18] has already made a comprehensive analysis on the vortex structures, scalar transportation by elastic instability in the curved channel. The present work further attempts to get a deep understanding on the scalar transportation evolution and its mechanism. The objective of this paper is to investigate the passive particles (scalar) transportation carried by the irregular viscoelastic fluids flow, combining the numerical simulations with the former experimental results. By direct numerical simulation (DNS) conducted in a curved microchannel unit geometry, the relationship between flow structures and the deformation of flexible molecules is elucidated. Also the mixing effect and mixing efficiency with elasticity variation are displayed, and the comparison with the experiments makes the mixing enhancement straightforward. Furthermore, simulations on various flow conditions would clarify our understanding of the essence of the mixing promotion in the confined curvilinear microchannel and the role of flexible microstructures in the elastically turbulent flows.
2. Numerical Simulation Procedures
Three-dimensional DNS was carried out for viscoelastic fluid flow in the curvilinear microchannel. For a liquid flow, to handle the microscaled flow as a continuum medium is to simulate the flow fluid from the macroscopic view; thus, a nondimensional strategy for the governing equations is considered during the conduction of DNS.
2.1. The Governing Equations
The viscoelastic fluid is usually considered as a solution composing of the flexible long-chain molecular part and simple Newtonian solvent part, for its mathematic description. To simulate an incompressible viscoelastic fluid flow, the governing equations read as follows.
The continuum equation is as follows:
The modified Navier-Stokes equation after introducing of elastic stress term is as follows:
where p is the pressure field and μ is viscosity of the Newtonian solvent. The elastic stress term τ p is to be presented in the constitutive equation of viscoelastic fluid, which depicts the relationship between the micromolecules’ deformation and the elastic stress. And this this relationship reads
where
where α≠0 is the mobility factor and a value of 0.001 is used in the numerical computation.
For the passive scalars θ being transported along with the flow, they are depicted by
where D is the diffusion coefficient of fluids.
2.2. Computing Domain and Nondimensional Procedures
The model for numerical simulation is one unit of curvilinear microchannel, as displayed in Figure 1, where periodic condition is used for both streamwise and spanwise directions, and radial direction is confined by the upper and lower no-slip walls. In the Cartesian coordinates, the dimensions of the geometry is nondimensionally set, with the range of [− 4, 4] in x direction, [− 4, 4] in y direction, and [0, 4] in z direction, respectively, and the computing domain is meshed by 64 × 64 × 64 grids. The computation is hereby carried out in a manner of that in Cartesian coordinates.

The computational model of curved channel flow.
The height of the channel is defined between the walls, as 2H. Hereby, parameters of the half-height H, mean friction velocity
where τ
w
is the mean wall shear stress on the entire wall and the nondimensionalized pressure p* is decomposed into average pressure part
For the nondimensional parameters appearing in as follows, such as Weissenberg number Wiτ, Reynolds number Reτ, Schmidt number Sc, Péclet number Pe, and ratio of solvent viscosity to solution viscosity β, are then defined as
For convenience, Weissenberg number based on the nondimensional mean velocity U m is also defined, read as Wi m = U m Wiτ. According to the previous work [18], Wiτ was regarded as the nominal Wi and Wi m was the effective Wi that could be comparable with experiments.
Therefore, dropping the symbol “*”, the nondimensional governing equations for viscoelastic fluid flow and scalar transportation read as
For the computation of scalar transportation, (12) is taken into calculation after the flow reached the numerically steady status. Specifically, the scalar θ is initially set to be 0 and 1 for the upper and lower half channel, respectively (also shown in Figure 5 later). And as time marching, it evolves with the deformation of carrier flow.
2.3. Numerical Procedures
An at least second-order difference method is adopted to discretize the differential equations. For spatial variables, the convective term and the conformation tensor transportation equation, the QUICK scheme, and MINMOD scheme are utilized, respectively. And for time-related terms, the second-order Crank-Nicholson method is used, and the time step is set as 10−5. Also an artificial diffusion term κ∇2C (κ is the relaxation factor) is added in (11) to eliminate the solution unsteadiness, where a sufficiently small artificial coefficient (less than 10−3) does not qualitatively change the flow behaviors [19].
For the specific flow conditions, the details are listed in Table 1.
Numerical cases and relative parameters.
3. Results and Discussions
In an early experimental study on the irregular flow field in the curvilinear microchannel, Li et al. [20] have shown the velocity distributions at the different planes along depth direction by micro-PIV measurement, indicating that a disordered flow state develops for the viscoelastic fluids even at an ultralow flow rate. The previous DNS study [18] on the elastic instability inception and development has clearly explained the chaotic state started at a critical Weissenberg number of about 2, after that the friction factor C f = 2/U m 2 grew sharply due to the generation of unstable secondary flow in the cross-section, where the secondary flow was perpendicular in the streamwise direction [18]. For the carried scalars, the temporal evolution of the ensemble root mean square (rms) of concentration decays exponentially, being excellently consistent with the experimental results [3]. Consequently, while the contribution of viscous effect on the kinetic energy (also mixing effect) maintained at low level, the elastic effect contribution predominated the flow and the inertial effect kept approximately zero. That could suggest that the flow almost depended on the motion of the flexible long-chain molecules of the solute, such as the stretching by the flow and curling or oscillating by the geometry. This would be the premise of the current work, since the scalar motion is completely dependent on the flow structures.
3.1. Flow Structures
To illustrate how the scalar behaves in the elastically chaotic flow, the flow structures should pay special attention at the outset. The substantial stretching of the flexible long-chain molecules will be attained as the flow state turns more irregular. This can be depicted by the trace of the molecular conformation tensor tr(

Extension of flexible molecules versus Wi; inset is the relative importance of elastic stress to viscous stress. The right is the extension contour for the flow field, tr(

Illustration of stretched stripes near the inner wall.
Along with the strong extension of viscoelastic fluid, the fluctuation of scalar concentration in the spatial field would be affected at various length scales. This can be illuminated by the decay of the two-dimensional scalar spectral function Eθ(k). According to Williams et al. [22], it is defined as
where
Therein,

Scalar power spectra by experiments in the wave number domain at different flow units (a), and the fling stripes in the burst events (b). Flow condition is about Pe = 4.78 × 107.

Scalar and scalar gradient for initial condition (t = 0).
Figure 4(a) displays the experimental results of scalar power spectra, where N represents the different flow units along the flow direction. As the flow gets fully developed, they show the similar decay in the space. One could identify the scalar spectra function decay with an exponential index of −0.95, which fairly agrees with the theoretical value of −1 predicted by Batchelor and Kraichnan [23, 24]. However, it should be noticed that the scalar spectra decay shows wiggles at the small wave number region, implying the stretched stripes with concentrated scalars fling in the bulk flow, according to relatively large-scale motion. This can also be explained by the fling motion as shown in Figure 4(b). In the experiments, the time-dependent motion of stripes is from strongly extended long-chain molecules where the fluorescent particles (scalars) concentrate. They are captured by using a high-speed CMOS camera. It can be seen that these stripes are mainly extended near the region of small curvature radius. The curved stripes are showing the fling motion from near-wall region to the bulk, and they have a quasiperiodic behavior. The above finding is consistent with the simulation results presented above.
3.2. Scalar Gradient Structures
To calculate the gradient for (5), one can get
where
In Figure 6, the contours represent the temporal evolutions of scalar concentration under some certain flow conditions, where the parameter t indicates the nondimensional computation time. As for a Newtonian fluid (Wiτ = 0), the carried scalars are in the motion of diffusion. The interface between the two layers of differently concentrated scalars becomes smearing, while the newly generated layers remain parallel along the flow direction. There is no spanwise motion found for all the time instant. In contrast, for a viscoelastic fluid flow, its state is dominated by the convection effect, and the scalars are promptly transported by the vortices to reach uniform distribution in the whole geometry. The flow components both in spanwise direction and radial direction are much stronger than Newtonian case, and scalar motion is quite fast.

Scalar contour versus time at different flow conditions.
Comparing with the previous work on the flow structures [18], it is obvious that the evolution of secondary flow structures show the similar features with that of scalar transportation. For the flow inside the curved channel, the new vortex structures with opposite rotating directions continuously generate near the wall and then move to the opposite walls. It would imply that the scalar is effectively transported with such irregular flow. Thus the scalar gradients in the flow geometry should relate to these special flow structures.
Accordingly, the scalar gradient evolution is visualized in the forms of isosurfaces of |ϕ| (Figure 7), both with full range of gradient value and a cut-off threshold of two times of average gradient value TH = 2〈|ϕ|〉 (TH represents the threshold). As expected, stratified patterns of scalar gradient form for the Newtonian fluid flow. In consequence of long-time decay, they transform into dispersively fine granular structures. This is apparently dissimilar to that in the developed turbulence flow and indicates that the decay of scalar gradient driven by diffusion seems much local in its motion. For viscoelastic fluid, the scalar gradient structures shape into closed curved surfaces or tubes, which subsequently change into flat tubes and then long stripes. These structures are different from that found in the developed inertial turbulence [17]. In inertial channel turbulence, the scalar gradients evolve into thin layer structures. Thus the motion and deformation of viscoelastic fluid are crucially responsible for the scalar gradient structures, and this will be addressed in the next section.

Scalar gradient structures versus time at different flow conditions. The isosurfaces are set to all range of scalar gradient magnitude and 2 times mean value of their magnitude as threshold (TH = 2〈|ϕ|〉), respectively.
As a whole, the scalar gradient evolves as follows in a confined elasticity-induced chaotic flow. The structures of the scalar gradient initiate from the closed surfaces in the center which are characterized by the vortices’ development (rolling) and then are extended in streamwise direction and compressed in radial direction by the carried flow. When the isosurfaces with large scalar gradient vanish in the bulk, some remain near the wall boundaries. That means that the scalars are transported universally by flow extension and rotation movements in the microchannel and finally concentrated at places where the shear effect dominates. Moreover, the thin threads of scalar gradient are probably a result of microstructure stretching effect, relating with the flow structures.
3.3. Characteristics of Microstructures
From the discussion above, the flow structures and the deformation of flexible molecules should be to a certain extent connected to each other. It is worth investigating such a relationship by using the probability distribution function (PDF) analysis between micromolecular conformation tensor
Figure 8 displays the PDFs of cosine value of the angles between eigenvectors

PDFs of cosine of the angles between eigenvectors
To further elucidate this viewpoint, the relationship between the vorticity

PDFs of cosine of the angles between eigenvectors
As an interesting side, we would like to investigate how the scalar gradients correlate with velocity vectors, as shown in Figure 10. It suggests that the direction of velocity is always orthogonal to the scalar gradient. Since the scalars are carried by the fluid motion, the scalars transportation speed changes at different velocity layers. Therefore, in this shear-dominated flow, the scalar gradients turn out the layer structures. As passively transported by the vortices of secondary flow, they form a closed tube appearance in the center of channel.

PDFs of cosine of the angles between velocity vectors
4. Conclusion
In summary, this work focuses on the passive particle motion in an irregular flow with viscoelastic fluid. Purely elasticity-induced chaotic flow and scalar mixing in the curvilinear microchannel has been investigated by a three-dimensional DNS study. The scalar evolution characteristics are well consistent with the experimental findings. For the structures of scalar gradients, their isosurfaces vary from “rolls” in the bulk to “threads” near the wall. The weak correlation between micromolecular deformation and flow field deformation tend to imply the importance of geometry layout. And the strong secondary flow again indicates stretching of viscoelastic micromolecules works as a medium for the enhancement of scalar transportation in the flow field. With the analysis on the PDFs of eigenvectors’ parallelism, the mechanism for scalar evolution, molecules’ stretching, and the scalar gradients’ structures are elucidated.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Footnotes
Acknowledgments
The authors gratefully acknowledge the financial support for this work by the National Natural Science Foundation of China (51076036), Foundation for Innovative Research Groups of the National Natural Science Foundation of China (51121004), China Postdoctoral Science Foundation (2013M541374), and Heilongjiang Postdoctoral Fund (LBH-Z12110) of China.
