Vortex generation in a finitely extensible nonlinear elastic Peterlin fluid initially at rest

It is well known that the mixing of two or more species in flows at low Reynolds numbers cannot be easily achieved since inertial effects are essentially absent and molecular diffusion is slow. To achieve mixing in Newtonian fluids under these circumstances requires innovative new ideas such as the use of external body forces (eg, electromagnetic mixers) or the stretching and folding of fluid elements (eg, chaotic advection). For non‐Newtonian fluids with elasticity, mixing can be achieved by enabling the emergence of elastic instabilities that results in chaotic flows in which mixing is significantly enhanced. In this work, our goal is to demonstrate that clearly identifiable vortical structures (eg, vortex rings) can be generated in a viscoelastic fluid initially at rest by the release of elastic stresses. In turn, these vortex motions promote bulk mixing by transporting fluid elements from one location to another more efficiently than diffusion alone. We demonstrate this first theoretically by using the finitely extensible nonlinear elastic Peterlin (FENE‐P) model to show that elastic forces can generate torque. Using this model, we derive an expression for the time rate of change of vorticity in an elastic fluid initially at rest caused by a sudden release of stored elastic stress. This process can be thought of as the release of elastic energy from a stretched rubber band that is suddenly cut at its center. We confirm this ansatz by performing a series of direct numerical simulations based on an in‐house pseudo‐spectral code that couples the FENE‐P model to the equations of motion for an incompressible fluid. The simulations reveal that a pair of vortex rings traveling in opposite directions, with Reynolds numbers on the order of one, is generated from the sudden release of elastic stresses. Secondary vortical structures are also generated. In the concluding section of this work, we address the potential for vortex motions generated by elastic stresses to promote mixing in microflows, and we describe a possible experiment that may demonstrate this effect.


INTRODUCTION
The primary objective of this work is to show that inherent fluid elasticity can lead in a straightforward manner to the formation of coherent vortex structures in microflows, despite the very low Reynolds numbers associated with these flows. In a nutshell, this can be achieved by recognizing that elastic forces can generate torques on fluid elements leading to well-defined vortex motions. It is important to place the current work in the context of the large body of work on microflow physics. For this reason, we offer in this Introduction a review of the physics of microflows in which the significance of the Reynolds and Peclet numbers are discussed. The need for innovative methods to promote mixing in these flows is emphasized in Section 1.1 of the Introduction. The effects of elasticity are discussed next in Section 1.2, and here we emphasize the emergence of purely elastic instabilities. Recent experiments in which chaotic flows were observed in a dilute polymer solution with rectilinear streamlines are also discussed. The objectives of the current work are expanded upon next in Section 1.3, followed by a review (Section 1.4) of the finitely extensible nonlinear elastic Peterlin (FENE-P) model which describes the flow dynamics of dilute polymer solutions. Finally, an outline of the body of the paper is given in Section 1.5.

