Numerical simulation of vortex rings in non-Newtonian fluids

The impulsive viscous flow through an orifice produces vortex rings that self-propagate until total dissipation of the vorticity. This work aims to study numerically such vortex rings for different types of non-Newtonian fluids, including the power-law, Carreau and simplified Phan-Thien-Tanner (PTT) rheological models, with particular focus on the post-formation phase. The simulations were carried out with the finite-volume method, considering small stroke ratios ( L / D ≤ 4) and laminar flow conditions ( Re G ≤ 10 3 ), while parametrically varying shear-thinning, inertia and elasticity. The vortex rings generated in power-law fluids revealed some peculiar features compared to Newtonian fluids, such as a faster decay of the total circulation, a reduction of the axial reach and a faster radial expansion. The behavior in Carreau fluids was found to be bounded between that of power-law and Newtonian fluids, with the dimensionless Carreau number controlling the distance to each of these two limits. The vortex rings in PTT fluids showed the most disruptive behavior compared to Newtonian fluids, which resulted from a combined effect between inertia, elasticity and viscous dissipation. Depending on the Reynolds and Deborah numbers, the dye patterns of the vortex rings can either move continuously forward or unwind and invert their trajectory at some point. Elasticity resists the self-induced motion of the vortex rings, lowering the axial reach and creating disperse patterns of vorticity. Overall, this work shows that the particular non-Newtonian rheology of a fluid can modify the vortex ring behavior typically observed in Newtonian fluids, confirming qualitatively some experimental observations.


