Asymmetric flows of complex fluids past confined cylinders: A comprehensive numerical study with experimental validation

Flow around a circular cylinder is considered to be one of the benchmark problems in non-Newtonian fluid dynamics. This fundamental geometry is representative of a wide range of systems involving flows of viscoelastic fluids^1,2 and has long been a popular choice for testing constitutive models and numerical schemes^3–7 against the abundant experimental data.^8–14 With the advent of microfluidics since the early 2000s,^15,16 the experimental study of the viscoelastic flow around a cylinder has entered an interesting regime in which inertia is almost negligible (i.e., Reynolds numbers Re ≪ 1), while elasticity can be simultaneously high (Weissenberg numbers Wi ≫ 1).^17–28 Such highly elastic flows combining the streamline curvature are well known for being prone to interesting and important elastic flow instabilities, and even “elastic turbulence” (a turbulent-like state that can arise in viscoelastic fluids even in the absence of inertia).^19,26–33 However, such flows remain as a great challenge to simulate numerically due in large part to the high Weissenberg number problem (HWNP) along with other numerical instabilities that cause solutions to diverge (see, e.g., Ref. 7).

A two-dimensional (2D) approximation of the geometry of a cylinder confined inside a channel can be described succinctly by the “blockage ratio” B_R = 2R/W, where R is the radius of the cylinder and W is the width of the channel [see Fig. 1(a)]. The flow in such a geometry presents a complex mix of kinematics. For an average flow velocity U inside the channel, the nominal wall shear rate of the plane Poiseuille flow upstream of the cylinder is $γ̇w=6U/W$⁠.³⁴ Approaching the cylinder, fluid elements undergo compression at the upstream stagnation point of the cylinder, where streamlines are caused to diverge into the gaps between the cylinder and the channel side walls. In these gaps (to an extent that depends on B_R), shear rates can be significantly higher than those in the upstream sections of the channel with a nominal value $γ̇w,gap=12U/W(1−BR)2$⁠.³⁵ If B_R is high, there may also be an appreciable transient extensional component in the flow as fluid elements approaching from the upstream are squeezed and subsequently expanded as they pass through the gaps. Finally, at the downstream stagnation point of the cylinder, fluid elements are subjected to an extensional deformation with a nominal elongational rate $ϵ̇=U/R$ and also a high residence time.³ 

The dynamics of viscoelastic flows around microfluidic cylinders are strongly affected by the blockage ratio. For high values of B_R ≳ 0.5, the gaps between the cylinder and the channel walls dominate and the dynamics resemble those of viscoelastic flows through microfabricated contractions, presenting instabilities and flow recirculations upstream of the cylinder.^{18,21–24,36,37} If B_R is reduced, thus reducing squeezing and shearing around the sides of the cylinder, a greater prominence is placed on the influence of the stagnation points. In this case, for $Wi=λϵ̇≳1$ (where λ is a characteristic relaxation time of the fluid), greatly elongated wakes are observed downstream of the cylinder due to the extensional stresses generated at the stagnation point.^3,25,38

Of particular note and relevant to the present work are recent experiments on flows around slender cylinders with low blockage ratios (B_R = 0.1).^28,35,38,39 For certain fluids tested in such devices, an intriguing new viscoelastic flow instability is observed when the Weissenberg number exceeds a fluid-dependent critical value Wi_c ∼ O(10–100). The instability, which can be quite steady in time over a range of Wi > Wi_c, is characterized by the preferential passage of fluid around one side of the cylinder (i.e., a lateral asymmetry of the flow field) and appears to be a pitchfork bifurcation with a random selection of the preferred side of passage.^35,39 In some cases, the degree of asymmetry is so extreme that nearly all the fluid is observed to pass on only one side of the cylinder, with a region of nearly stagnant fluid developing on the opposite side.^35,39 It is an experimental observation to date that only fluids exhibiting both shear-thinning and extensional elastic properties seem to be capable of developing this strongly asymmetric state.^35,38–40 The limited experiments conducted to date indicate that this flow state probably arises due to a combination of high extensional stress in the downstream wake and shear-thinning around the sides of the cylinder. It has been proposed that the flow imbalance may be initiated by some small random fluctuation of the high-stress wake. This could cause a small imbalance in the shear rate on either side of the cylinder, which becomes compounded if the fluid is shear-thinning.^35,39 However, gaining a clear insight into the problem by an experimental approach is complicated due to the difficulty in performing controlled rheological and/or geometric parameter studies allowing for independent variation of shear-thinning and extensional properties (or $γ̇w,gap$ and Wi).

In this work, we employ a recently developed finite element method (FEM) scheme that circumvents the HWNP⁷ in order to numerically examine flows of various viscoelastic constitutive models around cylinders, reaching previously inaccessible high Weissenberg numbers. Simulations are performed with a shear-thinning but inelastic fluid model (the Carreau–Yasuda, C–Y, model), an elastic but non-shear-thinning model (a finitely extensible non-linear elastic dumbbell model with Chilcott–Rallison closure, FENE-CR model), and a model that is both elastic and shear-thinning (the linear Phan–Thien–Tanner, l-PTT, model). The parameters of the fluid models are tuned so that their rheometric responses coincide with a series of experimental test solutions that are used to validate the simulations via complementary flow velocimetry measurements. Our aims are (1) to demonstrate the ability of the new FEM scheme to accurately predict viscoelastic flows with complex mixed kinematics up to high Wi, (2) to capture numerically the lateral flow asymmetry observed in experiments, (3) to validate the numerical method by comparison with complementary flow visualization experiments, (4) to verify and understand the dependence of the flow asymmetry on the nature of the fluid rheological properties, and (5) to numerically examine how the geometric parameter B_R affects the flow instability, which has not yet been studied (and would be difficult to study) experimentally.

In the numerical simulations, we consider the 2D creeping flow of a non-Newtonian fluid in a planar channel that features a circular cylinder located at the coordinate origin mid-way between the walls [see Fig. 1(a)]. The fluid is incompressible with constant density ρ, solvent viscosity η_s, zero shear rate viscosity η₀, and characteristic relaxation time λ. The width of the channel is denoted W and the radius of the cylinder R. The length of the channel is L = 250R, and in all simulations, the flow domain spans −125 ≤ x/R ≤ 125. The fluid is driven into the channel by the action of a pressure gradient generating a volumetric flow rate Q = UW per unit depth, where U is the average flow velocity of the fluid far from the cylinder. Figure 1(a) shows the flow geometry where the coordinates x and y correspond to the local streamwise velocity and its normal direction, respectively. We scale all lengths with the cylinder radius R, velocities with the average flow velocity U, and times with the characteristic flow time R/U. In addition, both the pressure and stress components are scaled with a viscous scale, η₀U/R. All nondimensional quantities will be indicated by the superscript “^*,” i.e., lengths (x^*, y^*) = (x/R, y/R) and velocities u^* = (u^*, $v$^*) = (u/U, $v$/U). The dimensionless groups that arise are the Weissenberg number, Wi = λU/R, the Newtonian solvent-to-total viscosity ratio β = η_s/η₀, and the blockage ratio B_R = 2R/W. In most of the simulations, inertia is neglected (i.e., the Reynolds number Re ≡ 0); however, we also confirm that this simplifying assumption is valid by including inertia in specific test cases.