Review of microflow physics
Microfluidics, a field of investigation devoted to understand transport phenomena in microdevices, has grown rapidly in recent decades and now constitutes a vast literature composed of thousands of research articles. Much of this effort is driven by a wide variety of important applications such as those in biotechnology, 1,2 biosensor development, 3 and in the design of micro-electro-mechanical-systems. 4,5 Many of these applications are based on fluid flow through microchannels or other similar configurations. As the length scale of these devices is reduced, the physical processes which were applicable at larger scales are reduced in importance, as other processes associated with viscous forces become more dominant. However, as pointed out in a recent review by Squires and Quake, 6 the reduction of inertial effects caused by the system scale reduction, does not render these flows uninteresting. Indeed, at these length scales other physical forces and processes (eg, capillarity, buoyancy, fluid elasticity, electro-chemical processes), previously masked by viscous forces, come into play. In the current work, we have chosen to focus primarily on certain aspects of fluid elasticity. More specifically, our focus is on the issue of mixing at low Reynolds numbers, as well as the possibility of understanding fluid instabilities in viscoelastic flows dominated by viscous effects. As a result, this work will concentrate primarily on the fluid physics of such flows, in the spirit of Squires and Quake, 6 and Stone et al, 7 rather than on the possible uses of the phenomena to construct specific microfluidic devices. We envision that a deeper understanding of the physics in such highly viscous flows is needed before applications based on new phenomena can be devised. In that vein, we are taking a first step along these lines here with respect to certain phenomena associated with fluid elasticity. An understanding of mixing in microfluidic flows involves considering several nondimensional numbers: Reynolds number (Re), Peclet number (Pe), and the ratio of these, the Schmidt number (Sc = Pe/Re).The effects of inertia in microfluidic flows can be understood by considering the Reynolds number, = , where U is a characteristic flow velocity, L is a characteristic dimension, and is the kinematic viscosity of the fluid. In a qualitative sense, the Reynolds number can be thought of as giving the ratio of inertial forces to viscous forces. Generally, inertial forces can be neglected when Re ≲ 1 (Squires and Quake, 6(p981) ). To get some idea of the size of the Reynolds number in typical microfluidic flows, we consider the flow of water in a channel flow of height L = 100 μm, with a flow velocity U = 1 cm/s, for which Re = 1. Although this is indeed a small Reynolds number, it is typically on the high end for microfluidic flows. Typical experiments employ aqueous solutions of glycol, whose viscosity is more than 100 times that of water, and use flow speeds as low as 1 mm/s, giving a Reynolds number on the order of Re = 10 −3 . These estimates show that, for all practical purposes, inertial effects in most microfluidic flows can safely be neglected. The near absence of inertial forces in such flows essentially eliminates the possibility of inertial instabilities, which can lead to turbulence in higher Reynolds number flows.
In microfluidics applications concerned with mixing, it is important to consider the rate at which substances diffuse into one another. For this purpose, the Peclet number = D , where D is the molecular diffusivity (diffusion coefficient) of a substance A into a substance B, needs to be considered along with the Reynolds number. For example, the diffusivity of the fluorescent dye fluorescein in water is significantly smaller than the diffusivity of momentum in water, which can be taken to be the kinematic viscosity of pure water at 25 • C. In this case, the Schmidt number is Sc ≈ 2500, which indicates that, assuming the flow remains laminar (which is always the case in most microfluidic flows), diffusion dominates over inertia in controlling the intermingling or mixing of fluorescein in water. In such cases, the ratio of the width of the diffusive boundary layer, D , to the downstream flow distance, x, can be shown by a simple scaling argument to be This estimate can be used to show that the distance downstream at which complete mixing can be expected, x m , is given by x m = Pe h h, where h = ℎ D and the channel width is h. Therefore, the mixing length scales linearly with the Peclet number. Under such conditions, it may be possible to design a microfluidic device such that two species remain relatively un-mixed for downstream distances which are large compared to the channel width. These considerations lead to the conclusion that if rapid mixing is not desired, devices should be designed such that Re ≲ 1, and Pe ≫ 1. An excellent example of this type of situation is given by the T-sensor. 8 On the other hand, many applications require the mixing between two or more substances in microfluidic systems. For example, if compounds are required to react on time scales much shorter than the diffusion time scale, rapid mixing of the reactants is essential. But at low Reynolds numbers, inertial instabilities and turbulence are no longer available for mixing. To overcome this fundamental problem, both passive and active methods have been devised. From a fundamental point of view, all mixing ultimately occurs on a molecular scale. For example, we can imagine that a droplet of liquid containing a dye (say fluorescein) falls into a container of pure water. If the flow is unstable, the droplet will be stretched and contorted to form a complex random structure consisting of many small tendrils whose thicknesses or diameters are many times smaller than that of the original droplet. The dye eventually diffuses into the surrounding liquid from each tendril through a purely molecular process. The breakup of the drop speeds up the diffusion of dye through a stretching and folding process, which is sometimes referred to as the Baker's transformation. 9 The original droplet is broken into smaller tendrils of length l, thereby increasing the surface area of the original droplet. From the theory of diffusion, we can understand why this process enhances diffusion by noting that the time scale for diffusion can be estimated as t D ∼ l 2 D , where l is the distance over which the substance has diffused and D is the diffusivity. As an example, this indicates that if the tendril length decreases by one-half, the diffusion time decreases by one-quarter, thus significantly decreasing the time for the complete homogenization of the diffusing species. [10][11][12][13][14][15][16][17][18][19][20] These ideas have been used recently by Stroock et al 21 to develop a passive mixing device called the staggered herringbone mixer. In this device a herringbone pattern was lithographically generated on the wall of a microchannel. This type of wall roughness produced a component of crossflow which was alternately reversed by a periodic change in the herringbone pattern, and ultimately gave rise to a chaotic pattern of tendrils, a process known as chaotic advection. 9,11,12,22,23 The rate at which tendrils were generated downstream of the channel entrance was shown to be proportional to the logarithm of the Peclet number, a significant improvement over the simpler case for nonchaotic advection, 8 and similar to the rate at which tendrils are expected to be produced in inertially driven turbulent flows. Other examples of passive mixing include the rotary mixer, 24 and other forms of micromixing. [25][26][27] Greater control over the mixing process can be gained by introducing active forcing. A particularly clear example of this is the use of magneto-hydrodynamics principles in a microcavity. 28 In that work, it was shown that the Lorenz body force can be generated with sufficient magnitude in a microflow to overcome viscous forces so that mixing could be achieved when needed. Other active methods such as thermal, 29,30 and acoustic [31][32][33] have been developed as well.