Introduction
Vortex rings are coherent toroidal structures of fluid spinning around a closed loop.These vortices can be observed in nature, e.g.those generated by dolphins and volcanos, and can be also generated artificially, as in virtually any impulsive jet flow.A simple, yet popular setup to generate vortex rings in the lab consists of a piston-cylinder arrangement (Fig. 1a), where the piston motion inside the cylinder drives fluid into a reservoir during a finite amount of time [1,2].In this setup, the ratio between the piston stroke and the cylinder diameter at the exit is known as the stroke ratio (L/D).
Vortex rings formed in Newtonian fluids, henceforth called Newtonian vortex rings, have been the subject of numerous experimental (e.g.[1,[3][4][5][6][7]), theoretical (e.g.[8][9][10][11][12][13][14][15][16][17][18][19]) and numerical studies (e.g.[20][21][22][23][24][25]), which analyzed the effect of several parameters on the different stages of their life, from generation to dissipation.Much of the underlying physics of these structures is now understood and interesting features were unveiled.Perhaps one of the most intriguing observations, shown experimentally [1] and later confirmed numerically [22], was the existence of a nearly constant formation number, i.e. a stroke ratio beyond which the circulation transported by vortex rings saturates.The numerical value of this nearly constant formation number was confirmed to be ~4 in several works [1,22], with little variation for a number of different conditions [22].Any additional circulation injected in the flow by a stroke ratio higher than the formation number is accumulated in a trailing wake of vorticity.Thus, the formation number also delimits the stroke ratio above which it is not possible to obtain a single vortex ring without a tail (at high Reynolds number).The reviews by Velasco Fuentes [26], Meleshko et al. [27], Lim and Nickels [28] and Shariff and Leonard [29] summarize most of the principal achievements made during the past years in the investigation of Newtonian vortex rings.
Contrarily to Newtonian vortex rings, the formation of vortex rings in non-Newtonian fluids (henceforth called non-Newtonian vortex rings) has received significantly less attention and much is still left to explore.Indeed, there is only a handful of studies devoted to non-Newtonian fluids and most of them emerged during the last decade.This is somewhat surprising, as in practice vortex rings are generated in non-Newtonian fluids in many cases, as for example the blood vortex rings produced during the left ventricle diastole [30,31].In addition, complex fluids often exhibit exotic responses, quite different from the corresponding Newtonian behavior (e.g.[32][33][34]).Therefore, there is both practical and scientific interest in studying non-Newtonian vortex rings.In the following paragraphs, we will briefly review the literature devoted to non-Newtonian vortex ringsall studies are experimental, although some include a theoretical analysis.For the sake of simplicity, we will use Re to refer generically to the Reynolds number reported in each work.Usually, Re uses as length and velocity scales the piston cylinder diameter (D) and the maximum piston velocity (U), respectively.However, the definition of this dimensionless number when applied to shear-thinning fluids is not trivial, because of the viscosity dependence on the shear-rate, and the reader is forwarded to each reference to check the specific definition that has been adopted in each case.It is worth noting that other works, not reviewed below, have also reported the formation of vortex rings in non-Newtonian fluids, though in different contexts.For example, Casanellas and Ortín [35] observed experimentally the formation of vortex rings, associated to flow reversal, in an oscillatory flow of micellar solutions inside a pipe, and Pillapakkam et al. [36] reported numerically the formation of vortex rings associated to a bubble rising in a viscoelastic fluid.
Palacios-Morales and Zenit [37] were among the first to study non-Newtonian vortex rings, using shear-thinning xanthan gum solutions with negligible elasticity and a power-law exponent in the range 0.48 < n < 0.61.Their experiments conducted for 138 ≤ Re ≤ 616 and 1 ≤ L/D ≤ 10, showed that shear-thinning (lower n) reduced the axial reach of the vortex rings.Moreover, such vortex rings expand abruptly in the radial direction as they slow down, in contrast to the gradual increase of diameter of Newtonian vortex rings.The effect of the stroke ratio on the vortex ring diameter was also negligible for L/D ≥ 4, suggesting the existence of a saturation value, as previously found for Newtonian fluids.Regarding vortex circulation, shear-thinning was shown to decrease its maximum value, whereas increasing the stroke ratio lead to a higher circulation.Depending on the stroke ratio, increasing Re did not always result in higher circulation, suggesting a non-monotonic relation between Re and vortex circulation.In the low Re range of the experiments, the authors were unable to observe a physical separation of the vortex ring from its tail (when it occurred, i.e. for L/D ≳ 4).Thus, it was not possible to ascertain about the existence of a formation number for purely viscous shear-thinning fluids, but a slower increase of the circulation with the stroke ratio was identified beyond L/D ≈ 6.Using a slug model approximation extended to power-law fluids, the authors showed that, under certain assumptions, such simple model could explain the decrease of the circulation for decreasing values of n.
The experiments by Bentata et al. [38] were also carried out with xanthan gum solutions (n = 0.65 and 0.76) of negligible elasticity, at even lower Re than in Palacios-Morales and Zenit [37], namely for 30 ≤ Re ≤ 300, and L/D = 3.This study confirmed some of the previous findings of Palacios-Morales and Zenit [37], for example, that shear-thinning reduces the axial reach of the vortex rings and leads to a sudden increase of diameter.In addition, it was shown that the late dissipation phase is Newtonian-like once the average strain-rate falls below the critical strain-rate delimiting the transition between the Newtonian plateau and the power-law regions of the shear viscosity exhibited by the xanthan gum solutions.The Carreau number controls such transition, therefore playing an important role in the type of dissipation phase.
In a subsequent study, Palacios-Morales et al. [39] replaced the xanthan gum solutions by polyacrylamide solutions, thus obtaining a viscoelastic shear-thinning fluid.The new set of experiments was performed for Re ≈ 20 and L/D = 4, corresponding to a Deborah number, Fig. 1.Schematic representation of (a) the idealized piston-cylinder setup and (b) the computational domain used to model it.The piston-cylinder setup is composed by a piston cased inside a long cylinder which pushes fluid into a large reservoir for a limited period of time, resulting in the generation of a (single) vortex ring which travels until total dissipation of the transported vorticity.The computational domain in (b) is axisymmetric and represents a wedge of the idealized piston-cylinder setup, where the piston motion is approached by a dynamic inflow boundary condition.The cylinder diameter is D, its length is 50D and Z B represents the length of the cylinder which penetrates inside the reservoir.The reservoir diameter is D R and the distance between the cylinder exit and the reservoir outlet is Z L .The drawings are not to scale.
F. Pimenta et al.De, of about 5. The main result of these experiments was that upon cessation of the piston translation, a secondary vortex, with negative vorticity, appeared in front of the primary vortex, annihilating it.Therefore, along with the annihilation of the primary vortex, the flow reversed its motion.A reduction in the circulation of the non-Newtonian vortex ring relative to the Newtonian case, at matched Re, was also reported by the authors, who hypothesized that the difference was stored as elastic energy.
More recently, Hegde et al. [40] conducted a series of experiments with viscoelastic polyacrylamide solutions at different concentrations and 3 ≲ L/D ≲ 7, where most of the tests were carried out at high Reynolds number (Re > 1000).When imposing the same piston velocity, and a necessarily varying Re among the different solutions, the authors observed a reduction in the circulation and kinetic energy of the vortex rings with increasing polymer concentration, as well as a late reduction in the propagation speed.The distribution of vorticity inside the vortices was markedly different between the Newtonian and non-Newtonian cases, being more disperse in the latter, but a formation number seemed to exist and remained independent of polymer concentration.In Re-matched experiments, the authors observed that some properties of non-Newtonian vortices, as enstrophy and peak vorticity, were not bounded between those of Newtonian vortices matching the low and high shear-rate viscosity (high and low Re, respectively) of the polymer solutions, which was attributed to differences in the vorticity distribution.The flow-reversal reported by Palacios-Morales et al. [39] was confirmed at low-Re experiments, although a secondary (negative) vortex was not reported, but the dye patterns clearly showed the unrolling process of the vortices while the flow reverses.
In a complementary study by the same group, Shashank and Sreenivas [41] looked at both polyacrylamide and polyethylene oxide solutions.The vortex rings in the mildly elastic constant viscosity polyethylene oxide solutions revealed a behavior very similar to Newtonian vortices and different from that of the strongly elastic and shear-thinning polyacrylamide solutions.The latter showed, again, a reduction of the propagation speed and an increased diameter, and the rollup of the ring was also inhibited and reversed at a later stage.The Reynolds number was high enough to trigger Widnall instabilities [11] in the Newtonian and polyethylene oxide solutions, which were absent in the polyacrylamide solution.However, it is possible that the low elasticity of the polyethylene oxide solution used (resulting in Deborah numbers below 1) was not high enough to highlight elastic effects.
These experimental studies reveal clear differences between Newtonian and non-Newtonian vortex rings, as well as some similarities.Moreover, some behaviors were cross-validated by different works, as for example the reduction of the axial reach in shear-thinning fluids [37,38] and flow reversal in viscoelastic fluids [39][40][41].Nevertheless, whereas the physics behind some observations can be reasonably explained based on simple and occasionally empirical arguments, others require a more complex analysis.Such analysis is sometimes challenging in experimental work for various reasons.For example, some variables may not be directly accessible, like polymer-induced stresses, it may be difficult to separate shear-thinning from elasticity, to change independently a given parameter while keeping others constant or to reach the limit of those parameters.Much of these limitations can be surpassed with the use of numerical simulations.In this work, we use computational fluid dynamics to investigate a variety of non-Newtonian vortex rings, with particular attention to the post-formation phase.
The numerical simulations are performed with a finite-volume method for Newtonian, power-law, Carreau and Phan-Thien-Tanner (PTT) fluids, thus covering different rheological behaviors (constant and shear-thinning viscosity; inelastic and elastic fluids).The Reynolds number is kept below 10 3 to keep our results framed with experiments and to avoid the onset of inertial instabilities, which justifies using an axisymmetric approach.The main purpose of this work is to reproduce and propose a physical explanation for some of the laminar non-Newtonian vortex ring behaviors reported in the experimental literature [37][38][39][40][41], while also extending their range of parameters.The vortex ring behavior in Newtonian fluids is taken as the basis of comparison and the rheology of the different fluids is the main effect under analysis.To the best of our knowledge, this is the first numerical attempt to model non-Newtonian vortex rings generated in a piston-cylinder setup.We note, however, that this study is not exhaustive regarding the parameters and constitutive equations under analysis, and there is still plenty of room for future works on this subject, particularly those related with viscoelastic vortex rings, which showed a particularly complex behavior and can be richer and represented by more complex constitutive equations than those used here.
The remainder of this work is organized as follows.The problem under analysis is formulated in Section 2, which provides some physical background and details about the numerical method implemented to solve it.Section 3 reviews several analytical solutions for laminar Newtonian vortex rings, that will be used to compare and extend the numerical results analyzed in Section 4. Finally, the main conclusions are discussed and summarized in Section 5.

Conceptualization
Vortex rings can be generated in different apparatus, both in experiments [7,42,43] and in numerical simulations [22].Each method has its own advantages and drawbacks, and can influence the vortex ring properties [22,42].Although relevant, studying the different methods is out the scope of this work.Instead, for our purpose we selected a setup mimicking the one used in the experiments with non-Newtonian vortex rings [37][38][39][40][41], the piston-cylinder arrangement (Fig. 1a), which is perhaps the most commonly used configuration.
In a piston-cylinder system, the piston cased inside a cylinder drives fluid into a reservoir for a limited time interval.The reservoir is typically large compared to the cylinder, such that its walls have a negligible effect on the vortex rings.Initially, the piston can be located at different positions inside the cylinder and it will be assumed here that this initial position is far enough from the cylinder exit, in order to have a negligible effect on the vortex ring properties.This simplification allows the use of a static mesh, since the piston motion can be adequately approximated by a dynamic inflow boundary condition.
The formation phase of the vortex ring corresponds to the time interval during which the piston moves inside the cylinder.The piston stroke is denoted by L, such that the stroke ratio (L/D) is the ratio between the piston stroke and the piston diameter (D).During the stroke phase, the shear-layer adjacent to the wall rolls-up and feeds vorticity to the vortex ring under formation at the cylinder exit.When the piston stops, the vortex ring detaches from the cylinder and propagates due the self-induced velocity.The vortex ring stops translating once all the vorticity transported is dissipated.In this work, we will focus the attention on the post-formation phase, i.e. the period of time just after the piston stops.

Geometry and mesh
The geometry of the computational domain adopted to model the piston-cylinder system is depicted in Fig. 1b with a cylindrical system of coordinates (r, φ, z) having its origin on the axis, at the cylinder exit plane.Since axisymmetry can be assumed for the flow conditions under investigation, the computational domain corresponds to a wedge of the idealized piston-cylinder setup (Fig. 1a), over which there is a single set of cells in the azimuthal direction.The cylinder has diameter D and length 50D, which is long enough to avoid entry effects, in agreement with the simplifications assumed in the previous section to model the piston translation.The length of the cylinder immersed in the open reservoir is denoted by Z B and it is assumed that the cylinder wall has a negligible thickness.The reservoir is open on its right side (outlet), its F. Pimenta et al. diameter is D R and Z L represents the distance between the cylinder exit and the reservoir outlet.
For the simulation of unbounded vortex rings, and unless otherwise stated, the following geometric parameters were used (denoted as geometry 1): Z B = 25D, Z L = 100D and D R = 100D.In order to reduce the computational time, these sizes were shortened to obtain the axial reach of the dye front more accurately, as well as to simulate viscoelastic fluids.In such cases, the following geometric parameters were used (denoted as geometry 2): Z B = 4D, Z L = 50D and D R = 24D.The reservoir outlet is far enough from the cylinder exit in both geometries, so that it has a negligible effect on the dynamics.The smaller values of Z B and D R in geometry 2 impose some degree of confinement, although it must be noted that any value of D R will likely affect the vortex dynamics once the vortex diameter exceeds ~DR /2, which should be accounted for during the analysis.As a reference, the following range of geometric parameters were used in non-Newtonian experiments in the literature [37][38][39][40][41]: 3.5 < Z B /D < 8, 20 < Z L /D < 60 and 9 < D R /D < 15 (the reservoirs used in the experiments typically have a rectangular/square cross section, such that D R represents in fact the largest dimension of the cross section).Thus, the experimental conditions typically impose a higher degree of confinement than any of the two sets of geometric parameters tested.
A preliminary mesh study was carried out to assess the effects of the finite resolution of the computational grids on the results (Appendix A).The details of the meshes used in this study with inelastic (Newtonian, Carreau and Power-law) and viscoelastic (PTT) fluids are presented in Table 1.The mesh characteristics in this table are for geometry 1, but the grid resolution for geometry 2 is identical, since this geometry is simply obtained by truncating the domain of geometry 1.The cells are uniformly distributed in the radial direction in region r ̃< 0.5 (where r ̃= r/ D), whereas they are compressed toward r ̃= 0.5 outside this region.There is also compression of cells in the axial direction toward z̃= 0 (where z̃= z/D).The smallest cell is always located at the lip of the cylinder, next to the wall, and its dimensions (normalized by D) are presented in Table 1.The number of grid divisions over different sections and directions of the mesh is also presented in the table, as well as the total number of cells (note that there is also grid for z̃< 0).Meshes M2 were used throughout this work according to the error analysis (Appendix A).

Governing equations
The equations solved in this work to simulate non-Newtonian vortex rings are mass conservation, momentum balance and the rheological constitutive equation for the fluid extra-stress tensor.
The mass conservation equation for isothermal, incompressible flows is written as where u is the velocity vector.Momentum balance can be expressed by ρ ( ∂u ∂t where t is the time, p is the pressure, τ' represents the extra-stress tensor and ρ is the fluid density.For the general case of a viscoelastic fluid, the extra-stress tensor can be decomposed into solvent (τ s ) and polymeric (τ) contributions, such that τ' = τ s + τ.A Generalized Newtonian fluid (GNF) can be obtained as a special case by setting τ ¼ 0. In this work, we investigate four different constitutive equations: Newtonian [44], power-law [44], Carreau [44] and PTT [45] (in its simplified form, with a zero second normal stress difference in homogeneous shear flow) fluids.The expressions of τ s and τ for each model are displayed in xy /γ) of these fluids in homogeneous shear flow is illustrated in Fig. 2.
A Newtonian fluid has constant viscosity η s (Fig. 2) and is inelastic.A power-law fluid is also inelastic, but has a shear-thinning viscosity controlled by the consistency index (K) and the power-law exponent (n) -lower values of n stand for more shear-thinning fluids (for n > 1 the viscosity shows a shear thickening behavior, but this type of fluids will not be analyzed here).The viscosity grows unbounded to +∞ as γ→0 and tends to 0 as γ→ + ∞(Fig.2).The Carreau model can be seen as a bounded power-law since the viscosity is bounded between η 0 at low strain-rates (the so-called first Newtonian plateau of viscosity) and η ∞ at high strain-rates, although the lower bound is neglected in this work through the use of a high η 0 /η ∞ ratio (Fig. 2).For intermediate strainrates, the Carreau model behaves as the power-law model with the same value of n and the power-law consistency index K = η 0 Λ n-1 .The Carreau model is typically a better option to model real fluids, as a lowshear Newtonian plateau is exhibited by many shear-thinning fluids [44].
The PTT model is one of several constitutive equations available to represent fluids that are simultaneously shear-thinning and elastic.The solvent component of the extra-stress tensor for this model is Newtonian-like, with a constant viscosity (η s ).On the other hand, the polymeric stresses need to be obtained implicitly from the solution of the partial differential equation in Table 2, which is controlled by the relaxation time (λ), the extensibility parameter (ε) and the polymer viscosity coefficient (η p,0 ).In the limit ε → 0, a PTT fluid loses its shearthinning behavior and behaves as an Oldroyd-B fluid.Note that we are considering a simplified form of the PTT model for the polymer stress

Table 1
Characteristics of the computational grids for geometry 1 (for geometry 2 the resolution is identical) used to assess the mesh dependency of the numerical results.
because the time derivative of the stress tensor is the upper-convected derivative, in contrast to the more general form using the Gordon-Schowalter time derivative [45].
The choice of the above constitutive equations, among the many available to simulate non-Newtonian fluids, was primarily based on the rheology of the fluids used in the experiments with non-Newtonian vortex rings [37,38,40,41], which were reviewed in the Introduction.
In particular, the rheology of inelastic or weakly elastic xanthan gum solutions can be well modeled by the power-law (at high strain-rates) and Carreau constitutive equations, whereas the PTT model is suitable for viscoelastic polyacrylamide solutions.The rheology of the viscoelastic polyethylene oxide solution used by Shashank and Sreenivas [41] could be approached, for example, by an Oldroyd-B model, but this particular class of constant-viscosity fluids will not be addressed in this study.These constitutive equations are adequate to model the rheological behavior of the above-mentioned fluids, even though the choice is not unique.
To help visualize and characterize vortex rings, the equation for the transport of a passive scalar (a dye) was added to the simulations.Denoting the concentration of the passive scalar by C, its corresponding transport equation is

∂C ∂t
where D eff is the effective diffusion coefficient, here assumed independent of the fluid rheology and flow conditions.

Numerical method
The set of governing equations was solved in the finite-volume framework of OpenFOAM® using rheoTool [46,47].The solver was described and validated in previous works [46,48], thus we will abbreviate our discussion on that subject and limit this section to the most important topics.
Convective terms in the governing equations are discretized with the CUBISTA high-resolution scheme [49], whereas centered-differences are used for interpolation in the remaining terms.Time-derivatives are discretized with the implicit Adams-Moulton scheme.The numerical method is formally second-order accurate in both space and time.The SIMPLEC algorithm [50] is employed for pressure-velocity coupling, since the equations are solved in a segregated approach [46].To retain the formal second-order accuracy of the time-derivative discretization scheme [46], the governing equations are solved twice inside each time-step.The time-step is updated dynamically during the simulation, using a double criterion based firstly on a maximum local Courant number (Co max ), which typically controls the time-step in the formation phase, and secondly on an advective time-scale, which typically assumes control in the post-formation phase.Three pairs of criteria were tested for inelastic (Δt ̃1 = min{0.25,Δt ̃Co,max=0.1 }, Δt ̃2 = min{0.5,Δt ̃Co,max=0.2} and Δt ̃3 = min{1, Δt ̃Co,max=0.4}) and viscoelastic (Δt ̃1 = min{0.05,Δt ̃Co,max=0.1 }, Δt ̃2 = min{0.1,Δt ̃Co,max=0.15} and Δt ̃3 = min {0.2,Δt ̃Co,max=0.2}) fluids (where t ̃= tU/D), in order to evaluate the results dependency on the temporal resolution.The error analysis is presented in Appendix A, leading to the choice of Δt ̃2 for both kinds of fluids.
A commonly called log-conformation procedure [51] is applied to solve the constitutive equation for the PTT model, which is currently a standard and widely used technique to improve the numerical stability of viscoelastic solvers [46].For viscoelastic flow calculations with the PTT model, the polymeric stresses are coupled to the flow velocity through the method described in Pimenta and Alves [46].
Regarding boundary and initial conditions, standard no-slip, nopenetration and no-flux conditions were assigned to the walls.For viscoelastic fluids, the polymeric stresses were linearly extrapolated at the walls [46].In order to satisfy continuity and flow incompressibility, the reservoir outlet is an outflow, where all the variables are assumed to be fully developed, except pressure, which is fixed (the large value of Z L ensures a minimal effect of this outflow boundary).Note, however, that the outflow could instead be assigned to any of the other boundaries of the reservoirin the experiments, one of the sides of the reservoir is left open to the atmosphere.As aforementioned, the piston motion inside the cylinder is approximated by a time-varying boundary condition assigned to the inlet velocity.A simple form of imposing a stroke ratio L/D is through the definition of a rectangular velocity program, i.e. the inlet velocity is kept fixed at u z = U for t ≤ L/U and zeroed right after.However, the discontinuity in time is not convenient from a numerical perspective as, for example, some quantities would go to infinity in the first time-step.Moreover, a perfect rectangular program is also unfeasible in the lab, as there is always a finite amount of time during which the piston accelerates from 0 to U (at the beginning) and decelerates from U to 0 (at the end).James and Madnia [25] proposed a smoother inlet velocity program that preserves a predefined stroke ratio and preserves continuity in time, where U is the maximum velocity attained during the program and t a and t b are variables controlling the shape of the program.If t b = L/U, the program tends to a rectangular function as t a → 0. In this work, we defined t a = t b /25 = L/(25U), which is small enough to have negligible acceleration/deceleration phases, but sufficient to keep the inlet velocity continuous in time.One can still consider with good approximation that u z (0 ≲ t ≲ L/U) = U and u z (t ≳ L/U) = 0.Although varying in time, the inlet velocity profile is uniform in space, i.e. over the whole inlet boundary, to simulate the piston displacement.The polymeric stresses are considered null at the inlet.At the beginning of the simulations, the fluid is quiescent, the stress field is null and the cylinder is entirely filled with the passive scalar (dye), at a concentration C 0 .It is worth noting Note that these parameters were selected arbitrarily, for the sake of the example (to show the viscosity dependence on the strain-rate), and do not aim to reproduce the fluid rheology of any particular case simulated.SI units without multiples or submultiples are assumed in the parameters of the models.
that the computational domain adopted in the simulations displays a singularity over the curve (r, z) = (D/2, 0), i.e. over the entire perimeter of the wall at the cylinder exit.However, neither the numerical stability, nor the accuracy of global parameters were significantly affected for the cases tested.The matrices resulting from the discretized set of equations are solved with a pre-conditioned Bi-CGStab solver [52].For the particular case of the pressure (continuity) equation, the parallel algebraic multigrid implementation of Hypre (BoomerAMG) [53] is used for pre-conditioning (it is available in rheoTool through the interface built by Pimenta and Alves [48]).The remaining matrices were pre-conditioned with the standard geometric algebraic multigrid algorithm of OpenFOAM®.

Dimensionless numbers
The dimensionless numbers controlling the generation and evolution of vortex rings depend on the specific rheological model.The stroke ratio and Reynolds number are two dimensionless numbers common to all fluid models, but the definition of the latter for fluids with variable viscosity depends on the specific rheological model.
The generalized Reynolds number for power-law and Carreau fluids can be defined as ) n ρU 2− n D n k [38], where k = K for the power-law model and k = η 0 Λ n− 1 for the Carreau model (the Reynolds number for Newtonian fluids arises as a special case for n = 1 and k = η s ).This is one of various possible definitions [54] and it is based on a constant Poiseuille number for laminar Poiseuille flow of power-law fluids.Its extension to Carreau fluids is an approximation based on the same constant Poiseuille number in the same pipe flow and on the similarity between the two models at high shear-rates, for small values of η ∞ .In this work, we set η 0 /η ∞ ≈ 10 10 , which allows in fact neglecting the effects of a viscosity plateau at high shear-rates.For Carreau fluids, a second Reynolds number can still be defined based on its viscosity at the low shear Newtonian plateau, Re 0 = ρUD η 0 .The power-law exponent (n) is a second dimensionless number common to both shear-thinning inelastic models.For the Carreau model, there is an additional dimensionless number to be taken into account, the Carreau number, defined D [38].The Carreau number compares a characteristic flow shear-rate with the characteristic shear-rate at the transition between the Newtonian plateau and the power-law region of viscosity decay (~1/Λ).The definition used here is again based on fully-developed Poiseuille flow.Low values of Cu correspond to a Newtonian-like behavior, whereas high values correspond to a power-law-like behavior.
The shear viscosity of a viscoelastic PTT fluid is also dependent on the strain-rate (Fig. 2), but a simple expression for a generalized Reynolds number is not straightforward.There are known analytical or semi-analytical solutions for the flow of PTT fluids inside pipes, from which a Reynolds number could be deduced based, for example, on the viscosity at the wall [55].One could alternatively fit (approximately) a Carreau model to the viscosity curve of the PTT model in homogeneous shear flow (Fig. 2) and use the above definition of Re G for Carreau fluids.However, for simplicity, the Reynolds number was defined based on the zero shear-rate viscosity (η 0 = η s + η p,0 ), Re 0 = ρUD The Péclet number is a dimensionless number relating the advective and diffusive transport of the passive scalar and is defined as Pe = DU/ D eff .This dimensionless number is O (10 10 ) in all the simulations, i.e. mass diffusion is made negligible compared to mass convection to obtain a well-defined interface, which does not necessarily follow the vortex bubble once the vorticity transport becomes itself diffusion-dominated (here related to momentum diffusivity).Note that the Péclet number does not have direct influence on the vortex dynamics (it may only affect our perception when analyzing the dye patterns).
In the remaining of this work, the following characteristic scales are used to normalize variables: D for length, U for velocity, D/U for time, U/D for vorticity, UD for circulation and C 0 for the passive scalar concentration.All dimensionless variables are over-scripted with a tilde, such that, for example, z̃c = z c /D and Γtot = Γ tot /(UD) represent the dimensionless axial coordinate of the vortex center and the dimensionless total circulation, respectively.