The non-Newtonian flow is described by the incompressible and isothermal Cauchy equations coupled with a constitutive equation, which accounts for the contribution of the non-Newtonian stresses. Neglecting inertia, the forms of the continuity, momentum, and constitutive equations are expressed, respectively, as

where u is the velocity vector, P is the thermodynamic pressure, I is the identity tensor, and τ is the non-Newtonian contribution to the total stress tensor. The deformation rate tensor, $γ̇$⁠, is defined as

where the superscript “T” denotes the transpose operator.

The usual no-slip and no-penetration boundary conditions (i.e., u = 0) are imposed on the cylinder surface and channel walls [see Fig. 1(b)]. At the inflow (x = −125R), we impose fully developed velocity and stress fields, i.e., one-dimensional (1D) equations for the velocities and the components of the stress tensor that are solved together with the 2D equations for the rest of the domain. More specifically, in each time-step, we solve the governing equations for 1D flows in the channel, under the constraint of a flow rate that is gradually increased from zero to UW,

The 1D flow profiles thus obtained (i.e., velocity and logarithm of the conformation tensor) are then imposed as Dirichlet conditions at the inflow boundary. Note that, in the case of direct steady state simulations, the term (1 − e^−tU/R) in Eq. (5) is dropped.

In order to eliminate any additional numerical error that could arise due to the truncation of the domain, the open boundary condition (OBC)⁴¹ has been applied along the outflow boundary at x = 125R. According to the OBC, the fluid velocities and stresses are not imposed at the outflow boundary but are calculated from the weak form of the equations for both velocity unknowns (extrapolated from the bulk).

The form of the operator g [given as a function of τ and $γ̇$ in Eq. (3)] depends on the choice of the constitutive model. Three constitutive models are selected for the numerical simulations, which represent three broad classes of fluids and allow us to understand the effects of different fluid rheologies: (1) shear-thinning but inelastic, (2) elastic but constant viscosity, and (3) both shear-thinning and elastic. The three models are discussed in detail in Subsection II C.

The rheological response of an inelastic, but shear-thinning fluid is modeled using the Carreau–Yasuda (C–Y) model.⁴² Fluids that obey this constitutive equation are classified as generalized Newtonian fluids (GNFs). The C–Y model can be used to describe shear-thinning effects in the steady flow curve but does not include memory effects and thus does not develop normal stresses related to microstructural deformation of a fluid. The C–Y constitutive equation is presented as

where $γ̇c$ is the characteristic shear rate for the onset of shear-thinning, n is the “power-law exponent” in the shear-thinning region, and a is a dimensionless fitting parameter that controls the rate of the transition between the constant viscosity and shear-thinning regions of the flow curve. The magnitude of the deformation rate tensor is denoted as $γ̇=0.5γ̇:γ̇$⁠. In our simulations, the adjustable parameters in Eq. (6) were selected by non-linear regression to the experimental flow curve measured with a very weakly elastic solution of xanthan gum in a glycerol:water mixture (fluid termed XG, see Sec. III B); see Fig. 2(a). We obtained η₀ = 8.1 Pa s, η_s = 0.1 Pa s, $γ̇c$ = 0.11 s⁻¹, n = 0.33, and a = 0.33. Figure 2(b) shows the results of capillary thinning measurements made with the experimental test fluids (see Sec. III B for details). In the case of the XG test fluid, the mid-filament diameter D decays almost linearly with time t until pinchoff occurs, as expected for a Newtonian (i.e., inelastic) fluid.^43,44 There is a slight departure to a fast exponentially decaying region immediately before pinchoff, which is the signature of the weak elasticity of the fluid.^45,46 In Fig. 2(c), we plot the steady uniaxial elongational viscosity (η_E) predicted by the C–Y model as a function of the elongation rate $ϵ̇$⁠. Since the GNF model is inelastic, the Trouton ratio in uniaxial elongation is Tr = η_E/η ≡ 3 and, over the range shown, the elongational viscosity thins monotonically with the increasing elongation rate, mirroring the flow curve in steady shear.

We model the rheological response of elastic fluids with constant shear viscosity using the finitely extensible non-linear elastic dumbbell model with the Chilcott–Rallison closure (FENE-CR model).³ The fluids that obey this constitutive equation are classified as Boger fluids. This means that these fluids do not exhibit shear-thinning effects (the shear viscosity η is independent of the shear rate), but they feature strong elastic effects and thus intense extension-rate-hardening. The FENE-CR constitutive equation is given in terms of the conformation tensor C as

Here, the inversed triangle over the conformation tensor denotes the upper convected derivative and f accounts for the finite chain extensibility effect,

where L² is the extensibility parameter. Additionally, the stress tensor is related to the conformation tensor as follows:

where G is the elastic modulus of the viscoelastic fluid.

The parameters of the FENE-CR model are selected to approximate an almost non-shear-thinning experimental test solution composed of 0.2 wt. % of a high molecular weight (M_$w$ = 4 MDa) poly(ethylene oxide) (PEO) dissolved in a glycerol:water mixture (termed PEO4, Sec. III B). The solvent viscosity is η_s = 8 mPa s and the elastic modulus is G = 0.175 Pa, which results in a steady shear flow curve in good agreement with that of the experimental test fluid [see Fig. 2(a)]. The relaxation time is set to match that measured for the fluid used in the experiments (λ = 0.044 s, Sec. III B), determined from the fit to the pronounced exponentially decaying portion of the capillary thinning measurements shown in Fig. 2(b). Note that this exponentially decaying “elastocapillary” thinning region is predicted by FENE-type fluid models.^45,46 For most of the simulations, the extensibility parameter is set to L² = 5000, which is appropriate for a 4 MDa PEO molecule,⁴⁷ although we also investigate the effect of varying this parameter in the model. As we can see in Fig. 2(c), the FENE-CR model predicts strong elastic behavior manifested by significant extension-rate-hardening.

Finally, we use the linear Phan–Thien–Tanner (l-PTT) model⁴⁸ to simulate the response of fluids that exhibit both shear-thinning and elasticity. The linear version of the PTT model is preferable over its exponential counterpart when considering viscoelastic solutions⁴⁹ because of its capacity to predict both shear-thinning and monotonic extension-rate-hardening. The l-PTT constitutive equation is given in terms of the conformation tensor C as