Effects of fluid elasticity
In the review above, we have been concerned entirely with the behavior of Newtonian fluids in microflows. However, non-Newtonian effects can be significant in these flows. In particular, we focus on fluids with inherent elasticity, otherwise known as viscoelastic fluids. More specifically, we focus on flows of dilute polymer solutions. However, before turning to microflows, it is important to note that these kinds of fluids have significant unexpected behavior at high Reynolds numbers. The early experiments of Toms 34 revealed that introducing only 10 ppm of high molecular weight polymer into a turbulent pipe flow can give rise to significant drag reduction (30%-40%). Since that time, many hundreds of papers have been written on this puzzling phenomenon, with much of the work summarized in several review articles through the years. [35][36][37] It might be expected that if inertial effects are removed, the effects of fluid elasticity would be significantly diminished. In fact, this does not turn out to be the case. Instabilities due solely to the inherent elasticity of the fluid have been discovered in a wide number of flows over roughly the past three decades. [38][39][40][41][42][43][44][45][46][47][48][49][50][51][52][53][54] In these flows, the streamlines are generally curved, and the instabilities are thought to occur as perturbations are amplified as they are coupled to the hoop stress associated with the mean flow. However, recent experiments 55,56 have revealed instabilities in the flow of a viscoelastic fluid within microchannels with rectilinear streamlines. In these experiments, the working fluid was a dilute solution of polyacrylamide in a 90%:10% glycerol-water solution that resulted in a Reynolds number of about 10 −2 , clearly indicating the absence of inertial effects. Large amplitude perturbations to the base flow were introduced by an array of cylinders mounted vertically near the entrance to the channel. It was observed that the flow became unstable for a sufficiently high Weissenberg number (defined in References 55 and 56 as =̇, where is the fluid relaxation time anḋis the mean flow strain rate), Wi ≥ 5.2, and for a sufficiently energetic flow disturbance (15 cylinders). In fact, the instability was observed far downstream (200 channel widths) from the channel entrance, where the flow transitioned to a chaotic, turbulent-like state. It was observed that the elongational strain rates were apparently amplified downstream of the cylinders, indicating the possibility of significant stretching of the polymer molecules. Experimental evidence thus far suggests that the observed chaotic flow is due to a nonlinear subcritical instability. 55,56

Primary objectives of the current work
The experiments which highlight purely elastic instabilities, as well as the difficulties discussed above regarding the need for enhanced mixing in microflows, have motivated us to investigate a somewhat different aspect of low Reynolds number viscoelastic flows. We can summarize these ideas in two related questions as follows: (a) What is the nature of the flows produced by the release of elastic energy in low Reynolds number flows? and (b) Can the release of elastic energy in microflows be converted into coherent vortex motions? If so, what is the nature of these vortical flows? These questions will be answered once the concept of converting stored elastic energy into flow kinetic energy is understood.
It is important to note that here we are not referring to the vorticity associated with the mean flow, as for example in the parabolic Poiseuillian velocity profiles observed in low Reynolds number pipe or channel flows. Instead we are referring to fluid motions with closed streamlines such as flows associated with vortex rings. Although the precise definition of what is meant by a vortex is not available, these motions can be reasonably well identified by using the 2 criterion, 57 or the Q-criterion, 58 which have been successfully used to identify coherent vortex motions embedded in complex flows, including turbulent ones. Such vortex motions have already been used in microflows 21,24,25 to enhance mixing. However, our ansatz asserts that such vortex motions in microflows may arise from the release of elastic energy alone, thus setting our work apart from previous ideas.
Recent works, 59,60 related to the questions posed above, have shown that torques arising from viscoelastic stresses can act to reduce the strength of hairpin vortices in high Reynolds number turbulent flows. They postulate that viscoelastic forces generate a torque opposite to the direction of rotation of drag producing vortices, which reduces the rotational energy of these structures. As far as we are aware, this is the first example in which the torque-producing potential of viscoelastic forces has been explored. In this work we explore the nature of such forces at the low Reynolds numbers associated with microflows, and in flow situations which do not involve other complex fluid motions such as those that arise in unstable or turbulent flows. In addition, we explore the ways in which viscoelastic stress can produce coherent vortex structures as opposed to ways in which these forces can affect already existing coherent vortical structures as in the above-mentioned work. 59,60 Although greater mathematical detail will be given below in the Problem Formulation Section, it is important to recall, in general terms, the characteristics of the viscoelastic force field, f , which allows for the generation of vortex motion. First, the force acting on a fluid due to a stress field given by the stress tensor T is given by the divergence of the stress f = ⋅ T. This implies that a spatially constant stress produces no force, and that only stress fields with spatial gradients can produce net forces on fluid elements. Furthermore, conservative force fields which possess a scalar potential, , such that f = , cannot generate torques on fluid elements since the curl of such a force field is zero. However, for flows of viscous liquids containing polymers, the tensor representing the polymer stress field does not in general give rise to a conservative force field. It follows that polymeric stress gradients can give rise to fluid torques and therefore also to vorticity. It will be shown in forthcoming sections that under certain conditions, polymer stress fields can be constructed to give rise to distinct and clearly identifiable vortex structures such as vortex rings.