Vortex ring characterization
In this section, we identify and define the set of properties used to characterize the vortex rings.
The vortex center in a given azimuthal plane is denoted by X c = (r c , z c ) and corresponds to the point of maximum vorticity magnitude inside the vortex bubble (the vortex diameter is D c = 2r c , according to the axisymmetry approximation).This is one of several definitions reported in the literature for the vortex center, since the definition and characterization of vortices is still a subject of current debate and investigation [56].However, although there are certainly differences between the various definitions of X c , the global trend can be often captured independently of the definition used.The velocity of the translational motion of the vortex rings is simply defined as U c = ∂z c /∂t.The axial position of the dye front (z d ) is defined as the z-coordinate (linearly interpolated) over the axis where the dye concentration reaches half its maximum value near the dye-no dye interface.For a dye concentration ranging between 0 (no dye in the computational cell) and 1 (the computational cell is full of dye), the interface is typically located where the concentration is equal to 0.5.
The total circulation (Γ tot ) is another flow property used to characterize vortex rings and corresponds to the surface integral of the azimuthal vorticity over a given azimuthal plane, outside the cylinder and only considering the faces with positive vorticity.Mathematically, the total circulation is given by where ω ϕ = n ϕ ⋅(∇ ×u) represents the azimuthal vorticity and A p corresponds to the group of faces in the azimuthal plane satisfying the aforementioned conditions.There is usually a distinction between total circulation and vortex circulation in the literature, where the latter corresponds to the vorticity integral associated with the vortex region only, i.e. excluding any possible trailing wake of vorticity.However, the laminar flow conditions simulated in this work, associated to small stroke ratios, made the identification of the vortex region not obvious.The same problem was recognized by Palacios-Morales and Zenit [37], who used the Q-criterion (simply put, in which the vorticity is only accounted for when its magnitude exceeds the magnitude of the rate of strain tensor) to delimit the vortex region.Nonetheless, we observed that the Q-criterion is too restrictive and excludes some of the vorticity inside the vortex bubble, only accounting for the vorticity inside the vortex core.Thus, given the absence of an objective and consensual criterion to delimit the vortex bubble, we preferred not to compute the vortex circulation and used instead the total circulation.