where Y is the linear PTT function,

The magnitude of the PTT parameter $ε$ governs the extension-rate hardening of the fluid: lower $ε$ implies stronger extension-rate hardening. For the range of values of $ε$ that we will consider, shear-thinning is mainly controlled by the solvent-to-total viscosity ratio β: smaller β implies stronger shear-thinning.^7,49 The stress tensor is related to the conformation tensor as follows:

In the l-PTT model, the elastic modulus G = 2.1 Pa, solvent viscosity η_s = 0.015 Pa s, PTT parameter $ε$ = 0.05, and solvent-to-total viscosity ratio β = 0.05. These parameters are selected so that both the flow curve in steady shear [Fig. 2(a)] and the exponentially decaying region of the capillary-thinning curve [Fig. 2(b)] are minimized with respect to an experimental test fluid composed of 0.5 wt. % of an M_$w$ = 8 MDa PEO in water (termed PEO8, Sec. III B). Specifically, the capillary-thinning data of the experimental sample are fitted in the elastocapillary regime by non-linear regression to the relation originally proposed by Anna and McKinley, ^43,46,50,51

where A and B are constants and σ = 0.07 N m⁻¹ is the surface tension of the fluid sample. Equation (13) provides the relaxation time λ = 0.145 s of both the experimental fluid sample and the l-PTT model fluid, and also an estimate for the upper plateau of the steady uniaxial elongational viscosity η_E,∞ = 15.12 Pa s. ⁵¹ This value of η_E,∞ is used to determine the parameter $ε$ = 0.05 such that the extensional response of the l-PTT model reaches a similar plateau value at high extensional rates [Fig. 2(c)]. Clearly, as shown in Fig. 2(c), this l-PTT model predicts quite strong extension-rate hardening in uniaxial elongational flows. In our simulations, we also investigate the effects of varying the degrees of shear-thinning and elasticity by manipulating $ε$ and β.

The Petrov–Galerkin stabilized Finite Element Method for Viscoelastic flows (PEGAFEM-V), which was proposed by Varchanis et al.,⁷ is used to solve the governing equations. The aforementioned finite element (FE) method makes use of linear interpolants for all variables and combines classical finite element stabilization techniques^52–54 with the log-conformation representation of the constitutive equation.⁵⁵ The variational formulation along with a detailed explanation of the FE method is given by Varchanis et al.⁷ 

The mesh is generated by the quasi-elliptic mesh generator introduced by Dimakopoulos and Tsamopoulos.⁵⁶ In all simulations, we use triangular elements. To check the mesh convergence of our numerical solutions, we used three consecutively doubled meshes, whose characteristics are quoted in Table I. Mesh M2 was used in all other simulations. Mesh M1 in a region of the flow domain surrounding the cylinder is represented in Fig. 1(b). Time integration is performed using a fully implicit second order backward finite difference scheme preceded by a quadratic extrapolation step for the prediction of the solution at each new time instant. Consequently, the numerical method features second order accuracy in space and time.

For a given set of flow parameters (Wi, β, etc.), we perform transient simulations until a steady state is reached. Subsequently, the steady state obtained by the transient simulation is used as the initial guess for direct steady state simulations at the same values of the flow parameters, and we perform continuation to a certain flow variable (e.g., Wi).

In order to trace the families of the steady solution branches of a parameter k (from the group [Wi, β, $ε$⁠, and B_R]), we employ the pseudo-arc-length continuation algorithm,⁵⁷ as implemented in the FEM framework by Varchanis et al.⁵⁸ According to this algorithm, given a steady-state solution vector (S) and an initial solution (S₀) at an initial parameter value (k₀) on a solution branch, the value of k for the next computational step is not defined a priori but is a function of the pseudo-arc-length Δs_a. The inner product between the derivative of the initial solution with respect to the pseudo-arc-length (⁠$Ṡ$⁠) and (S − S₀) is added to the derivative of k₀ with respect to the pseudo-arc-length (⁠$k̇0$⁠) and (k − k₀) and set to be equal to Δs_a. The resulting equation

is added to the total system of equations that are solved by the Newton–Raphson method. This handling enables efficient tracking of pitchfork and saddle-node bifurcations in the parameter space. Furthermore, using the bifurcation theory and the fact that the parameter continuation starts from a steady state obtained by transient simulations, we can determine whether a solution branch is stable or unstable.

The microfluidic cylinder device, which is the same as that employed in several prior experimental studies and described in detail therein,^35,38,39 was fabricated by selective laser-induced etching in fused silica glass using a commercial “LightFab” three-dimensional (3D) printer (LightFab GmbH, Germany).^59–61 This method enables production of extremely slender cylinders that retain rigidity owing to the high material modulus (≈75 GPa). The device has a width W = 400 μm and a cylinder radius R = 20 μm. The finite height of the channel (hence also the length of the cylinder in the direction normal to the plane of the page in Fig. 1) is H = 2000 μm. The channel aspect ratio is relatively high A_R = H/W = 5, which has been shown to provide good uniformity of the flow profile through much of the channel height,^38,39 i.e., a good approximation to the 2D flow assumed in the simulations. The blockage ratio B_R = 0.1 also matches that used in the majority of the simulated flows.

Three viscoelastic test solutions are studied in experiments that serve as comparisons to validate the numerical simulations with each of the three constitutive models. The first test fluid, which will be referred to as XG, is a very weakly elastic, but strongly shear-thinning solution of 0.1 wt. % of xanthan gum dissolved in an 82 wt. % mixture of glycerol and water (ρ = 1210 kg m⁻³, η_s = 56 mPa s). The second fluid, referred to as PEO4, is a viscoelastic solution with almost constant shear viscosity composed of 0.2 wt. % of an M_$w$ = 4 MDa PEO dissolved in a 50 wt. % mixture of glycerol and water (ρ = 1124 kg m⁻³, η_s = 5 mPa s). The final fluid, referred to as PEO8, is a viscoelastic and shear-thinning solution composed of 0.5 wt. % of an M_$w$ = 8 MDa PEO dissolved in pure water (ρ = 997 kg m⁻³, η_s = 0.9 mPa s). All reagents were obtained from Sigma-Aldrich. The rheology of the test fluids in steady shear was measured at 25 °C using an Anton Paar MCR 502 stress-controlled rotational rheometer fitted with a 50 mm diameter 1° cone-and-plate geometry, and the resulting flow curves are shown in Fig. 2(a). The shear rheology of the XG fluid is well described by the shear-thinning but inelastic C–Y model discussed in Sec. II C. The PEO4 test fluid is reasonably approximated by the constant viscosity but elastic FENE-CR model (Sec. II C) with the extensibility parameter set to L² = 5000. The final PEO8 test fluid is represented by the l-PTT model (Sec. II C) with the shear-thinning and elasticity parameters set to β = 0.05 and $ε$ = 0.05, respectively.