The FENE-P model and its applications
To model the presence of polymers in microflows, we will employ the FENE-P model of polymer-solvent dynamics. 61,62 This model is derived using a kinetic theory approach, and results in the field equations for the evolution of the polymeric conformation and related stress field. In this particular model, the solvent is viewed as being populated with large numbers of linear polymer molecules, each of which is thought of as a dumbbell consisting of two massless spheres, connected by a spring. The solvent can stretch the dumbbells via viscous forces, which then react back on the fluid, thus coupling polymeric effects directly to fluid motions. In this approach, the dumbbells are prevented from stretching to infinite lengths by introducing a nonlinear spring model through the Peterlin function. The FENE-P model has proven effective in recent years in exploring complex phenomena such as the effect of polymer additives in fully turbulent flows, although some unrealistic aspects of the model have been discussed in recent work. 63 High resolution direct numerical simulations (DNS) 64-68 of such flows have confirmed the experimental observations of Toms 34 and others. 35,36 These simulations have shown that the drag reduction phenomenon can be replicated by coupling the polymeric forces modeled as Peterlin dumbbells with the equations of fluid motion. In addition, for use in simulating turbulent high Reynolds number polymeric flows, the FENE-P model has recently been incorporated into Reynolds-Averaged Navier-Stokes models, [69][70][71] and k-epsilon turbulence models, [72][73][74] and has been used to study the effects of polymers on heat transfer. 75 The model has also been used with success in exploring the physics of low Reynolds number flows. 76-79

Outline
The remainder of this paper contains the following sections: Problem Formulation, Numerical Methods and Results, and Summary and Discussion. In the Problem Formulation section, we describe the details of the mathematical model, which involves coupling the FENE-P model for the evolution of a nonlinear dumbbell system to the momentum equation, constrained by the condition of incompressibility. It is shown from these equations that when elastic forces are released in a localized spatial region in an impulsive manner, vortex-ring like structures can be generated in a fluid initially at rest. In the Numerical Methods and Results section, the theory is confirmed by performing a series of three-dimensional DNS using a pseudo-spectral code. The results show that in fact, when elastic energy is released in the manner described, two vortex ring-like structures are generated, which travel in opposite directions as predicted by the theory.
In the current work, we show that forces due to polymeric stresses can generate coherent vortex motions in fluids. This mechanism itself is novel and worthy of in-depth exploration independent of its potential applications. The Summary and Discussion section presents a feasible scenario in which this mechanism may play a role in creating elastic instabilities in flows with rectilinear streamlines, as well as how it may be used to enhance mixing in low Reynolds number flows. To strengthen our argument, we also briefly describe a physical experiment which could potentially confirm in a laboratory the theory and simulations described herein.