Analytical solutions for Newtonian fluids
In spite of its apparent complexity, vortex rings are sufficiently simple vortical flow structures amenable to analytical handling.There is a vast theory developed for this type of vortices in Newtonian fluids F. Pimenta et al. [19], from which some (exact) analytical and asymptotic solutions exist.While it is out of scope of this work to verify their applicability to our numerical results, some will be used to check the consistency of the simulations with Newtonian fluids, which represent the basis of comparison for non-Newtonian fluids.Thus, it is worth to briefly introduce some analytical solutions framed with this work, i.e. valid for the post-formation phase of laminar (low Re) Newtonian vortex-rings, which will be mainly used to assess the dynamic evolution of some vortex properties.The interested reader is directed to the original works for more details.
The most comprehensive formulas describing the transient behavior of low-Re vortex rings are perhaps those developed by Fukumoto and Kaplanski [14] and Fukumoto [13], which are capable to describe both the early and late time dynamics.They are based on hypergeometric functions, but can be approached with good accuracy by simple asymptotes.Incidentally, these asymptotes are corrections to analytical solutions developed years before [16][17][18][19]57] and that will be used in this work for comparison purposes due to their inherent simplicity.The early-time (υt/R 2 0 << 1) evolution of the translational motion velocity of the vortex rings can be approached by Saffman's solution [18,19], where υ is the kinematic viscosity, Γ 0 is the initial circulation entrained by the vortex ring, R 0 is the initial vortex ring radius and ∊ represents the error arising from the omission of higher-order terms.On the other extreme, i.e. for late times (υt/R 2 0 >> 1), when most of the vorticity has already been dissipated, the translational velocity can be approximated by Rott-Cantwell's solution [16,17], where I 0 represents the specific impulse initially transmitted to the vortex ring.Kambe and Oshima [57] developed the following asymptotic solutions for the vortex radius and circulation in the late-time limit t → +∞ (vanishing Re), Eqs. ( 7) and ( 8) also correspond to the late-time behavior of Saffman's solution for the matured stage [14,18], where k and k' are constants to be fitted to the data.The integration of U c up to t →+∞ leads to an estimate of the maximum axial reach of a vortex ring [14], considering k ≈ 10.15 and k' ≈ 8.909 [14].

Results
The results will be first presented for power-law fluids, followed by Carreau and PTT fluids.

Power-law fluids
The simulations with power-law fluids were carried out for powerlaw exponents between 0.4 and 1 (Newtonian), covering a wide range of shear-thinning fluids that can be experimentally tested.The generalized Reynolds number was kept between 10 and 1000, as in most non-Newtonian laminar flow experiments [37,38].However, it is worth noting that the experiments by Maxworthy [58] showed that Newtonian vortex rings were stable only up to Re = 600, i.e. the vortices lose axisymmetry beyond that Re.The Newtonian simulations in this work were for Re ≤ 200, with the exception of the results plotted in Fig. 9, which contains data up to Re = 900, although the corresponding simulations were only run for a short time (t ̃≈ L/D).For shear-thinning fluids (n < 1), there are two simulations in the range 600 < Re G ≤ 1000, obtained for fluids with n ≤ 0.6 (the same exceptions apply to Fig. 9).Thus, most of the simulations fall in the stable regime identified in Maxworthy's experiments [58], where the axisymmetry assumption should be valid.Nevertheless, interpretation of our results outside this range should be done cautiously, although we consider that the increased viscous stabilization effect over time (owing to a decrease of the strain-rate) should keep the resulting vortex rings axisymmetric and stable.The stroke ratio was held fixed at L/D = 4, which is also within the range observed in experiments [37,38].Therefore, the independent parameters under study are n and Re G .

Kinematics
The trajectories described by the vortex rings at different Re G are plotted in Fig. 3.In general, the trajectories exhibit two different phases.Initially, the vortex rings move essentially in the axial direction, with only minor changes in their radial coordinate.As Re G is increased, the initial vortex ring radius seems to converge to an approximately constant value ( R0 ≈ 0.77) independent from both n and Re G , and the same seems to happen to the trajectories, i.e. the vortex rings follow the same path independently of these two dimensionless numbers.This is in qualitative agreement with experimental results [37].
In the second phase, when the total circulation enters a power-law decay (see Section 4.1.4),the vortex bubble grows both radially and axially.This expansion of the vortex bubble, more abrupt as n decreases, is characterized by the radial motion of the vortex center away from the axis (r̃c increases), and a less intense axial diffusion of the vortex center, which first diffuses backward and then forward.For shear-thinning fluids (n < 1), there is also an inward motion (r̃c decreases) of the vortex center in the transition between the two phases, as clearly evidenced in Fig. 3a and Video 1 (Supplementary material).The theories developed for idealized Newtonian vortex rings [14] do not predict a final backward motion of the vortex center in the axial direction.Such behavior, observed in the numerical simulations is attributed to the boundaries effect (including the cylinder outer wall) and to the specific criterion used to define the vortex center, which is more sensitive to the diffusion of the vorticity peak inside the vortex bubble than other criteria, as for example Helmholtz's definition [59] of the vortex centroid (based on an integral of vorticity).It is worth noting that due to the unboundedness of the power-law viscosity at low shear-rates, i.e. as time advances, the background flow develops spots of noisy vorticity, which hinders the automatic tracking of the vortex rings based on a vorticity criterion.In general, this did not restrict our analysis, since these spots appeared at late times, well after the radial coordinate of the vortex rings exceeded 15D in geometry 1, which was defined as the (safe) maximum radial distance to obtain results unaffected by the side walls of the reservoir.The exception was the set of simulations carried out for n = 0.4, as evidenced by the short trajectories depicted in Fig. 3.
The radial coordinate of the vortex rings is plotted in Fig. 4 as a function of time, for n = 0.6 and n = 1.Although this is not clear from this figure, the visualization of the corresponding video (Video 3, Supplementary material) shows that the radial displacement of the vortex center corresponds to a radial increase in the size of the vortex, i.e. to vortex growth.In general, one can identify three different phases for both fluids, but the three are not present for all Re G and their start/end occur at different instants of time for each case.Taking the curve at Re G = 900 (n = 0.6) as example, the initial phase holds approximately until t ≈ 8 and is characterized by an almost constant vortex radius.The intermediate phase is observed during the period 8 ≲ t ̃≲ 50 and displays a radius growth approximately proportional to tm0 , with m 0 ≈ 0.14.We    verified that this power-law exponent (m 0 ) is nearly constant at sufficiently high Re G for the other values of n simulated (including the Newtonian case, as shown in Fig. 4b), suggesting that in this phase the radius (or vortex) growth rate is independent of shear-thinning and saturates with Re G .At low Re G , there is no clear separation between the initial and intermediate phases (Fig. 4).The final phase (t ̃≳ 50) is present in all curves plotted in Fig. 4 and also exhibits power-law growth (r c ∝ tm1 ), but significantly faster than the previous one.The final powerlaw exponent in this phase depends on the shear-thinning level, being lower for the Newtonian fluid (m 1 ≈ 0.5), than for the power-law fluid with n = 0.6 (m 1 ≈ 1.8).We have further estimated m 1 ≈ 0.8 for n = 0.8 (the trajectories for n = 0.4 are too short to obtain an accurate fit), but the corresponding plot is not shown for conciseness.Notwithstanding the dependence on n, this power-law exponent (m 1 ) shows little variation with Re G for both fluids plotted in Fig. 4, which reflects the dominant role of viscous effects over inertia in this phase.For the Newtonian fluid the power exponent in the late phase (m 1 ≈ 0.5) compares well with Eq. ( 8).The inward motion of the vortex center in the transition to the last phase can be clearly seen in Fig. 4a for the shear-thinning fluid.
The translational velocity of the vortex rings center is plotted in Fig. 5 as a function of time for n = 0.6 and n = 1.The initial velocity is slightly higher for the Newtonian fluid and is close to half the piston velocity (the value given by the slug model prediction for Newtonian fluids [60]), in agreement with the experiments of Palacios-Morales and Zenit [37].The velocity starts decaying right after the vortex rings detach from the cylinder lip in the Newtonian case, whereas in the shear-thinning fluid the velocity reaches a maximum a few instants later, after which it starts decaying.The velocity enters a power-law decay in the Newtonian case, and the corresponding power-law exponent shows little dependency on Re G .All viscous vortex rings experience a vanishing local/instantaneous Re as t ̃→ +∞, which explains the Re G -independent behavior in Fig. 5b, although vortex rings with lower Re G develop this decay period earlier.The final decay is, however, slightly faster than the theoretical prediction of Eq. ( 7), which is also depicted in Fig. 5.In contrast, the velocity decay for the shear-thinning fluid plotted in Fig. 5a is clearly faster and does not seem to exhibit a power-law evolution.In fact, the curves at different Re G can be better described (at least locally) through expressions like , and this could be generalized to the other values of n not shown in Fig. 5.This faster decay was already expected due to the enhanced viscous effects associated with the combined (and interrelated) velocity decrease and viscosity increase over time (there is only a velocity reduction in the Newtonian case).