The relaxation times (λ) of the fluids at 25 °C were determined in a uniaxial extensional flow by measuring the diameter as a function of time of the liquid bridge generated in a capillary thinning extensional rheology device [Haake CaBER 1, Thermo Scientific, see Fig. 2(b)].^45,46 The CaBER device was fitted with plates of diameter D₀ = 6 mm, the initial gap between the plates was 1 mm, and the plates were separated to a final gap of 6 mm by linear displacement at a rate of 0.1 m s⁻¹. The measured relaxation times are λ = 0.006 s (XG), λ = 0.044 s (PEO4), and λ = 0.145 s (PEO8).

The test fluids are driven through the microfluidic cylinder device at precisely controlled volumetric flow rates Q using two neMESYS low pressure syringe pumps (Cetoni GmbH) with 29:1 gear ratios. One of the pumps is used to infuse fluid into the device, while the second pump withdraws fluid at an equal and opposite rate from the downstream outlet. The pumps are fitted with Hamilton Gastight syringes of appropriate volumes such that the specified “pulsation free” dosing rate is always exceeded even at the lowest imposed Q. Connections between the syringes and the microfluidic devices are made using non-compliant poly(tetrafluoroethylene) tubing. The average flow velocity in the microchannel is U = Q/WH.

The flow of the polymeric test solutions around the cylinder is quantified using microparticle image velocimetry (μ-PIV, TSI Inc., MN). The test fluids are seeded with a low concentration (c_p ≈ 0.02 wt. %) of 2 μm diameter fluorescent tracer particles (PS-FluoRed particles, Microparticles GmbH) with excitation/emission wavelength 530/607 nm. The midplane of the geometry (i.e., the plane bisecting the axis of the cylinder) is brought into focus on an inverted microscope (Nikon Eclipse Ti) with a 5×, numerical aperture NA = 0.15 Nikon Plan Fluor objective lens. Under these conditions, the measurement depth over which microparticles contribute to the determination of the velocity field is δ_m ≈ 125 μm (or ≈H/16),⁶² and the flow field is expected to be 2D over this region.³⁸ Particle fluorescence is induced by excitation with a dual-pulsed Nd:YLF laser with a wavelength of 527 nm and a time separation between laser pulses Δt. A high speed camera (Phantom MIRO) operating in the frame-straddling mode captures pairs of particle images in synchrony with the laser pulses. At each flow rate examined, the time Δt is set so that the average displacement of particles between the two images in each pair is ≈4 pixels. The flows examined appear steady on the spatiotemporal scale of the experiment, so at each flow rate, 50 image pairs are processed using an ensemble average cross-correlation PIV algorithm (TSI Insight 4G) in order to reduce noise. A recursive Nyquist criterion is employed with a final interrogation area of 16 × 16 pixels² to enhance the spatial resolution and obtain 2D velocity vectors u = (u, $v$⁠) spaced on a square grid of 12.8 × 12.8 μm. Further image analysis, generation of contour plots and streamline traces, is performed using the software Tecplot Focus (Tecplot, Inc., WA).

The Reynolds number, describing the relative strength of inertial to viscous forces in the flow experiments, can be defined in a number of different ways due to the various length scales in the flow domain. Considering the flow around the cylinder itself, we define Re_cyl = ρUR/η₀, where the fluid density ρ is assumed to be equal to that of the solvent and we use the zero-shear-rate viscosity η₀ as a representative viscosity since the cylinder is located on the channel centerline where the shear rate vanishes. By this definition, Re_cyl ≤ 0.05 at the highest applied values of U. However, since some of the fluids are significantly shear-thinning [Fig. 2(a)], it is also important to consider the Reynolds number in the regions between the cylinder and the channel walls, $Regap=0.5ρUW/η(γ̇w,gap)$⁠, where $γ̇w,gap$ is the nominal wall shear rate in the gap and $η(γ̇w,gap)$ is the shear-rate-dependent viscosity. Based on this definition, Re_gap approaches a maximum value of 0.25 at the highest imposed flow rates. Based on these two estimates of the Reynolds number (both <1) computed using contrasting length scales and fluid viscosities, we consider that inertial effects can be ignored, and the assumption of creeping flow in the numerical simulations is reasonable.

The Weissenberg number, describing the relative strength of elastic to viscous forces in the flow around the cylinders, is defined identically in both the experiments and the numerical simulations, i.e., $Wi=λU/R=λϵ̇$⁠, where $U/R=ϵ̇$ is a characteristic deformation rate for the flow, defined by the nominal velocity gradient in the cylinder wake. Strong elastic effects and extension-rate-hardening can be expected in the wake for Wi ≳ 1.

For ease of comparison between the numerical and experimental results that follow, lengths and velocities in the experimental data are also nondimensionalized by R and U, respectively. All nondimensionalized quantities in the results that follow are indicated by the superscript “^*,” as defined in Sec. II A.

First, we report on the flow of a shear-thinning but inelastic fluid around the confined cylinder. The fluid is modeled numerically by the C–Y GNF model (Sec. II C) and represented in experiments by the very weakly elastic XG solution (relaxation time λ = 6 ms, Sec. III B). Velocity fields around the cylinder determined by the respective simulations and experiments performed over a range of U are shown in Fig. 3. In the simulations, the velocity fields appear almost self-similar and the flow is, in fact, perfectly symmetric about both the x = 0 and y = 0 axes. In the experimental flow fields (right column of Fig. 3), we note the development of a minor degree of fore-aft asymmetry as the flow rate is increased, with a more extended region of low flow velocity observed downstream of the cylinder, which is not seen in the corresponding numerical simulations. Since this fore-aft asymmetry is already observed at flow rates where the Reynolds numbers in the experiments are very low [in Fig. 3(b), Re_cyl ≈ 10⁻⁴ and Re_gap ≈ 0.04] it is unlikely to be due to inertia and is most likely explained by the slight elasticity of the experimental fluid, which is not accounted for in the numerical model. However, in the experiments, as in the simulations, the flow remains essentially perfectly symmetric about the y = 0 axis over the whole range of accessible flow rates. The lateral asymmetry in the flow field around the cylinder, which has been reported for flows of wormlike micellar (WLM) solutions and for shear-thinning viscoelastic polymer solutions,^28,35,38,39 is not observed here for this shear-thinning but essentially inelastic fluid.

The features of the flow fields discussed above are presented in a more quantitative form by the flow velocity profiles provided in Fig. 4. There is a rather good agreement between the streamwise numerical and experimental profiles along y = 0 upstream of the cylinder [i.e., for x^* < 0, Fig. 4(a)] and between the general forms of the numerical and experimental streamwise velocity profiles taken along x = 0 [Fig. 4(b)]. Of course, the slightly higher maximum values of the flow velocity measured experimentally are expected since (in contrast to the simulation) the experimental geometry is not perfectly 2D but has a finite height. As mentioned, the small deviation between simulation and experiment downstream of the cylinder [i.e., for x^* > 0, Fig. 4(a)] is most likely explained by the slight elasticity of the XG test fluid.