General formulation
The problem of interest is governed by the momentum equation and the continuity equation for an incompressible fluid given by Equations (1) and (2), respectively: and where D/Dt is the material derivative, V = (u, v, w) = (u 1 , u 2 , u 3 ) is the fluid velocity in the x, y, z , or x 1 , x 2 , x 3 directions, p is the pressure, T is the stress tensor, and is the density. For a dilute polymer solution, the stress is decomposed into a Newtonian component n , and a polymeric component p , as follows: where n = 2 0 S, 0 is the solution viscosity, is the ratio of the solvent viscosity to the solution viscosity, and S is the rate-of-strain tensor. For a FENE-P fluid the polymeric stress 61,62 is given by: where is the polymer relaxation time, C is the conformation tensor defined as the average over all possible molecular configurations of the product of the end-to-end vectors associated with the polymer molecular length R = √ tr(C), I is the unit tensor, and f(R) = (L 2 − 3)/(L 2 − R 2 ) is the Peterlin function, where L is the maximum allowable molecular extension. In the equations above, and in all subsequent ones, C, L, and R are made nondimensional by the rest length, or square of the length as appropriate, of the polymer molecule.
The evolution of conformation tensor C is governed by: where D p is the polymer diffusivity. The first two terms on the right-hand side of Equation (5) represent the stretching and reorientation of polymer molecules by flow gradients, which are balanced by the third term, the spring relaxation term. The balance between stretching and relaxation implies that with no stretching, the polymers will relax back to their equilibrium length, where the molecules undergo only Brownian motion.

Generation of vorticity from elastic stresses
To demonstrate the generation of vorticity by purely elastic stresses, we first take the curl of Equation (1), while taking into account the decomposition given by Equation (3). For an incompressible fluid of constant density, the evolution equation for the vorticity becomes: where 0 = 0 / is the kinematic viscosity of the solution. Proceeding from left to right, the three terms on the right-hand side of Equation (6) represent vortex-stretching and tilting, diffusion of vorticity, and the generation of vorticity by elastic stresses. According to Equation (6), the last term can be interpreted as a source of vorticity caused by the force per unit volume F = ⋅ p , which in turn is a result of gradients in the polymeric stress p . If this force field has a nonzero curl, it will create local spin (vorticity) in the fluid. More simply, this term represents a net torque on fluid elements, creating local spin. To place the current work in context, we note that there are other well-known nonconservative forces which can serve as sources of vorticity such as Coriolis and buoyancy forces, as well as more recently discovered mechanisms. 80 In the following, a right-handed x, y, z coordinate system will be used interchangeably with the labels x 1 , x 2 , and x 3 , and the terms streamwise, vertical or wall-normal, and spanwise. Using this notation, for a fluid which is initially at rest (V = 0), Equation (6) becomes: where = 0 (1 − )/ , and both f and C ji are to be evaluated at time t = 0. We note that the Peterlin function f varies in space and time through tr(C), so that a detailed expansion of the right-hand side of Equation (7) would involve four terms, three of which depend on the spatial gradients of f . Choosing the maximal allowable polymer length L to be such that L 2 ≫ tr(C), where here we evaluate tr(C) at t = 0, then f ∼1 initially. This will in fact be the case in all simulations to be presented in this work. With this initial condition, Equation (7) becomes: To see how this result can be used to generate a spatially concentrated region of vorticity, consider a rectangular domain (eg, a channel) with domain lengths in the x, y, and z directions of L x , L y , and L z . For purposes of simplicity, we additionally consider here only the effects of the diagonal components of the conformation tensor, C 11 , C 22 , and C 33 in the generation of vorticity, even though the off-diagonal components of the conformation tensor can also be used.
Within the larger rectangular domain, we define a smaller rectangular region b When started from a state of rest, the fluid will be in equilibrium if F I G U R E 1 Schematic view of the geometry used in the simulations. F signifies the force generated by the polymeric stress field that in turn generates vorticity as shown in the x 1 − x 3 plane. The situation shown is that which results from a deficit (K < 0) in polymeric stress in a region of scale size h -- no elastic stresses exit, which is the case when C 11 = C 22 = C 33 = 1, C ij = 0 for i ≠ j.Now consider a different initial state in which at t = 0, C 22 = C 33 = 1, C ij = 0 for i ≠ j, and C 11 = C 0 11 outside the smaller region, and C 11 = C 0 11 + K, inside that region. Since C 11 > 0 and C 0 11 > 0, it should be kept in mind that if K < 0, then it must be true that |K| < C 0 11 . This scenario will create unbounded gradients in C 11 at the boundaries of the smaller region which would not be computationally feasible. However, for the purposes of this example, it should be considered as an idealized case whose purpose is solely illustrative.
In the above scenario, Equation (8) indicates that the vertical and spanwise vorticity components, Ω 2 and Ω 3 , are created initially as follows: and The generation of vorticity in this manner can be thought of as having been generated by two oppositely directed impulsive forces acting on the fluid in the x-direction, as shown in Figure 1, in which the case K < 0 is shown. We also note that such a force is similar to the impulsive spatially localized body force fields that have been previously used to generate vortex rings theoretically, 81, (pp. 205-208) and numerically. 82 However, it is worthwhile noting that in Swearingen et al, 82 a vortex ring was generated using a body force, whereas in our case, vorticity is generated by the gradients of internal elastic stresses.
Since the vortex lines associated with these components of vorticity must be connected, we can expect the resultant flow to consist of two vortex ring-like structures, one propagating to the right and the other to the left. When K > 0, the impulsive forces are directed toward each other and colliding vortex ring-like structures can be expected. These two scenarios, which are created either by a deficit or an excess of elastic stress in a small spatial region, can be interpreted as due to the release of stored elastic energy in the fluid. An analogy, though not perfect, would be to view these vortex generation cases as due to the cutting at the midpoint (K < 0), or the release of the ends (K > 0) of a stretched rubber-band.
If the domain has no-slip walls at y = ∓ L y /2, then vortex ring-like structures can also be created in such a way as to propagate toward the walls. This can be realized through the generation of spanwise (z-direction) and streamwise (x-direction) vorticity using the initial conditions, C 11 = C 33 = 1, C ij = 0 for i ≠ j, and C 22 = C 0 22 outside the smaller region, and C 22 = C 0 22 + K inside that region with K < 0. In this case, the analog of Equation (9a) and (9b) is: and This situation may be of interest in cases where it is desired, for example, to increase heat flux in a channel whose walls differ in temperature. We will not discuss this case further in this paper except to note that the generation of vorticity in this manner clearly illustrates the possibility as stated in the Introduction regarding the goal of mixing enhancement in low Reynolds number flows.