Dye patterns
The dye patterns developed by the vortex rings at three different instants of time are exemplified in Fig. 6 for n = 0.6 and n = 1, at Re G = 100, and can be seen in Video 2 (Supplementary material) for other combinations of n and Re G and over their whole lifecycle.The transient and final dye patterns are naturally different for each pair (n, Re G ), but, more importantly, they all are qualitatively similar in the sense that they all show the formation of a mushroom-like structure that moves continuously forward in the axial direction.These structures (dye patterns) are wider, more elongated and more detached from the axis, and travel larger distances as n and/or Re G increase.

Vorticity
The dynamic evolution of the vorticity patterns is also plotted in Fig. 6 for the same two cases (n = 0.6 and n = 1, at Re G = 100), but more combinations of n and Re G can be seen in Video 3 (Supplementary material).In general, the patterns are qualitatively similar between the different conditions, as they all show a single and well-defined region where the vorticity is more intense, as well as the absence of any significant vorticity trailing wake.The isocontour lines in the vortex core are nearly elliptical in the early instants of the post-formation phase, with the major axis slightly tilted with the z-axissome models of Newtonian vortex rings predict a Gaussian vorticity distribution, i.e. circular isocontours in a cut plane [2,8], which have been observed experimentally in the vortex core (e.g.[2]).However, there is no symmetry around the vortex center (point of maximum vorticity) in both radial and axial directions.The vortex bubble grows continuously over time, although this is more evident near the end of life of the vortex rings, when the vorticity transport becomes diffusion-dominated.The video also shows the growth of a counter-rotating vortex near the edge of the cylinder exit, that interacts with and has some influence on the main vortex, especially near the end of life of vortex rings with a low axial reach.
It is interesting to note in Fig. 6 and Video 3 (Supplementary material) that there is a clear relation between the evolutions with time of vorticity and dye concentration while the transport of both fields is dominated by convection.Once the vorticity transport becomes diffusion-dominated (late-time behavior), the vorticity patterns keep changing by diffusion (controlled by an increasingly large kinematic viscosity for shear-thinning fluids), whereas the dye patterns become frozen due to the very large Péclet number of order O(10 10 ) in all simulations.
Comparing the numerical results with the experiments by Bentata et al. [38] and Palacios-Morales and Zenit [37], both with  shear-thinning (Carreau) fluids, we observe a qualitative agreement in the vorticity patterns.Palacios-Morales and Zenit [37] reported a discontinuity between the tail and the core of the vortices for Newtonian fluids, which was apparently absent in shear-thinning fluids at similar Re G .However, this observation was made for a high stroke ratio (L/D = 8), for which a trailing jet of vorticity is noticeable and expected.

Total circulation
The evolution of the total circulation normalized by the maximum value attained during the whole vortex lifecycle ( Γtot,max ) is plotted in Fig. 7 for one of the shear-thinning fluids (n = 0.6) and for the Newtonian reference (n = 1).Generally speaking, the decay of circulation is slower as inertia increases and faster for the power-law fluid.However, after an initial period, the vorticity enters a power-law decay with an exponent nearly independent of inertia (for sufficiently high Re G ), but sensitive to shear-thinning.In fact, the power-law decay exponent is approximately − 1 for n = 1 (in agreement with the theoretical prediction of Eq. ( 9)), − 1.6 for n = 0.8 and − 3.59 for n = 0.6.The data for n = 0.4 was insufficient to obtain an accurate fit.
Denoting α as the power-law exponent for the late-time growth of r ̃c with t ̃(Fig.4) and χ as the power-law exponent for the decay of the total circulation with t ̃(Fig.7), then it follows naturally that Γtot / Γtot,max ~ r ̃cχ/α for the period of time during which the two relations hold true.
From Eqs. ( 8) and ( 9), for a Newtonian fluid χ/α = − 2 is expected and this is also confirmed from Figs 4b, 7b and 8b.In addition, this ratio also holds true for n = 0.6 and 0.8, as shown in Figs.8b and c.The data collected at different Re G collapse relatively well for sufficiently high values of Re G .However, there is a clear deviation from this scaling for n = 0.4, although a power-law relation between Γtot / Γtot,max and r ̃c still seems valid in this case (Fig. 8a).
The total circulation reaches a maximum ( Γtot,max ) close to the instant when the piston motion stops and the relation between Γtot,max and the power-law exponent of the shear-thinning fluids is plotted in Fig. 9 as a function of Re G .Independently of n, the variation is nonmonotonic, with the total circulation fed to the flow reaching maxima in the range 100 ≤ Re G ≤ 150.For a given Re G , Γtot,max decreases with increasing shear-thinning, but seems to tend to an n-independent asymptote (~3) as Re G increases.This is probably due to the increase of   the ratio of diffusion over convective time scales with Re G inside the cylinder, which results in flatter and more homogeneous (less developed and also less dependent on n) velocity profiles of the fluid leaving the cylinder during the limited period of time when the piston is moving.It is interesting to note that although Γtot,max starts to decrease beyond a given Re G , the axial reach increases continuously with Re G (Fig. 10) due to the scaling of z̃c being approximately proportional to the product Γtot, max Re G , as will be shown in the next section .The numerical results are in agreement with Palacios-Morales and Zenit [37], where it was shown both experimentally and theoretically (through the application of the slug-flow theory to power-law fluids) that shear-thinning reduces the circulation injected in the flow.

Axial reach
The final stages of axial diffusion of vorticity hamper the definition of an axial reach based on z̃c, which has been replaced by two alternative parameters.The first is z̃c ,1 and represents the maximum axial reach of the vortex center before any significant change of its trajectory (e.g.backward motion).This new parameter is identified in Fig. 3b through circles and plotted in Fig. 10a for different combinations of (n, Re G ).It is noteworthy that z̃c ,1 could not be defined for some cases, as for example for n = 1 and Re G = 10, where a first peak of axial reach is not obvious.The second parameter proposed is the final axial position of the dye front  measured on axis, i.e. z̃d ,f = z̃d(t ̃→∞).The front of the dye (interface) on the axis does not share, in general, the same axial coordinate as the vortex center (z̃d ∕ = z̃c), but there is a relation between both variables (not necessarily constant over time).In general, when one of the variables moves forward in time, the other follows the same trend, except, possibly, near the end of life of the vortices, since the dye concentration is very weakly diffusive compared to the vorticity.The final position of the dye front is arguably more suitable to characterize the final state of the vortex rings, since z̃c becomes undefined once all the vorticity is dissipated (t ̃→ +∞).
The two quantities are plotted in Fig. 10.z̃c ,1 increases with inertia and, for the same Re G , it decreases with shear-thinning (Fig. 10a).This agrees with experiments [37,38] and can be attributed to the shear-dependent viscosity increasing on moving downstream.Indeed, as the vortex rings translate axially, their enstrophy decays, the vortex rings slowdown and expand as they entrain fluid, and the strain-rate decreases leading to an increase of the viscosity which feeds this positive feedback loop.The same trend is observed for z̃d ,f , as shown in the inset of Fig. 10b, where z̃d ,f scales approximately with Re G 0.7 .Eq. ( 10) predicts a linear scaling between the axial reach and Re G for Newtonian vortex rings and Re G -independent R0 and Γ 0 .However, in general, both R0 and Γ 0 depend on n and Re G , such that a linear relation should not be expected when plotting the axial reach as a function of Re G .When extending this analysis to the numerical results, one needs further to consider that Eq. ( 10) does not account for the axial distance travelled during the formation phase (t ̃< L/D), that originates from the piston motion.In an attempt to test the possible compliance of Eq. ( 10) with the numerical results, some approximations were made.We consider that z̃d ,f is proportional to the effective axial reach of the vortex, that Γtot,max is proportional to Γ 0 and that the piston-contributed axial distance is given by z̃d ,0 ≈ z̃d (t ̃= L/D), which represents the inflation of the vortex ring as the piston moves.In addition, because it was seen in Section 4.1.1 that R0 seems to converge to an n-and Re G -independent constant value as Re G increases, R0 is here assumed to be constant and its effect ignored.Under these assumptions, (z̃d ,f -z̃d ,0 ) is plotted against Γtot,max Re G in Fig. 10b, where it can be seen that a power-law scaling with exponent ~1 (a linear relation) is observed at the high-Re G limit (where the assumption of constant R0 holds approximately).Therefore, one may conclude that Eq. ( 10) is also likely valid for power-law fluids after replacing the constant coefficient by an appropriate function of n.Cu and Re 0 evolve inversely for n < 1.Thus, for fixed n and Re G , a Newtonian-like behavior is obtained in the limit of high Re 0 or (Cu→0), whereas in the low Re 0 limit (Cu→+∞) one can expect a power-law-like behavior.In particular, for the set of simulations carried out, we have the following correspondence at the lowest Cu simulated in each case: The vortex ring trajectories for n = 0.8 are plotted in Fig. 11.In each panel, there are two additional curves plotted.The black dashed line corresponds to the power-law case at matched n and Re G and the red dotted line is for a Newtonian vortex ring at the same Re 0 as the curve with the lowest Cu plotted.As expected, the trajectories for the Carreau fluid are intermediate to these of the two limiting fluids, tending to the power-law trajectory as Cu increases and to the corresponding (matched Re 0 ) Newtonian curve as Cu decreases.
The radial vortex growth, as measured by the radial position of the vortex center, and the total circulation decay follow trends that at each instant of time are in accordance to how its viscosity varies with shearrate.At a sufficiently high Cu, i.e. when the initial space-averaged shearrate is such that the Carreau viscosity lies in the power-law region, the vortex growth and the decay of the total circulation are similar to those of a power-law fluid with the same n.As time advances, the velocity and shear-rate decrease, leading to an increase of the viscosity, which enters the Newtonian plateau.Once this happens, the vortex grows and the total circulation decays as for a Newtonian fluid.This behavior can be clearly seen in Fig. 12, where the time evolutions of r ̃c and Γtot / Γtot,max are plotted for n = 0.4 and n = 0.8 at different Cu and fixed Re G = 100.Taking the curve Cu = 1000 for n = 0.8 as an example (Figs.12b and d), the power-law-like behavior is observed for t ̃≲ 100, whereas beyond that time the evolution is Newtonian-like.All parts of the plot include, as solid straight lines, the asymptotic time variations discussed in the previous sections.Naturally, the curves with low values of Cu only exhibit Newtonian-like behavior.To support the link between the dynamic evolution of viscosity and the variation of r ̃c and Γtot / Γtot,max , the insets in Figs.12a and b show the variation with time of the corresponding space-averaged viscosity (η), calculated as η = ∫ B ηω ϕ dS/ ∫ B ω ϕ dS, where B corresponds to the group of cells (faces in an azimuthal plane) outside the cylinder that verify ω ϕ > ω ϕ,max /100.Thus, η represents the vorticity-weighted average viscosity inside the vortex ring region, which is delimited by the isocontour having 1/100 of the maximum vorticity in the core.As expected, the changes of slope in the main plots of Fig. 12 occur approximately at the same time when η transits to a constant value, i.e. when the spaceaveraged viscosity enters the low shear Newtonian plateau.
The results for the axial reach of the dye front are plotted in Fig. 13.The dashed lines represent the value attained by the corresponding power-law fluid (at matched n and Re G ) and they define the low Re 0 (high Cu) asymptote.The black solid lines represent the data for Newtonian vortex rings , plotted in the inset of Fig. 10b (note that Re 0 = Re G for Newtonian fluids), and they define the high Re 0 (low Cu) asymptote.Again, the behavior of z̃d ,f for Carreau fluids is effectively bounded between power-law (high Cu) and Newtonian (low Cu) fluids.Although not tested, we expect the two asymptotic limits to hold for larger Re G , as long as the simulations remain laminar and free of instabilities.
Based on these results, we can conclude that, for Carreau fluids, the relevant form of the Reynolds number is Re G at high Cu and Re 0 at low Cu.For intermediate values of Cu, none of the forms is adequate to uniquely characterize the rich behavior under all conditions, for a single value of n.Instead, Re G is more adequate to define the early instants after the piston stops (power-law-like behavior), whereas Re 0 is better at characterizing and help understand the end of life of the vortex rings (Newtonian-like behavior).These conclusions are in line with the experimental study of Bentata et al. [38] discussed next.