In Fig. 5, we present additional numerical data from the C–Y model simulations in the form of the normalized principle stress difference, PSD^* = PSD(R/η₀U), where $PSD=(τxx−τyy)2+(2τxy)2$⁠. The range of average flow velocities extends significantly higher than that shown in Fig. 3, but clearly, the flow continues to remain perfectly symmetric about both the x = 0 and y = 0 axes. In addition, since the fluid is inelastic, there is no intense stress growth in the cylinder wake, which has been observed for strongly viscoelastic fluids and extensible models such as FENE-CR (e.g., Refs. 3 and 38), and this will be seen in Subsection IV B.

Here, we report results on the flow of non-shear-thinning but elastic fluids. The numerical simulations are performed using the FENE-CR model (Sec. II C), while experiments are conducted with the PEO4 polymer solution of almost constant shear viscosity (see Sec. III B).

Figure 6 shows a comparison between simulated and experimental velocity magnitude fields (where the extensibility parameter of the FENE-CR model is set to L² = 5000), for a range of flow rates and hence Wi. As Wi is increased, both the simulations and experiments show the growth of a greatly extended downstream wake of low velocity, not seen in the inelastic fluid (Fig. 3). This is attributed to the stretching of polymers (or dumbbells) in the wake, with an increase in the extensional stress, which locally viscosifies (extension-rate-hardens) the fluid and modifies the flow field.^3,38 In the present simulations, this strong stress growth with increasing Wi > 1 in the cylinder wake is illustrated in Fig. 7. Although the development of this downstream wake means that the flow becomes strongly asymmetric about x = 0 as Wi becomes high, laterally about y = 0, perfect symmetry is maintained in both the experiments and the simulations up to the highest Wi = 46 tested. There is no clear preference for the fluid to flow around either one or the other side of the cylinder. This observation is also consistent with earlier experiments using almost constant viscosity but highly elastic polystyrene solutions in a similar flow geometry.³⁸ We note that, in the numerical simulations presented here, the solutions become transient for Wi ≥ 10; in these cases, the fields presented in Figs. 6 and 7 are time averaged. This transience is not evident in the experiments. As discussed in Sec. III D, the flows examined experimentally appeared essentially steady given the accessible spatiotemporal scales and noise levels of the velocimetry method, and all the presented experimental velocity fields are time-averaged. We note that, for Wi ≥ 20, simulations with the FENE-CR model were extremely challenging, and sometimes, the numerical scheme diverged. In such cases, simulations were restarted using progressively smaller time steps until convergence was achieved.

Further details of the time-dependence for Wi ≥ 10 are provided by the results of transient simulations (see Fig. 8). In Fig. 8(a), the drag force on the cylinder [$FD=−∫ΓCexn:(−PI+τ+ηsγ̇)dΓC$⁠, where Γ_C represents the cylinder surface, e_x is the unit vector in the x-direction, and n is the unit normal to the cylinder surface] is plotted as a function of the time from the startup of flow at a few different Weissenberg numbers. At each imposed Wi, the drag force increases with time up to a plateau value. Although the plateau in F_D appears quite steady in each case, in fact, for Wi > 10, there are small fluctuations with time. The time dependence at higher values of Wi is more apparent when considering the streamwise stress τ_xx in the wake of the cylinder. This is extracted from the simulations at a point on the y = 0 axis 5.5 radii downstream of the cylinder centerpoint and plotted as a function of time in Fig. 8(b). For Wi = 20 and Wi = 46, τ_xx clearly varies erratically with time. Power spectra of the fluctuations do not show any characteristic frequencies (data not shown). The fluctuation in the stress arises due to a small side-to-side motion of the downstream wake about the y = 0 axis, which does, in fact, cause a very small time-dependent lateral asymmetry in the flow around the sides of the cylinder. Similar to prior works,^35,39 the degree of lateral asymmetry I is assessed using a dimensionless parameter that varies between −1 and 1,

where u₁ and u₂ are the streamwise fluid velocities measured near the midpoints between the cylinder and the top and bottom side walls, respectively (note: in simulations, we take exact values from the nearest mesh nodes: $u1=u|x*=0,y*=5.3846$ and $u2=u|x*=0,y*=−5.3846$⁠). By this measure, I = 0 implies a perfectly symmetric flow, while I = ±1 implies that all of the fluid passes on just one side of the cylinder.

The time-dependence in the lateral asymmetry for the FENE-CR model fluid at Wi = 46 is shown in Fig. 8(c). Note that the level of asymmetry is very small; here, |I| ∼ 10⁻⁴, whereas values |I| ∼ 1 have been measured in some experiments.^35,39 In addition, in Fig. 8(c), the asymmetry parameter varies about a value of zero and on a time average becomes almost completely negligible, as evident from the images provided in Figs. 6(c) and 7(c).

We further investigate the effect of the extensibility parameter L² on the response of the FENE-CR model around the cylinder. As shown in Fig. 9, higher L² results in an increase in the drag force on the cylinder [Fig. 9(a)] and, for a given Wi, an increasingly long downstream wake [Figs. 9(b)–9(d)]. Since higher L² implies increased deformability of the dumbbells, both effects are attributed to the resulting increase in elastic stress generated in the wake of the cylinder. However, varying L² does not qualitatively change the response of the model and no significant lateral asymmetry is observed even at the highest extensibility of L² = 10 000.

In this section, we report on the flow of fluids that exhibit both shear-thinning and elastic rheological responses. The simulations are performed using the l-PTT constitutive model (Sec. II C), and experiments are carried out using the PEO8 test fluid (Sec. III B).

Figure 10 provides a comparison between flow velocity fields obtained from the numerical simulations with the l-PTT model and from experiments performed with the PEO8 fluid over a range of imposed Wi. For the simulations, the shear-thinning and elasticity parameters have been set to β = 0.05 and $ε$ = 0.05, respectively. Clearly, there is good agreement between the experiments and the simulations, which both show rather symmetric flow fields at lower values of Wi [i.e., Wi = 1.5 and Wi = 8, Figs. 10(a) and 10(b)], but distinctly asymmetric flow fields (about both y = 0 and x = 0) at the two higher Wi shown [i.e., Wi = 46 and Wi = 61, Figs. 10(c) and 10(d)]. Experimentally (to date), only shear-thinning viscoelastic fluids have been observed to develop such laterally asymmetric flow states.^35,39 Using the newly developed FEM formulation to reach sufficiently high Wi,⁷ the l-PTT model, which incorporates both shear-thinning and elasticity in the constitutive equation, successfully replicates the experimentally observed lateral asymmetry. It is important to mention that the example flow fields presented in Fig. 10 have been selected in order to show consistent sequences of behavior in which the fluid passes the cylinder preferentially with high velocity on the side of positive y. However, the selection of this preference is quite random in both the experiments and the simulations. As will become evident during the presentation of later numerical results, there is equal probability that the fluid will flow preferentially around either side of the cylinder such that the asymmetric state is bistable, and there also appears to be no dependence of the favored configuration on the choice of model parameters.