NUMERICAL METHODS AND RESULTS
To illustrate the generation of vorticity by the mechanism described above, we have performed simulations for the case in which only C 11 varies initially as described in Equation (9a) and (9b). For this purpose, Equations 1 to 5 were solved using a pseudo-spectral code 83 in a channel geometry in which all field variables are expanded in Fourier modes in the horizontal (x − z) plane and Chebyshev modes in the vertical (y) direction. Spectral methods exhibit exponential convergence 84 and have been used with success in simulating fully developed turbulent flows. 85 The resolution was 64 × 65 × 64 in the x, y, and z directions, and the corresponding domain lengths normalized by the channel half-width were L x = 4 , L y = 2 and L z = . On the top and bottom walls of the channel, no-slip conditions were imposed on the velocity field and no flux conditions, C ij / x 2 = 0, were applied to the conformation.

Effects of forcing strength and Weissenberg number
To determine the effects of the strength of the forcing, K, and the Weissenberg number, on the generation of vorticity, two sets of simulations were performed. In the first set, the strength of the forcing, K * = |K|/L 2 , was varied while keeping the Weissenberg number fixed, Wi = 0 /h 2 = 4. In the second set, K * = 10 −1 was kept fixed and Wi was varied. In these simulations, the following parameters were fixed as follows: L = 100, = 0.5, and 0 /D p = 2.The polymer diffusivity D p (see Equation (5)) is used to insure numerical stability 86 for all simulations presented in this work. The initial conditions were u = v = w = 0 for the velocity. For components of the conformation we set C 22 = C 33 = 100, C ij = 0 for i ≠ j, and C 11 = C 0 11 = 1200 except in a small cubical region of size h = 0.5 in units of the channel half-width as shown in Figure 1. Inside the cube C 11 = C 0 11 + K, with K < 0 and |K| < C 0 11 . The center of the cubical region corresponds to the center of the computational domain. Smoothing was applied to the edges of the cube using hyperbolic tangent functions.
The overall flow strength is indicated by the time-dependent volume-averaged kinetic energy k(t), and enstrophy e(t), which are defined as follows: where V is the domain volume and Ω x , Ω y , and Ω z are the components of the vorticity. These are used to define a time-dependent Reynolds number R(t), and a time dependent nondimensional vorticity magnitude Ω(t), as follows: R = √ 2k h∕ 0 and Ω = √ e h 2 ∕ 0 . In the flows considered here, the kinetic energy reaches a maximum value k m , which defines a peak Reynolds number R * = √ 2 k m h∕ 0 . A nondimensional time is defined by t * = t 0 /h 2 where t is time.
The temporal evolution of the Reynolds number R, and Ω, are shown for different forcing magnitudes K * , in Figures 2 and 3. Both R(t) and Ω(t) show similar behavior which includes a rapid rise to a peak value followed by a viscous decay. Such trends are indicative of not only the production of flow kinetic energy, but also vorticity, as predicted by Equation (9a) and (9b). The rapid viscous decay of flow energy and vorticity is certainly to be expected, since the maximum Reynolds number for the flow for which K * = 0.1 is R max = V max h/ 0 ≃ 3, where V max is the maximum magnitude of the velocity that occurs at t * = 6 × 10 −2 . Such low Reynolds numbers ensure that viscous forces dominate over inertial and elastic forces in the later stages of these flows.
The physics of the formation and evolution of vortex rings is relatively simple and can be understood in the context of Equation (6). The initial conditions for the elastic stress ensure that at time t = 0 a spatially localized transient impulsive force field, F = ⋅ p , acts on a motionless fluid. This force is specified to generate torques in the fluid by ensuring that × F is nonzero. The geometry of the force-field generates two oppositely directed vortex rings as discussed earlier. All elastic and kinetic energy generated in this manner must eventually be dissipated by viscosity, a result which is clearly observed in Figures 2 and 3.
In Figures 4 and 5, the effects of forcing strength K * , and Weissenberg number Wi, on the peak Reynolds number R * ,are shown. If we consider R * as an indication of flow strength, it is apparent that the flow strength increases almost linearly with the forcing strength, and roughly inversely with Wi. These results are consistent with Equation (9a) and (9b) since the spatial gradients of C 11 are proportional to the forcing strength, and is inversely proportional to the polymer relaxation time. In addition, in all of these simulations, including the case corresponding to the lowest forcing magnitudes, vortex rings were formed as indicated by flow visualizations, and as predicted by the theory described in Section 2.