Comparison with experiments
The purpose of this section is to reproduce numerically some of the trajectories reported in the experimental investigation of Bentata et al. [38].The numerical setup was adjusted to reproduce as closely as possible the experimental conditions, including the geometric dimensions of the reservoir (here assumed cylindrical) and fluid rheology.The stroke ratio is 3 and the simulations were performed for their fluid X2, a xanthan gum solution displaying a Carreau-Yasuda shear viscosity function, whose parameters can be found in Bentata et al. [38].The In each panel, the two dashed lines represent the power-law fluid limit at matched (n, Re G ) and the black solid line represents the Newtonian limit (data points were omitted and only the connecting line is shown; see the data points in the inset of Fig. 10b).The red and blue arrows indicate the points with Cu = 10 for each of the two curves (n = 0.4 and n = 0.8), which is approximately the Carreau number marking the transition between the power-law behavior and the Newtonian limit.Note that the Carreau number increases as Re 0 decreases (see the text).The data displayed in this plot was obtained in geometry 2. Carreau model is a particular case of the more general Carreau-Yasuda model [44], when its additional coefficient in the exponent is set to 2. The Carreau number ranges between 13 and 61 for Re G varying between 30 and 200 (our definitions of Cu and Re G match those used by Bentata et al. [38]).
The numerical results are plotted in Fig. 14 along with the corresponding experimental data [38].The time evolution of the axial coordinate of the vortex center is plotted in Fig. 14a, where the symbols represent the experimental data and the dashed lines represent the numerical results obtained with the vortex center definition presented in Section 2.6.There is qualitative agreement between both sets, but the comparison can be further improved by using the same definition for the vortex center as in Bentata et al. [38], i.e. if z c in the simulations is taken as the point on the axis where the velocity magnitude reaches its peak and r c as the radial coordinate of the stagnation points inside the vortex bubble, in a reference frame attached to the vortex ring.The evolution of z̃c according to this definition is plotted as solid lines in Fig. 14a.The trajectory described by the vortex rings is plotted in Fig. 14b, where the symbols, dashed lines and solid lines have the same meaning as in Fig. 14a.Again, there is qualitative agreement between simulations and experiments, which improves with the use of matched definitions for the vortex center.
It is worth noting that the local velocity field used to locate the stagnation points inside the vortex bubble needs the determination of the translational velocity of the vortex rings, which is obtained by simple differentiation of the data in Fig. 14a.Therefore, the small deviations between experiments and simulations observed in Fig. 14a propagate and become amplified in Fig. 14b.Moreover, it was observed that the definition used by Bentata et al. [38] to locate the radial coordinate of the vortex center is highly sensitive to the translational velocity of the vortex center, which is itself a quantity that vanishes over time.For example, small variations of the translational velocity were seen to be sufficient to ensure a monotonic variation of the solid line for Re G = 30 in Fig. 14b, which is the behavior that would in fact make sense.Still, there are certainly other reasons that contribute to the deviation observed, as for example differences in the characterization of the fluids (xanthan gum solutions are not absolutely inelastic), differences in the final position of the piston (at the nozzle exit in the experiments and far from the nozzle exit in the simulations), in the nozzle geometry (finite wall thickness in the experiments and zero wall thickness in the simulations), not to mention experimental uncertainty.Therefore, the level of agreement between simulations and experiments is, in general, good for the cases simulated.The results presented in this section also serve to draw attention to the differences that may result from the simple use of unmatched criteria to characterize vortex rings.

PTT fluids
The simulations with PTT fluids were performed fixing β = 0.01, ε = 0.05 and L/D = 2 in geometry 2. The two variable parameters were the Deborah number, De = {0, 0.2, 0.5, 1, 2, 5, 10, 20}, and the Reynolds number, Re 0 = {10, 25, 50}.These conditions do not reproduce any particular fluid/flow reported in the experiments with viscoelastic fluids [39][40][41], but they do offer a reasonable approximation to that ensemble in the laminar regime.While it is worth exploring higher values of L/D, Re 0 and De in future studies, the ranges considered in this work are wide enough to meet our intents and straddle the experimental range.We remember that a distinguishing feature observed experimentally in laminar viscoelastic vortex rings was the flow reversal [39][40][41], accompanied by the unwinding of the vortices [40,41].The first and major purpose of the simulations of viscoelastic vortex rings is to reproduce this behavior.