The numerical simulations are used to look in greater detail at the flow fields and corresponding stress fields around the cylinder (Fig. 11) as Wi is increased beyond the onset of lateral asymmetry. We observe that, for lower Wi [e.g., Fig. 11(a)], the laterally symmetric flow is accompanied by the development of an elastic downstream wake along the y = 0 axis (similar in form to the wake observed for the non-shear-thinning FENE-CR model, Fig. 7). However, for Wi ≳ 20 [e.g., Fig. 11(b)], the flow field begins to exhibit lateral asymmetry, with a higher flow speed at negative y than at positive y. The asymmetry is also manifested by a distortion in the associated stress field, with the downstream wake being bent toward positive y. At even higher Wi = 46 [Fig. 11(c)], the lateral asymmetry in the flow and stress fields becomes significantly stronger.

In Fig. 12, we take the absolute value of I as a function of Wi in order to compare the asymmetry parameter obtained from the l-PTT simulations (lines) and the experiments (symbols) with the PEO8 fluid. Clearly, the numerical curves display a sudden transition to the asymmetric state, whereas for the experimental data, the growth of asymmetry with Wi is more gradual. This may be explained by the polydispersity of the experimental polymer sample. Such commercially available PEO samples typically have molecular weight distributions spanning around three orders of magnitude,⁶³ which could be expected to smear the transition (note that the experimental Weissenberg number is calculated based on the single value of the relaxation time determined via the CaBER measurement, but in reality, there is a wide spectrum of relaxation times present in the test fluid due to the molecular weight distribution). We would expect to obtain a better match to the experimental result by employing a multi-mode constitutive model in the simulations (see, e.g., Ref. 64). However, we choose a simple single mode model since it enables a more precise understanding of the role of the various material parameters on the instability mechanisms. In any case, Fig. 12 shows a reasonable agreement between the experiment and the simulations in terms of both the value of Wi_c for the onset of asymmetry (Wi_c ≈ 20–30) and also the maximum plateau value of |I| (around 0.6–0.7). Figure 12 also serves to demonstrate the numerical mesh convergence (see Table I for details of the various meshes employed in this study). The results are almost mesh-independent, in terms of both the asymmetry parameter and the stress in the downstream wake of the cylinder (Fig. 12, inset). In Table II, we give exact absolute values of the asymmetry parameter (|I|) estimated using the three meshes at various values of Wi. Since we reach up to very high Wi (i.e., relatively high flow velocities, U) in these simulations, we also confirm that inertia remains negligible by performing numerical simulations in mesh M2 that include the inertial terms in the momentum equation. As can be seen from Fig. 12, the result is almost indistinguishable from that under the creeping flow assumption. Mesh M2 was chosen for the majority of the simulations due to its slightly superior accuracy over mesh M1 but lower computational cost compared with mesh M3. Note that the reduction in the stress in the downstream wake for Wi > Wi_c shown in the inset of Fig. 12 is due to the deflection of the high stress region away from the y = 0 axis as the flow develops an increasing level of lateral asymmetry.

In this subsection, we investigate the influence of varying the shear-thinning (β) and elasticity (⁠$ε$⁠) parameters in the l-PTT model.

Figure 13 shows the effect of varying β for a fixed value of elasticity $ε$ = 0.05. In Fig. 13(a), we show the asymmetry |I| as a function of Wi for four different β values, and in Fig. 13(b), we show the corresponding flow velocity fields in the vicinity of the cylinder at Wi = 25. It is clear that as β is increased (i.e., shear-thinning is reduced), the onset Wi_c for asymmetry increases. In addition, for a given value of Wi, increasing β tends to reduce the magnitude of the asymmetry. For sufficiently weak shear-thinning (i.e., β = 0.1), the asymmetric flow state cannot be observed. Increasing the degree of shear-thinning in the fluid model tends to enhance asymmetric flows.

In Fig. 14, we show the effect of varying the elasticity of the l-PTT model for a fixed level of shear-thinning with β = 0.05. Figure 14(a) shows the asymmetry parameter |I| as a function of Wi for four different $ε$ values, and Fig. 14(b) shows the corresponding flow velocity fields in the vicinity of the cylinder at Wi = 25. In this case, it is clear that increasing $ε$ (i.e., reducing the elasticity) results in an increase in the onset Wi_c for asymmetry and, for a given Wi, a decrease in the magnitude of I. For very weak elasticity (i.e., $ε$ = 0.1), the asymmetric flow state is not observed. Increasing the elasticity of the model fluid tends to enhance asymmetric flows.

In Fig. 15(a), we present a stability diagram for the l-PTT model in the Wi–β parameter space. Stability boundaries between symmetric and asymmetric flow states are marked along contours of $ε$⁠. Following lines of constant $ε$ from high to low β (recall lower β implies stronger shear-thinning), Wi_c decreases (i.e., shear-thinning is destabilizing). As $ε$ is decreased for fixed β (recall lower $ε$ implies higher elasticity), Wi_c also decreases (i.e., elasticity is destabilizing). However, for a given value of Wi, the flow destabilizes at lower β as $ε$ is increased (or at lower $ε$ as β is increased). Clearly, there is an interplay between the shear-thinning and the elasticity. Low levels of shear-thinning can be compensated for by high levels of elasticity, and vice versa, although a certain amount of both shear-thinning and elasticity is required in order for the asymmetric flow state to be inducible. We collapse all our data from the l-PTT model onto a single curve in the plot shown by Fig. 15(b). Here, for each set of model parameters, we have extracted the values of U/R and $γ̇w,gap$ and computed the corresponding values for the extensional viscosity in the cylinder wake η_E,$w$ake and the shear viscosity on the side of the cylinder η_side. In Fig. 15(b), these are plotted against each other in the normalized form and fitted with the empirical relation $ηE,wake=36.4η0/((ηside/η0)−0.83−1)$⁠. We note that this function asymptotes to infinite extensional viscosity as η_side → η₀ (i.e., when there is no shear-thinning), and to zero extensional viscosity as η₀ → ∞ (i.e., as the fluid becomes infinitely shear-thinning).