Effects of maximal polymer extension length (L) and the ratio of viscosities ( )
Additional simulations were performed to determine the effects on the flow of the maximal polymer extension length, L, and the ratio of viscosities, . In a first series of simulations, L was increased from 70 to 100, while keeping the following parameters fixed as follows: K * = 10 −1 , Wi = 4, = 0.5, and 0 /D p = 2. The results, shown in Figures 6-8, indicate that the flow kinetic energy and vorticity decrease as L increases, although the changes in flow energy are relatively small. This is entirely consistent with Equation (9a)

Three-dimensional visualizations of the formation of vortex rings generated by the release of elastic energy
Visualizations of the flow for the case in which K * = 10 −1 , Wi = 4, and 10 −1 ≤ t * ≤ 2.5 × 10 −1 are shown in Figure 12. At each time instant, a surface of constant vorticity magnitude, as well as color maps of the velocity magnitude and velocity vectors on a horizontal plane at the center of the flow are shown. At t * = 10 −1 it is evident that two primary vortex rings are formed soon after the release of elastic energy, each traveling in opposite directions along the x − axis. This is in agreement with Equation (9a) and (9b) which predicts that only vertical (y) and spanwise (z) vorticity should be produced in this case, generating oppositely directed vortex rings. In the last flow snapshot at t * = 2.5 × 10 −1 , two secondary vortex rings

F I G U R E 13
Visualization of the three-dimensional velocity vector field associated with vortex rings formed by the release of elastic energy. In addition to the vector field, isosurfaces of the vorticity magnitude are used to visualize the right and left traveling vortex rings. This image corresponds to the case K * = 10 −1 , Wi = 4, and t * = 2 × 10 −1 . This is the same case shown in Figure 12C. Each vector is color coded according to its magnitude as indicated in the color bar in the image. Velocity is made non-dimensional by h/ 0 appear behind the two primaries. An extensional flow is created in the wake of the primary rings, which generate excess polymeric stress at the center of the domain. The secondary rings, which rotate in a direction opposite to the primary, are formed when this excess elastic energy is released. Other flow visualizations not shown here indicate that vortex rings are formed in all cases, even for the case of lowest forcing magnitude (K * = 2 × 10 −2 ). In Figures 13 and 14, three-dimensional velocity and vorticity vector fields are shown for time instant t * = 2 × 10 −1 for the case K * = 10 −1 and Wi = 4. It is evident, that these fields are consistent with the vector fields that are associated with vortex rings. In particular, the vorticity vectors F I G U R E 14 Visualization of the three-dimensional vorticity vector field associated with vortex rings formed by the release of elastic energy. In addition to the vector fields, isosurfaces of the vorticity magnitude are used to visualize the right and left traveling vortex rings. This image corresponds to the case K * = 10 −1 , Wi = 4, and t * = 2 × 10 −1 . This is the same case shown in Figure 12(c). The vorticity magnitude isosurface is 0.78. Each vector is color coded according to its magnitude as indicated in the color bar in each image. Note that red colored vectors emanate from the interior of the vortex rings where the vorticity is expected to be largest, and blue vectors near the isosurface are indicative of the smaller vorticity near the surface of the vortex rings. Vorticity is made non-dimensional by h 2 / 0 near the isosurfaces of vorticity magnitude are oriented in directions in agreement with the theory (see Equation (9a) and (9b)) which predicts that vertical and spanwise components of vorticity will be dominant.