Dye patterns
The time evolution of the dye front (on the axis) is plotted in Fig. 15 for different combinations of Re 0 and De.For Re 0 = 10 and De up to about 0.5, the dye front evolves forward continuously even if the axial reach decreases with De.For De above 0.5, the evolution of z̃d exhibits oscillations over time, including a backward motion of the front.For the other higher values of Re 0 , the curves still exhibit oscillations beyond a critical De, but they become weaker as inertia increases (note that El = De/Re).Thus, inertia seems to have a stabilizing effect and promotes (or reduces the departure from) a Newtonian-like behavior.The backward motion of the dye front is an elastic effect, driven by normal stresses that develop and relax near the axis.Due to the fading memory of viscoelastic fluids, such polymeric stresses develop and relax over time scales of the order of the relaxation time.In addition, the magnitude of the polymeric stresses is also proportional to the relaxation time.If this characteristic fluid time scale is shorter than or of the order of the other flow time scales (diffusion and advection), the polymer stresses will impact the vortex dynamics early, but if the relaxation time is too large the impact may be smaller because the vortex may dissipate or become too weak  before the polymer stresses have a chance to interfere.As we see in Fig. 15, for Deborah numbers of order 1 these polymer stresses counteract the self-induced forward motion of the vortex rings, which is also affected by viscous effects, but as De exceeds 1 the effect of polymer stresses seems to change in nature.Simultaneously, a higher De also has implications on the shear viscosity with the fluid becoming less viscous due to its shear-thinning at lower strain rates (both the zero-shear polymer and solvent viscosities are kept constant).Therefore, the occurrence or not of a backward motion of the dye front depends not only on the balance between inertial, viscous and elastic forces, but also on the relation between relevant time scales.It is important to note that, in some cases, there is a local backward flow (e.g. in the tail of the vortex ring) even in the absence of a visible backward motion of the dye front.Also note that the backward motion of the dye front of viscoelastic fluids is unrelated to the backward diffusion of the vortex center near the end of their trajectories (Fig. 3) reported earlier for power-law and Newtonian fluids.
The dye patterns developed over time are illustrated in Fig. 16 for two vortex rings at De = 20, for Re 0 = 10 and 25.The dye patterns for these and other combinations of Re 0 and De were also recorded in Video 4 (Supplementary material) during the whole lifecycle of the vortex rings.Both the figure and the video show the backward motion of the dye front for the cases detected earlier in Fig. 15.Moreover, for De = 10 and De = 20, it is also possible to observe another interesting behavior, which was absent for the inelastic fluids analyzed in the previous sections.Indeed, even though the vortex rings expand radially once they are released from the cylinder exit (t ̃= 2), as observed for inelastic fluids, they start to unwind a few instants later, while at the same time they contract and decelerate (eventually inverting their motion).The unwinding and contraction of the vortex rings can be better seen in Videos 5-7 (Supplementary material) for De = 10 and 20, where the computational domain was seeded with massless particles that move with the flow.The initial anti-clockwise motion of the particles inside the vortex ring (winding) progressively slows down and, eventually, transforms into a clockwise motion (unwinding).The trajectories described by the particles also evidence an extensional flow along the axis, which can eventually break the dye filament due to diverging streamlines.
The dye patterns obtained corroborate the flow reversal reported in experiments with viscoelastic shear-thinning fluids [39][40][41], as well as the unwinding of the vortex core [40,41].However, it is worth noting that whereas Palacios-Morales et al. [39] based their observation on velocity measurements, Hegde et al. [40] and Shashank and Sreenivas [41] reported the phenomena based on the visualization of dye patterns.The two methods have different sensitivities in the detection of flow reversal, as a local and possibly temporary reverse flow (detectable by probing the local velocity field) does not necessarily pull the dye front backwards if diffusion (either physical or numerical) is non-negligible.On the other hand, the backward motion of the dye front necessarily points out to flow reversal.
The vortex unwinding and flow reversal are manifestations of a popular non-Newtonian effect known as elastic recoil [44], which can be found in different types of viscoelastic fluid flows.For example, the main vortex that develops in the lid-driven cavity flow of a viscoelastic fluid can reverse its rotation after the moving wall stops [61] and the Poiseuille flow of a viscoelastic fluid can also reverse its motion when the external pressure gradient is relieved [44].The relaxation of stresses (partial structural recovery of polymer molecules) is responsible for such phenomena and does not require a variable shear viscosity (e.g. a Boger fluid can also show elastic recoil).

Vorticity
The vorticity decay in power-law fluids was shown to be faster than in a Newtonian fluid (cf.Fig. 7), but the vorticity patterns/isocontours were qualitatively similar among the two classes of fluids (Video 3, Supplementary material).The transport equation of vorticity in viscoelastic fluids has an additional source term proportional to the polymeric stress tensor and this introduces significant changes in the vorticity patterns.Video 8 (Supplementary material) shows the time evolution of the vorticity contours for some of the PTT fluid cases simulated.After the piston stops, the vorticity field develops multiple peaks, as well as neighboring regions of vorticity with opposite sign (this renders the vortex tracking algorithm ineffective when it is based on the identification of the point of maximum vorticity).This was an expected result considering the backward motion and vortex unwinding observed in the dye patterns.However, after a given period of time, which seems to increase with De and decrease with Re 0 , the vorticity patterns rearrange to form Newtonian-like patterns, i.e. quasi-elliptic isocontours of vorticity around a single well-defined peak.This behavioral transition reflects the vanishing of elastic effects, as shear-thinning is unlikely the cause behind the disperse patterns of vorticity, otherwise that would be also observed in the simulations with power-law and Carreau fluids.
Although the conditions imposed in the numerical simulations do not match exactly those used by Palacios-Morales et al. [39], it is worth noting the differences in the vorticity patterns between both cases.The maps of vorticity in the experiments show the formation of a well-defined Newtonian-like vortex ring that is quickly annihilated by a secondary vortex of opposite vorticity that grows in front of it and that originates the flow reversal [39].The presence of the secondary vortex of opposite vorticity was not evident in the simulations, since the vorticity is significantly more dispersed in space than was seen in the experiments.However, the numerical maps of vorticity clearly show the existence of neighboring regions with opposite rotation (c.f.Video 8).

Axial reach
The analysis of the axial reach of the dye front reveals the interplay between elasticity and shear-thinning effects.Such analysis is carried out in Fig. 17, where z̃d ,f is plotted against Deeach curve represents a different Re 0 .For De ≲ 1, z̃d ,f decreases with De and this is essentially an elastic effect.In fact, in the low-De range the viscosity is almost constant and independent from De (Fig. 2), such that the space-averaged Reynolds number of the flow can be effectively approached by Re 0 .In terms of viscous effects, this is analogous to the low Cu limit previously analyzed for Carreau fluids, corresponding to the Newtonian plateau of viscosity.Therefore, when fixing Re 0 while varying De, but keeping its value below ~1, the level of inertia can be effectively considered constant and the increased elastic effects felt by the vortex rings when De is incremented leads to a lower axial reach.On the other hand, for De ≳ 1 the axial reach starts to rise with De and this is mainly due to the predominant effect of shear-thinning over elasticity.Elastic effects continue rising with De and keep resisting the vortex motion, but the effective Reynolds number of the flow is no longer Re 0 .This range of De already falls in the power-law-like region of viscosity decay of the PTT fluid (Fig. 2) and, consequently, the space-averaged Reynolds number of the flow is higher than Re 0 .This decrease of viscous effects, not reflected in the constant value of Re 0 (in contrast with the analysis for power law fluids, that used a generalized Reynolds number accounting for such effects of shear-thinning), is behind the increase of z̃d ,f with De.Note that besides the effect on the average Reynolds number, the De-dependent viscosity also influences the vortex rings via spatial gradients of viscosity.Both effects co-exist and contribute to the results observed, but the effect of De on the average value of Re, i.e. through the viscosity, is likely more important to the trends shown in Fig. 17.
In this analysis, an additional characteristic of viscoelastic fluids (already introduced in Section 4.3.1)needs to be considered.In these fluids, the relaxation time governs the growth of the polymeric stresses over time, which in turn influences the development of the velocity profile inside the cylinder.The Deborah number in the simulations was varied by changing the relaxation time, such that a higher De also means a higher relaxation time and, consequently, a slower elastic response.This effect is visible in Fig. 15, where the first oscillations appear later as De increases.Considering that the piston is only active up to t ̃≈ L/Dthis is the time left for the stresses and velocity profiles to develop inside the cylinder -, then one can conclude that a higher De may not necessarily increase the amount of elastic energy fed to the flow if λ >> L/U.A new dimensionless number can be formed to assess this effect, De L = λU/ L, which is related to De by De/De L = L/D.This relation should be explored in future works through the variation of the stroke ratio.

Discussion and conclusions
This work presents a numerical study of viscous and viscoelastic vortex rings generated in non-Newtonian fluids, with special emphasis on the post-formation phase.The fluid rheology and flow conditions were selected in accordance with the few experimental studies available.Power-law, Carreau and PTT rheological models, encompassing both inelastic and viscoelastic shear-thinning fluids, were tested and compared against the Newtonian behavior.The vortex rings were formed under laminar flow conditions (10 1 ≤ Re ≤ 10 3 ) assuming flow axisymmetry, as found in the experiments, and the stroke ratio was kept low (L/D ≤ 4) in order to avoid a long trailing wake of vorticity.The range of conditions addressed in this work, along with a summary of the main findings are presented in Table 3 and briefly discussed below.
The variable viscosity of power-law fluids introduces three major differences compared to Newtonian fluids: (i) a smaller amount of circulation is produced due to differences in the velocity profile at the cylinder exit; (ii) the occurrence of spatial gradients of viscosity in the flow; (iii) an enhanced reduction of the local space-averaged Reynolds number over time.These effects contributed together and up to different extents to the features observed in power-law vortex rings.At sufficiently high Re G , shear-thinning has shown little influence in the earlystage of the post-formation phase, which is advection-dominated and where the diffusive term accounting for viscosity effects in the vorticity transport equation has a relatively small importance.The radial growth of the vortex rings in this phase is nearly independent of n and there are only little differences in the trajectories between fluids.As time advances, the vorticity contained within the vortex bubble diffuses and dissipates, the vortex rings expand while the surrounding fluid is entrained into the vortex bubble and the vortices slowdown, leading to a progressive reduction of the local space-averaged Reynolds number.The flow enters a diffusion-dominated phase and shear-thinning effects become evident.The reduction of the local space-averaged Reynolds Fig. 17.Axial reach of the dye front for PTT fluids.The lines are only a guide to the eye and the data displayed in this plot was obtained in geometry 2.