Having identified a suitable fluid model with which to simulate the occurrence of the laterally asymmetric flows beyond critical Weissenberg number conditions, we may gain further insight into the problem by considering the effects of geometrical changes through varying the blockage ratio B_R. All of the simulations presented in this subsection are performed using the l-PTT model with β = $ε$ = 0.05. If the hypothesis regarding the instability mechanism presented in the Introduction and in previous papers^35,39 is correct, then a variation of B_R should modify the onset conditions for asymmetry. For example, for a fixed Wi > Wi_c, a manipulation of B_R could be made in order to shift the gap wall shear rate $γ̇w,gap=12U/W(1−BR)2$ onto different parts of the flow curve. The a priori expectation is that asymmetry will only occur if $γ̇w,gap$ lies within the shear-thinning region of the flow curve, and the fluid is sufficiently elastic. So far, a systematic study of this flow phenomenon in terms of the geometric parameter space has not been attempted, and it would be difficult to perform properly by an experimental approach, due to requiring a large number of different devices with fabrication challenges.

In Fig. 16(a), we report the effect of varying the blockage ratio by showing the magnitude of the asymmetry parameter |I| as a function of Wi for various values of B_R. Figure 16(b) shows flow velocity fields for fixed Wi = 25 at four different values of B_R. For increasing blockage ratios, Fig. 16(a) shows a reduction in Wi_c. In fact, in all four B_R cases, at Wi_c, the gap wall shear rate is almost the same at $γ̇w,gap≈100$ s⁻¹, which is close to the middle of the shear-thinning region of the flow curve (see Fig. 2). For a fixed Wi = 25, increasing the blockage ratio in the range shown (0.08 ≤ B_R ≤ 0.2) results in an increase in the magnitude of the lateral flow asymmetry. At first thought, this may be attributed to the resulting increase in the shear rate. However, the scenario becomes more complicated if the blockage ratio is continuously increased to higher values.

In Fig. 17, we track the solution branch by keeping Wi = 30, β = 0.05, and $ε$ = 0.05 constant and performing pseudo-arc-length continuation to B_R. At low B_R ≲ 0.07, the flow is stable in the laterally symmetric state (I = 0), even at this fairly appreciable value of Wi. As the blockage ratio is increased beyond a critical value B_R,c ≈ 0.07, the flow undergoes a forward (supercritical) bifurcation, the symmetric state destabilizes, and the asymmetry parameter becomes non-zero following the bistable asymmetric branch, where I can be either positive or negative. As B_R is further increased, |I| increases through a maximum, before decreasing and turning backward through a saddle-node bifurcation. The symmetric flow state restabilizes on a backward (subcritical) pitchfork bifurcation at the second higher critical blockage ratio B_R,c2 ≈ 0.26.

This interesting resymmetrization of the flow at high B_R as elasticity is kept constant is reminiscent of a resymmetrization observed in earlier experiments with a range of shear-thinning viscoelastic solutions of hydrolyzed poly(acrylamide) (HPAA).³⁵ In those experiments, the blockage ratio was fixed at B_R = 0.1 and the independent variation of Wi and $γ̇w,gap$ was not possible since both quantities scale with U. However, several of the HPAA solutions tested exhibited flow asymmetry over a limited range of Wi > Wi_c and recovered the symmetric state at even higher Wi. It was shown that the recovery of symmetry occurred as $γ̇w,gap$ approached the high shear rate Newtonian-like plateau of the flow curve. These observations were interpreted as follows: despite the fact that elastic effects were increasing due to the incrementing Weissenberg number, the flow regained symmetry because shear-thinning effects around the cylinder became too weak to sustain the imbalance.

In Fig. 18(a), we present the results of additional simulations carried out for smooth variations of B_R at four different values of Wi. At all four values of Wi, the fluid model behaves similarly, exhibiting both upper and lower critical values of B_R. A minor detail worthy of note is that, for Wi = 15, the pitchfork at B_R,c2 is not subcritical; the saddle-node bifurcation is absent at this Wi and the flow is resymmetrized via a second supercritical transition. Despite this detail, some general trends are evident: as Wi is increased, both B_R,c and B_R,c2 decrease, the range of B_R over which the symmetric flow is unstable tends to increase, and the maximum value of |I| also increases. Since the shear-thinning and elasticity parameters are fixed at β = 0.05 and $ε$ = 0.05, respectively, a variation of B_R at fixed Wi is equivalent to a simple variation of the shear rate. The curves shown in Fig. 18 can therefore be interpreted as meaning that, for each Wi, there is a certain range of shear rates (not too high and not too low) over which the symmetric flow can become unstable. According to our a priori stated hypothesis, this range of shear rates should lie on the shear-thinning portion of the flow curve. We computed the lower critical shear rate $γ̇w,gap,c$ and the upper critical shear rate $γ̇w,gap,c2$ corresponding to the bifurcation points B_R,c and B_R,c2 (respectively) shown in Fig. 18(a). The lower and upper critical shear rates are marked in Fig. 18(b) on the flow curve of the l-PTT model fluid. Clearly, by this nominal measure of the shear rate, the lower critical shear rates all coincide at $γ̇w,gap,c≈100$ s⁻¹, about midway into the shear-thinning region of the flow curve. For all Wi, the shear rates over which the symmetric flow is unstable all lie on the most shear-thinning part of the flow curve. The upper critical shear rate reaches a progressively higher value as the Weissenberg number is increased. This happens because at higher values of Wi, elastic effects are stronger, and asymmetric flows can be supported for weaker shear-thinning effects (see Fig. 15). Apparently, regardless of the value of Wi and B_R, as the wall shear rate approaches the high shear rate viscosity plateau, the flow will necessarily resymmetrize as shear-thinning effects become negligible.

The observations described above are largely consistent with the proposed instability mechanism (discussed in the Introduction) based on an initial random lateral perturbation of the downstream elastic wake for Wi > 1 causing a small asymmetry in the flow around the sides of the cylinder that can be reinforced only if $γ̇w,gap$ is in the shear-thinning regime. However, so far, we have not identified the nature of the perturbation of the elastic wake that is necessary to initiate the instability. The fact that the flow loses symmetry for an almost constant value of $γ̇w,gap,c≈100$ s⁻¹, regardless of the imposed Wi, strongly suggests a universal scaling for the onset of instability. A well-known criterion for the onset of purely elastic (i.e., inertialess) flow instabilities proposed by McKinley and co-workers considers the generation of elastic tensile stresses along curving streamlines and can be expressed as^29–31 