SUMMARY AND DISCUSSION
The present investigation was motivated by the observation that elastic stress gradients in viscoelastic flows of dilute polymer solutions can give rise to coherent vorticity. Moreover, the mathematical structure of the system of equations governing these flows suggests that coherent vortex motions (eg, vortex rings) may be produced by choosing special conditions for the initial stress distribution in the fluid. Coherent vortex motions have the ability to mix regions of the flow raising the possibility of using such motions to promote mixing in low Reynolds number, high Peclet number flows in microdevices. We are additionally motivated by recent experiments 56 in which instabilities and chaos were observed in polymeric microflows with rectilinear streamlines. The theory presented here, which is based on the well-validated FENE-P model, indicates that initial stress distributions can be chosen which will generate vortex ring-like structures. As we have shown in the numerical simulations, when the initial elastic stress distribution is chosen to contain a localized deficit of stress, two vortex rings traveling in opposite directions are generated upon the release of stored elastic energy. The theory suggests that choosing more complex initial stress distributions will result in coherent vortex motions which could enhance heat transfer and mixing. A series of DNS using pseudo-spectral numerical methods were performed in order to test the theoretical predictions. The simulations show that in fact, two vortex rings traveling in opposite directions are generated when the initial stress distribution is associated with a local stress deficit. In addition, the strength of the vortex motions followed the trends predicted by the theory.
Our results raise several open questions not fully addressed in this work. One concern is the extent to which these coherent vortical structures can be used to enhance mixing and transport in low Reynolds number flows. This is related to the issue concerning whether vortex ring-like structures can be made to propagate in any desired direction. For example, the theory suggests (see Equation (10a) and (10b)) that vortex rings traveling from wall-to-wall can be generated. Their interaction with a wall whose temperature differs from the surrounding fluid could be used to modify the heat transfer rate. It is clear from the theory and simulations that the intensity of the vortical structures is controlled fundamentally by the extent to which elastic energy can be stored initially in a solution. This raises the possibility of increasing the Reynolds numbers of the vortical structures so that instabilities might emerge.
Another issue concerns the potential effects of releasing elastic energy in shear flows with rectilinear streamlines. In the experiments of Qin et al, 56 an array of cylinders placed in a microchannel was shown to produce a chaotic flow downstream of the cylinders. We envision that flow around these cylinders creates an elongational flow component which would then stretch incoming polymer molecules. If true, this stretching effect will create spatially complex polymer stress gradients which could lead to fluid torques and to vorticity. The mean velocity gradient associated with the parabolic profile in the microchannel could also intensify and reorient (stretching and tilting of vorticity as indicated in Equation (6)) this vorticity. This scenario would result in a complex three-dimensional vorticity field and possibly chaos, which was in fact observed. This process can be summarized in three steps as follows: (a) Localized stretching of the polymer molecules, (b) generation of vorticity via polymer stress gradients, and (c) generation of a three-dimensional vorticity field by vortex stretching and tilting of the existing vorticity by the mean flow resulting in a chaotic flow. The ideas outlined above, and related ones, are currently under investigation.
From the perspective of the FENE-P model and the parametric study undertaken in our work, we are also actively working on determining parameters that pertain to specific polymers and fluids used in experiments. 56 We are undertaking a molecular all-atom approach to determine how a prototypical linear polymer such as polyacrylamide behaves when dissolved in a fluid such as glycerol. Of particular interest is the determination of the fluid environment, and temperature on the polymer conformation tensor C. 87 The possibility of experimentally verifying the phenomena described in this work is also being pursued. For example, Sousa et al 88 generated an extensional flow of a dilute polymer solution using a convergent-divergent flow field to stretch polymers at a stagnation point. If stored elastic energy in this region of the flow could be released through the cleavage of the extended polymers, then fluid motions similar to the ones observed in our simulations should be generated. Laser energy acting upon specially designed photocleavable molecules 89 are now under consideration for achieving the release of elastic energy within the fluid. We are currently considering the design of such experiments in future work.