Table 3
Summary of the range of fluid and flow conditions tested and of main vortex ring features observed.Although not specifically listed below, Newtonian fluid flows were also investigated and used for comparison, as they are limiting cases of all fluid types.number is boosted in shear-thinning fluids, since it is not only due to a decrease of the velocity (common to Newtonian fluids), but also to an increase of the viscosity, as local shear-rates vanish.This is the cause, for example, for the faster decay of velocity and vorticity (circulation), as well as the reduced distance travelled by power-law vortex rings and the faster radial expansion observed in this phase.Although no real fluid exhibits a truly power-law rheology over the whole range of shear-rates due to the unphysical behavior near the extremes, this rheological model is still useful as the limit of the more realistic Carreau inelastic fluid model.Besides, for intermediate shear-rates the power-law model should provide an acceptable approximation to real fluids.
The Carreau model offers a better approach to real fluids through the addition of lower and upper bounds to the variable viscosity, although the lower bound (high shear-rate plateau) was not considered in this work.Its rheology is therefore bounded between that of a power-law fluid at high shear-rates and that of a Newtonian fluid at low shearrates.In practice, the viscosity boundedness is reflected on the vortex ring properties.A Newtonian-like behavior is reached in the limit Cu→0, whereas a power-law-like behavior is recovered in the limit Cu→+∞.Interestingly, for intermediate but sufficiently high values of Cu, a vortex ring can start exhibiting a transient power-law-like evolution and evolve after some time to a Newtonian-like behavior, when the local space-averaged viscosity enters the Newtonian plateau.This was observed, for example, in the time evolutions of the radial expansion rate and circulation decay.The Carreau number and its relation to a characteristic Carreau number of order 1, corresponding to the shearrate at the transition from the Newtonian plateau of viscosity to the power-law region of viscosity decay, is therefore the tuning parameter between the two extremes, in agreement with the experimental study by Bentata et al. [38].This should be taken into account in the selection of the most adequate form of the Reynolds number to describe the flow, either the form based on the Newtonian plateau (Re 0 ) or the one based on the power-law region of viscosity decay (Re G ).
In a viscoelastic fluid, vortex rings experiment inertial, viscous and elastic forces.Upon cessation of the piston motion, inertia convects the vortex rings away from the cylinder exit.The elastic force resulting from tensile stresses builds up over time and opposes this motion, whereas viscous forces damp the vortex motion irrespective of its direction.The balance between these three forces, measured by De and Re 0 , leads to different dye patterns, including a backward motion of the dye front and the unwinding of the vortex coretwo visible effects of an elastic recoil that were observed in experiments [39][40][41] and that were absent in the simulations with inelastic fluids.The use of Re 0 for viscoelastic shear-thinning fluids does not account for viscosity changes while varying De (this is inherent to the constitutive equation and not to the particular flow).Therefore, increasing De while keeping Re 0 fixed contributes to increase the relevance of inertia through the reduction of viscosity, which, in general, is not proportional to the increase of elasticity, since the buildup and relaxation of elastic forces also depend directly on De.This explains, for example, the increase of the axial reach with De, at large De, for fixed Re 0 , notwithstanding the opposing elastic force, which was seen to effectively reduce the axial reach for De ≲ 1 (when the shear viscosity still lies in the Newtonian plateau and elastic forces build up quickly).The vorticity maps presented several peaks (with both positive and negative signs) during the post-formation phase, but the vorticity patterns rearrange after some time (once the elastic effects vanish) and give rise to Newtonian-like (quasi-elliptic) isocontours with a single and well-defined peak.Some similarities can be drawn between PTT and Carreau fluids, as both ultimately transit to a Newtonian-like behavior when t ̃→+∞ due to their constant shear viscosity at low shear rates and vanishing (local) elasticity for the viscoelastic case.
To conclude, this work has shown that a non-Newtonian rheology can influence, in many ways, the behavior of vortex rings generated in a piston-cylinder setup and that the finite-volume method is a viable option to capture and study those changes.The results of this work contribute to a better understanding of vortex rings in complex fluids and can find application in processes where impulsive jets occur in non-Newtonian fluids.This numerical study can be extended in numerous ways, as the investigation of non-Newtonian vortex rings is still incipient.Viscoelastic vortex rings displayed clearly a more disruptive behavior in relation to the Newtonian case, which should justify a more detailed study for such kind of fluids, exploiting other constitutive equations, widening the space of design parameters (Reynolds and Deborah numbers, solvent viscosity ratio, stroke ratio, etc.) and assessing possible loss of axisymmetry, especially when the flow instabilities appear due to elasticity.A judicious choice of the constitutive equations should allow a clear splitting between elastic and shear-thinning effects, which is not always easy in experimental studies [39][40][41].Another interesting research direction is to assess the effect of a non-Newtonian rheology on the existence of a formation number, as observed for Newtonian fluids [1,22] (it was already confirmed experimentally by Palacios-Morales and Zenit [37] for inelastic shear-thinning fluids).This forwards the suggestion left by Palacios-Morales et al. [39] related with the conservation of vorticity in viscoelastic fluids.A more challenging endeavor should be the numerical simulation of the turbulent viscoelastic vortex rings produced in the experimental study by Hegde et al. [40] and Shashank and Sreenivas [41].An equally challenging and not less important task would be the development of analytical solutions for non-Newtonian vortex rings, similar to those discussed in Section 3 for Newtonian fluids.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

η 0 .
The elastic effects can be quantified through the Deborah number, here defined as De = λU D , or alternatively by the elasticity number, El = De/Re 0 .The remaining dimensionless numbers for PTT fluids are the solvent viscosity ratio, β = η s η 0 , and the extensibility parameter (ε).

F
.Pimenta et al.

Fig. 3 .
Fig. 3. Trajectories described by the vortex ring center for power-law fluids with (a) n = 0.4; (b) n = 0.6; (c) n = 0.8; (d) n = 1 (Newtonian case).The circles plotted in (b) aim to exemplify variable z̃c ,1 (maximum axial reach of the vortex center before any significant change of the trajectory) for n = 0.6, which is then plotted in Fig. 10 for all values of n.

F
.Pimenta et al.

Fig. 5 .
Fig. 5. Evolution with time of the translational velocity of the vortex center for (a) n = 0.6 and (b) n = 1.Note that the piston only stops its motion at t ̃≈ 4, which marks the beginning of the post-formation phase.The power-law exponent of − 1.5 is that of the theory of Eq. (7).

F
.Pimenta et al.

Fig. 6 .
Fig. 6.Dye (lower-half of each panel) and vorticity (upper-half of each panel) patterns at three different instants of time for n = 0.6 and n = 1, at Re G = 100.The vorticity is normalized by its instantaneous maximum value outside the cylinder (ω φ,max ).

Fig. 7 .
Fig. 7. Decay of the total circulation with time for different levels of inertia and two power-law exponents: (a) n = 0.6 and (b) n = 1.The data are normalized by the maximum total circulation registered during the whole life cycle of the vortex rings ( Γtot,max ), which is plotted in Fig. 9.This maximum value usually occurs just after the vortex leaves the cylinder lip.

F
.Pimenta et al.

Fig. 8 .
Fig. 8.Total circulation plotted as a function of the radial coordinate of the vortex ring center for (a) n = 0.4; (b) n = 0.6; (c) n = 0.8; (d) n = 1.The data is normalized by the maximum total circulation ( Γtot,max ), which is plotted in Fig. 9.

Fig. 9 .
Fig. 9. Maximum total circulation ( Γtot,max ) as a function of Re G for different power-law exponents.The lines are only a guide to the eye.

4. 2 .
Carreau fluids 4.2.1.Generic behavior The simulations with Carreau fluids were performed for two values of n (0.4 and 0.8), two values of Re G (50 and 100) and multiple Cu values, ranging from Cu = 0.005 up to Cu = 1000.While fixing n and Re G , different values of Cu were simulated by varying Λ.The following relation between Cu and Re 0 can be deduced, Re 0 = Re G 3n+1 4n Cu n− 1 , i.e.

Fig. 10 .
Fig. 10.Axial reach of vortices: (a) first peak of axial reach for the vortex center (z̃c ,1 ) plotted as a function of Re G and (b) corrected axial reach of the dye front plotted as a function of Γtot,max Re G for different power-law exponents.The inset in (b) plots the axial reach of the dye front without being corrected by its position at t = L/D (z̃d ,0 ) as a function of Re G.The lines are only a guide to the eye and the results represented in (b) were obtained with geometry 2.

Fig. 11 .
Fig. 11.Vortex ring trajectories for Carreau fluids: (a) n = 0.8, Re G = 50; (b) n = 0.8, Re G = 100.In each panel, the black dashed line represents the corresponding power-law fluid case at matched (n, Re G ) and the red dotted line represents the Newtonian case at matched Re 0 for the curve with the lowest value of Cu.

Fig. 12 .
Fig. 12.Time evolution of the radial coordinate of the vortex center (a and b) and of the normalized total circulation (c and d) for Carreau fluids at Re G = 100: (a) and (c) n = 0.4; (b) and (d) n = 0.8.In each panel, the black dashed line represents the corresponding power-law fluid case at matched (n, Re G ) and the red dotted line represents the Newtonian case at matched Re 0 for the curve with the lowest value of Cu.The insets in (a) and (b) show the corresponding evolution of the spaceaveraged viscosity inside the vortex (see the main text).

Fig. 13 .
Fig. 13.Axial reach of the dye front (z̃d ,f ) as a function of Re 0 (Cu) for two Carreau fluids (n = 0.4 and n = 0.8), at Re G = 50 (lower panel) and Re G = 100 (upper panel).In each panel, the two dashed lines represent the power-law fluid limit at matched (n, Re G ) and the black solid line represents the Newtonian limit (data points were omitted and only the connecting line is shown; see the data points in the inset of Fig.10b).The red and blue arrows indicate the points with Cu = 10 for each of the two curves (n = 0.4 and n = 0.8), which is approximately the Carreau number marking the transition between the power-law behavior and the Newtonian limit.Note that the Carreau number increases as Re 0 decreases (see the text).The data displayed in this plot was obtained in geometry 2.

Fig. 14 .
Fig. 14.Comparison between our numerical results and the experimental results of Bentata et al. [38] with a xanthan gum solution (fluid X2) displaying a Carreau-Yasuda rheology.Panel (a) plots the time evolution of the axial coordinate of the vortex center and panel (b) represents the trajectory described by the vortex center.The symbols represent the experimental data, the dashed lines represent the numerical results obtained with the vortex center definition presented in Section 2.6 and the solid lines represent the same numerical results, but with the vortex center as defined in Bentata et al. [38].

Fig. 15 .
Fig. 15.Time evolution of the dye front for PTT fluids at (a) Re 0 = 10, (b) Re 0 = 25 and (c) Re 0 = 50.The data displayed in this plot was obtained in geometry 2.

F
.Pimenta et al.

Fig. 16 .
Fig. 16.Snapshots of the dye patterns obtained in geometry 2 at four instants of time for two PTT vortex rings (De = 20) at Re 0 = 10 and Re 0 = 25.

Table 2
Solvent and polymeric stresses for the different constitutive equations used in this work.