where $R$ is the characteristic radius of curvature of the streamline and τ₁₁ is the tensile stress in the stream direction. When, at some locality in the flow field, the dimensionless product of parameters on the left-hand side of Eq. (16) exceeds a critical value $Mcrit2$ (which can depend on the geometrical configuration), then the flow becomes prone to instability originating from that location. For flows around a confined cylinder, simple scaling arguments for the value of $R$ near the downstream stagnation point indicate that instability will arise for 1/Wi_c ∼ B_R (with linear stability analysis giving $Mcrit2≈37$⁠).²⁹ As shown by the inset of Fig. 18(b), our numerical data for B_R,c vs 1/Wi follow the predicted linear scaling very well. This suggests that the initial perturbation to the flow field that drives the onset of flow asymmetry is due to the accumulation of elastic tensile stress along the strongly curving streamlines passing near the downstream stagnation point of the cylinder.

We have used a recently developed finite element flow solving scheme that circumvents the high Weissenberg number problem in order to simulate 2D inertialess flows of various non-Newtonian fluid models around circular cylinders, reaching up to high deformation rates. Our objective in modeling this flow was to better understand how different aspects of fluid rheology (i.e., shear-thinning and viscoelasticity) affect the manifestation of a lateral flow asymmetry around the cylinder that has been observed in some recent experiments.

Using a shear-thinning but inelastic model (a Carreau–Yasuda generalized Newtonian fluid), we were unable to reproduce the lateral flow asymmetry and the flow around the cylinder remained both fore-aft and laterally symmetric at all tested flow rates. Complementary experiments with a shear-thinning and only weakly elastic xanthan gum solution with a similar rheology to the model also remained stable and symmetric up to the highest achieved Weissenberg number Wi = 30, and the measured flow fields were in good agreement with the model predictions.

A non-shear-thinning but highly elastic fluid model (the FENE-CR model) was then employed in order to understand the effect of elasticity in isolation. The simulation results were compared with those of experiments performed using an elastic and almost constant viscosity PEO solution. In this case also, neither the model nor the experimental fluid showed the onset of the steady lateral flow asymmetry observed in earlier experimental works. However, in both the simulations and the experiments, we observed the growth of a strong elastic wake downstream of the cylinder due to the polymer (or dumbbell) extension in the region of the stagnation point for Wi ≳ 1. For higher Weissenberg numbers Wi ≥ 10, the simulations yielded time dependent solutions (not apparent in the time-averaged experimental results) in which the downstream wake fluctuated weakly from side to side, resulting in a very small, time-dependent lateral flow asymmetry. However, the fluctuation was about the symmetric state, and on a time-average, the flow appeared completely symmetric. We also tested the effect of varying the extensibility parameter in the FENE-CR model, but with no qualitative change in the behavior.

Finally, we examined the flow of a shear-thinning and elastic fluid model around the cylinder using the l-PTT model. We compared the flow response of this model with a shear-thinning and elastic PEO solution of similar rheology. In this case, the model showed the growth of an elastic wake downstream of the cylinder for Wi ≳ 1 (similar to the elastic, non-shear-thinning fluid) and the onset of a steady lateral flow asymmetry above a critical value of the Weissenberg number Wi_c ∼ 20. We also observed good agreement between the flow patterns and degrees of flow disparity around the sides of the cylinder resulting from the experiments and the simulations.

These results strongly support the notion that both shear-thinning and elasticity of the fluid are necessary for the flow asymmetry to occur. Having obtained a model that could reproduce the experimental observations accurately, we were then able to use the l-PTT model to systematically vary the strength of the shear-thinning and elasticity parameters and to study their effects on the onset conditions of the asymmetry and its development with increasing Wi.

For a fixed degree of elasticity $ε$ in the l-PTT model, decreasing the level of shear-thinning caused an increase in the value of Wi_c and a reduction in the degree of flow asymmetry for Wi > Wi_c. Similarly, for a fixed value of the shear-thinning β, reducing the elasticity also caused an increase in Wi_c and a reduction in the degree of asymmetry for Wi > Wi_c. For very low values of either shear-thinning or elasticity, the onset of asymmetry was completely suppressed. Our analysis shows that there is a complex interplay between the shear-thinning and the elasticity of the fluid. A degree of both is required in order for the instability to occur, but to some extent, weak shear-thinning can be compensated for by strong elasticity, and vice versa.

Our final set of simulations using the l-PTT model involved fixing the values of β and $ε$ while the blockage ratio of the channel B_R was varied. This effectively allowed us to vary only the shear rate in the channel, and specifically around the sides of the cylinder, which we assume is the most relevant region of shearing in the flow domain. Such a parameter study has not been performed previously for this flow. For a given Wi, we found the flow remains stable and symmetric for low values of B_R and becomes bistable and asymmetric as B_R is increased beyond a point of supercritical bifurcation at B_R,c. At B_R,c, the wall shear rate in the gaps either side of the cylinder is always around $γ̇w,gap,c≈100$ s⁻¹, close to the middle of the shear-thinning portion of the flow curve. Our analysis shows that B_R,c scales linearly with 1/Wi, as expected for the onset of a purely elastic instability driven by elastic tensile stresses along the tightly curved streamlines passing close to the downstream stagnation point. We believe that this elastic instability provides the initial perturbation that allows the large scale lateral flow asymmetry to develop (if the fluid is also sufficiently shear-thinning). As B_R is increased beyond B_R,c, the degree of flow asymmetry grows through a maximum, before symmetry is again recovered at another bifurcation point at B_R,c2. By computing the shear rate on the sides of the cylinder corresponding to the second bifurcation point, we showed that the region of bistable asymmetric flow lies on the most shear-thinning portion of the flow curve. Due to the increasing influence of elasticity, the bistable region spans a greater range of shear rates as Wi is increased but never impinges on the high shear rate pseudo-Newtonian viscosity plateau where the shear-thinning vanishes.

Our results help to paint a complete picture of the nature of this novel flow instability, which may be relevant to understanding a number of phenomena observed in shear-thinning viscoelastic flows (e.g., selection of preferred paths through complex geometries such as porous media). With reference to experiments carried out with wormlike micellar solutions, in future work, we plan to extend our numerical study to constitutive models with non-monotonic flow curves and more complex relaxation dynamics in order to understand how shear-banding and the breaking/reformation of micelles affect the instability.

The data that support the findings of this study are available within the article.

S.V. gratefully acknowledges the support of Greece and the European Union (European Social Fund, ESF) through the operational program “Human Resources Development, Education and Lifelong Learning” in the context of the project “Strengthening Human Resources Research Potential via Doctorate Research” (MIS-5000432), implemented by the State Scholarships Foundation (IKY). J.T. acknowledges funding from the LIMMAT Foundation under the project “MuSiComPS.” C.C.H., A.Q.S., and S.J.H. gratefully acknowledge the support of the Okinawa Institute of Science and Technology Graduate University (OIST) with subsidy funding from the Cabinet Office, Government of Japan, and also funding from the Japan Society for the Promotion of Science (JSPS, Grant Nos. 17K06173, 18K03958, 18H01135, and 20K14656) and the Joint Research Projects (JRPs) supported by the JSPS and the Swiss National Science Foundation (SNSF).