1. Introduction
Rotational flows with curvilinear streamlines are routinely encountered in many applications, such as transport and handling of fluids, dispensing, spin coating, flows in centrifuges, extruders, flows around objects, lubrication and journal bearing flows, polymer processing, as well as rheometry. The torsional flow between co-axial parallel disks or between a cone and a plate are canonical examples in which a unidirectional shear flow with curved streamlines is generated; hence, they are commonly employed in rheometry and as canonical systems in which to study viscoelastic flow stability. Even in the limit of small cone angles, steady secondary motions and flow instabilities can develop (Larson Reference Larson1992; Shaqfeh Reference Shaqfeh1996) owing to nonlinearities stemming from interactions among streamline curvature, fluid inertia and non-Newtonian fluid properties.
Specifically, in a cone-and-plate geometry with cone radius $R$, a cone angle $\theta$ and rotating at speed $\varOmega$, a Newtonian fluid with a constant viscosity $\eta _s$ and density $\rho$ experiences the onset of inertially driven toroidal secondary motions beyond a critical rotation rate (Olagunju Reference Olagunju1997), and at higher rates, inertial turbulence, when the centrifugal force becomes dominantly larger than the viscous force (Sdougos, Bussolari & Dewey Reference Sdougos, Bussolari and Dewey1984). However, highly elastic ‘Boger fluids’ with a constant viscosity $\eta _0$ exhibit a time-dependent transition to an unstable state characterized by an enhanced and fluctuating shear stress (initially interpreted as anti-thixotropic behaviour: Jackson, Walters & Williams Reference Jackson, Walters and Williams1984), when sheared beyond a critical value of the dimensionless rotational speed, or Deborah number $De=\tau _s\varOmega$, where $\tau _s$ is the shear relaxation time (McKinley et al. Reference McKinley, Byars, Brown and Armstrong1991), which is constant for a Boger fluid. The unstable flow manifests as radially propagating logarithmic spiral vortices (Öztekin, Brown & McKinley Reference Öztekin, Brown and McKinley1994). Furthermore, linear stability analysis of the cone-and-plate flow using different constitutive models reveals the critical role of the dimensionless geometry parameter $1/\theta$ (Olagunju Reference Olagunju1997) on the onset of elastic instabilities, and this sensitivity has also been confirmed by experiments (Öztekin et al. Reference Öztekin, Brown and McKinley1994; McKinley et al. Reference McKinley, Öztekin, Byars and Brown1995). Thus, the instability is a function of both the dimensionless rotational speed, i.e. $De=\tau _s\varOmega$, and a dimensionless form of the shear rate $\dot {\gamma } = \varOmega /\theta$, i.e. the Weissenberg number $Wi=\tau _s\dot {\gamma }$. A comparison of critical conditions in the $Wi$–$De$ state diagram with linear stability results indeed validates this codependency (Öztekin et al. Reference Öztekin, Brown and McKinley1994; McKinley et al. Reference McKinley, Öztekin, Byars and Brown1995). Thorough reviews of these earlier studies of viscometric flow instabilities can be found from Larson (Reference Larson1992) and Shaqfeh (Reference Shaqfeh1996).
Pakdel & McKinley (Reference Pakdel and McKinley1996) argued that this type of instability generically arises due to the non-uniform stretching of the polymer molecules in a curvilinear shearing flow, which amplifies the elastic ‘hoop stress’ in the fluid. This enters the radial momentum balance and amplifies radial velocity perturbations, making the torsional flow of viscoelastic fluids unstable under sufficiently strong driving conditions. They proposed an instability criterion in terms of the product of $Wi$ and $De$ exceeding a critical magnitude denoted $M_{c}$, and showed that this could be written in the generic form $De Wi = (\tau _s\varOmega ) ({\tau _s\dot {\gamma }} ) > M_{c}^2$. More generally, we can write $De=\tau _s U /\mathcal {R}$, with $\mathcal {R}$ being the (geometry-dependent) characteristic radius of the streamline curvature and $U$ being the characteristic streamwise fluid velocity. The ratio $\mathcal {R}/U$ gives the characteristic convective time of a flow experiment. Thus, to take into account, more generally, the curvature of a two-dimensional flow and the action of tensile stress difference along the streamlines, this criterion for the onset of instability can be rewritten as (McKinley, Pakdel & Öztekin Reference McKinley, Pakdel and Öztekin1996)
where $U$ is a characteristic streamwise fluid velocity, $\mathcal {R}$ is the geometry-dependent characteristic radius of curvature of the streamline, $N_1$ is the first normal stress difference in the fluid and $\sigma =\eta _0\dot {\gamma }$ is the shear stress. Furthermore, linear stability analysis of the steady base flow of an Oldroyd-B fluid in a cone-and-plate geometry (Olagunju & Cook Reference Olagunju and Cook1993; Olagunju Reference Olagunju1995, Reference Olagunju1997), as well as in a Taylor–Couette geometry (Schaefer, Morozov & Wagner Reference Schaefer, Morozov and Wagner2018), shows that the viscosity ratio parameter $\beta _P=\eta _P/\eta _0$ (where $\eta _P$ is the polymer contribution to the fluid viscosity and $\eta _{\infty }$ is the Newtonian plateau viscosity at high shear rates, such that $\eta _0 = \eta _P+\eta _{\infty }$) significantly affects the nature of the elastic instability. The effect of $\beta _P$ can be readily incorporated in the condition for purely elastic instability ((1.1)) by substituting the Oldroyd-B result $N_1 = 2\eta _P\tau _s\dot {\gamma }^2$ and $\sigma =(\eta _P+\eta _{\infty })\dot {\gamma }$, which gives
In dilute polymer solutions, as $\beta _P \to 0$, the critical shear rate $\dot {\gamma }_c=\varOmega _c/\theta$ required for the onset of instability diverges, in accordance with experiments. More recently, Schiamberg et al. (Reference Schiamberg, Shereda, Hu and Larson2006) have studied the effect of changing the viscosity ratio $\beta _P$ on the onset of secondary motion and the evolution towards a fully developed nonlinear state commonly referred to as ‘elastic turbulence’ in a parallel plate geometry using a polyacrylamide Boger fluid.
Theoretically, it should be possible to modify the criterion in (1.1) for the onset of purely elastic instability in shear-thinning viscoelastic fluids by allowing $\tau _s$ and $\eta _P$ to both be shear-rate dependent and writing them as functions of the applied shear rate $\dot {\gamma }$. However, the reduction in the viscosity due to shear thinning means a concomitant increase in inertial effects, which can systematically modify the purely elastic instabilities observed in Boger fluids at very low Reynolds numbers, $Re \ll 1$. This inherent nonlinear coupling between inertia and elasticity in shear-thinning viscoelastic fluids makes the corresponding torsional flows more challenging to understand, especially the critical conditions for the onset of elasto-inertial instabilities. Hence, there have been relatively few studies elucidating the effect of shear thinning on viscoelastic flow stability. Dutcher & Muller (Reference Dutcher and Muller2013), Schaefer et al. (Reference Schaefer, Morozov and Wagner2018) and Lacassagne, Cagney & Balabani (Reference Lacassagne, Cagney and Balabani2021) have investigated shear-thinning-mediated elasto-inertial transitional pathways in Taylor–Couette geometries. Similar observations have also been made in the case of the flow of a shear-thinning viscoelastic fluid through a tube (Chandra et al. Reference Chandra, Mangal, Das and Shankar2019). Each of these studies reveals the inherent coupling of fluid elasticity (which tends to destabilize the base shearing flow) and inertia (which tends to restabilize the unsteady flow of viscoelastic fluids). The resulting stability diagrams are best represented by three-dimensional plots in terms of dimensionless parameters characterizing the geometry, fluid elasticity and the inertia of the flow (Dutcher & Muller Reference Dutcher and Muller2009).
Very recently, Datta et al. (Reference Datta2022) reviewed our current understanding of the broader topic of viscoelastic flow instabilities and elastic turbulence. They suggest representing the critical conditions for the onset of viscoelastic flow instabilities in the $Wi$–$Re$ plane. In this plane, a set of exploratory experiments with a given rate-independent viscoelastic fluid in a fixed geometry trace a line with a slope given by the elasticity number $El=Wi/Re$ (see figure 4 of Datta et al. Reference Datta2022) eventually intersecting with a corresponding stability boundary demarcating the critical conditions for the onset of instability. Tracing the critical conditions using different viscoelastic fluids and flow geometries enables exploration of the entire phase space, and the various flow states observed will showcase the transitions in the instability mechanisms as inertial, elastic and geometric effects are varied. In the present work, we investigate the development of time-dependent instabilities in a family of shear-thinning viscoelastic polyisobutylene (PIB) solutions of different PIB concentrations ($c_P$) using a set of cone-and-plate geometries of different cone angles ($\theta$) and radii ($R$) to elucidate the combined effects of fluid elasticity, shear thinning, inertia and geometry on the onset of torsional viscoelastic flow instabilities.
2. Materials and methods
We use polymeric solutions of PIB (MW $\approx 10^6\ {\rm g}\ {\rm mol}^{-1}$) dissolved in a paraffinic oil (GADDTAC, Lubrizol Inc.). We perform all our measurements at a constant temperature $T=20\,^\circ$C. The Newtonian base oil has a solvent viscosity $\eta _s=18.07\ {\rm mPa}\ {\rm s}$ at $20\,^\circ$C. Dilute solution viscometry determines the polymer intrinsic viscosity to be $[ \eta ] = 3.69\ {\rm dL}\ {\rm g}^{-1}$ and this gives the critical overlap concentration of the polymer solute as $c^* \simeq 0.77/[ \eta ] = 0.23$ wt.% (Graessley Reference Graessley1980). The solutions were all measured to have a constant density $\rho =873.1\ {\rm kg}\ {\rm m}^{-3}$ and surface tension $\varGamma = 29.7\ {\rm mN\ m}^{-1}$. We vary the dissolved concentration of polymer in the solution to change the viscoelastic properties. We work with three semi-dilute ($c_P > c^*$) solutions: 3 wt.%, 2 wt.% and 1 wt.%, respectively, and one close to $c^*$ with $c_p= 0.3$ wt.%. In table 1, we summarize the key material properties for the family of polymeric fluids used in the study. These polymeric solutions are all shear thinning to various extents, with their viscoelasticity decreasing at lower concentrations, as shown in figure 1 and table 2.
The inelastic Cross model, $\eta (\dot {\gamma })=\eta _\infty +(\eta _0-\eta _\infty )/[1+(\lambda \dot {\gamma })^n]$, does an excellent job at describing the viscosity as a function of shear rate for all solutions (for best-fit parameter values, see figure 1a and table 2). Here, $\eta _0$ and $\eta _\infty$ are viscosities in the limit of zero and infinite shear rate, respectively, $\lambda ^{-1}$ is a measure of the characteristic shear rate for the onset of shear thinning, and $n$ is the shear-thinning index. The extent of shear thinning in these viscoelastic solutions can be quantified in several ways. One direct metric is the ratio $\eta (\dot {\gamma })/\eta _0$. A more useful ratio that also provides consistency with earlier analysis, as we show later, is the dimensionless function $\beta _P(\dot {\gamma }) \equiv \eta _P(\dot {\gamma })/\eta (\dot {\gamma })=[\eta (\dot {\gamma }) -\eta _{\infty }]/\eta (\dot {\gamma })$, which quantifies the relative polymer contribution to the total (rate-dependent) solution viscosity at a given shear rate.
Shear-thinning rheology also means that the first normal stress coefficient $\varPsi _1 = N_1/\dot {\gamma }^2$ and the relaxation time of the solutions are no longer material constants, but become rate-dependent functions denoted by $\varPsi _1(\dot {\gamma })$ and $\tau _s(\dot {\gamma })$. From the functional form of the upper-convected derivative, the first normal stress difference can be related to the shear stress and the relaxation time of the fluid through the expression $N_1(\dot {\gamma }) \simeq 2\tau _s\dot {\gamma }\sigma$ (Bird, Armstrong & Hassager Reference Bird, Armstrong and Hassager1987), which, after a slight rearrangement, gives $N_1(\dot {\gamma })/\sigma (\dot {\gamma }) \simeq 2\tau _s\dot {\gamma }$. This ratio of the first normal stress difference $(N_1)$ to the shear stress in the sheared fluid $(\sigma)$ is termed the stress ratio $S_R$ and is a direct quantitative measure of nonlinear viscoelastic effects in these polymeric solutions (see figure 1b). The importance of nonlinear elastic effects at the onset of flow instability in each fluid is evident as the stress ratio $S_R=N_1/\sigma \gtrsim 10$. One can also eliminate $\dot {\gamma }$ by substituting $\dot {\gamma } = \sigma /\eta$ in the expression above, which gives $N_1 \approx 2(\tau _s/\eta )\sigma ^2$. So, one can anticipate on theoretical grounds that $N_1$ will vary quadratically with $\sigma$. The quadratic dependence of the first normal stress difference $N_1$ on the shear stress is well known for polymer solutions (Lodge, Al-Hadithi & Walters Reference Lodge, Al-Hadithi and Walters1987; Binding, Jones & Walters Reference Binding, Jones and Walters1990) and has recently been documented in viscoelastic emulsions as well (Kibbelaar et al. Reference Kibbelaar, Deblais, Briand, Hendrix, Gaillard, Velikov, Denn and Bonn2023). This key observation that $N_1 \sim \sigma ^2$ will be crucial later in enabling us to incorporate shear-thinning rheology effects in a unified critical instability criterion (cf. § 3.4).
Furthermore, in semi-dilute polymer solutions, $\tau _s$ and $\eta$ are both typically power law functions of the polymer concentration $c_P$ (Heo & Larson Reference Heo and Larson2005). As a result, the modulus $G_c(c_P) \approx \eta (c_P)/\tau _s(c_P)$ only increases weakly with increasing the polymer concentration $c_P$. Hence, $N_1 \approx 2(\tau _s/\eta )\sigma ^2 \approx \sigma ^2/G_c(c_P)$ may be expected to only decrease weakly with the polymer concentration $c_P$ in the fluids used. We plot the first normal stress difference $N_1$ as a function of the shear stress $\sigma$ in figure 2 to test this prediction. We find that not only $N_1 \sim \sigma ^2$, but this relationship also holds irrespective of the polymer concentration except when the flow becomes unsteady (hollow symbols). Thus, we can estimate the steady-state values of $N_1(\dot {\gamma })$ from measurements of steady-state shear stress alone as
where $A$ is a constant coefficient obtained from the rheological measurements shown in figure 2. For the fluids used in this study, we find $A=0.29$ Pa$^{-1}$.
To incorporate rate-dependent rheological effects in the criterion for the onset of elastic instabilities (cf. (1.1)), we use the more general definition of the Weissenberg number $Wi=N_1/2\sigma$ (White Reference White1964), which can be simplified using (2.1) as
Shear-thinning rheology also results in a rate-dependent characteristic shear relaxation time $\tau _s(\dot {\gamma })$, which enters (1.1). This relaxation time can be directly evaluated from experimental measurements of the first normal stress difference and the polymer contribution to the shear stress $\sigma _P=\eta _P(\dot {\gamma })\dot {\gamma }$. We can thus define $\tau _s(\dot {\gamma })=N_1(\dot{\gamma})/(2\eta _P(\dot {\gamma })\dot {\gamma }^2)$. After a slight rearrangement along with substituting $\sigma = \eta (\dot {\gamma })\dot {\gamma }$ and $\beta _P(\dot {\gamma }) = \eta _P(\dot {\gamma })/\eta (\dot {\gamma })$, we can thus write
This general expression for the rate-dependent characteristic relaxation time of a viscoelastic fluid can be combined with (2.1) if it is validated experimentally for a given material system. For completeness, we note that this definition correctly reduces to the expression used in the literature (derived from the Oldroyd-B formulation) for the case of a non-shear-thinning Boger fluid; $N_1= 2\tau _s \sigma _P \dot {\gamma }=2 \tau _s \eta _P \dot {\gamma }^2$, where $\eta _P=\eta _0-\eta _s, \lim _{\dot {\gamma } \to 0}\varPsi _{1} = \varPsi _{1,0} = 2 \tau _s \eta _P$ and $\lim _{\dot {\gamma } \to 0}\tau _s(\dot {\gamma }) = \tau _s$. The rate-dependent shear relaxation time $\tau _s(\dot {\gamma })$ defined in (2.3) can also be used (if desired) to define a rate-dependent Deborah number $De=\tau _s(\dot {\gamma })\varOmega$, which incorporates shear thinning in the relaxation time. In addition, the normal stress difference ratio $\psi = -\varPsi _{2,0}/\varPsi _{1,0}$ for the PIB solutions used in this study fall in the range $0.205\unicode{x2013}0.243$ as measured from rod-climbing rheometry (More et al. Reference More, Patterson, Pashkovski and McKinley2023). Here, $\lim _{\dot {\gamma } \to 0 }{\varPsi _{2}} = \varPsi _{2,0}$ with $\varPsi _{2}$ being the second normal stress coefficient. We have included additional rheological characterization of the various PIB solutions used in this study, including small amplitude oscillatory shear measurements, in the supplementary material available at https://doi.org/10.1017/jfm.2024.254.
3. Results and discussion
As shown in figure 1(a), a sudden jump in the measured apparent steady-state shear viscosity is observed for all solutions beyond a sufficiently high shear rate when sheared in a cone-and-plate geometry. These measurements are not erroneous but are manifestations of the onset of time-dependent secondary motion beyond a critical (composition and geometry dependent) shear rate $\dot {\gamma }_c$. As a result, the stresses in the samples also suddenly increase due to large velocity disturbances. This is the origin of the so-called ‘anti-thixotropic’ transition (Jackson et al. Reference Jackson, Walters and Williams1984), observed and analysed previously for constant viscosity PIB Boger fluids (McKinley et al. Reference McKinley, Byars, Brown and Armstrong1991; Öztekin et al. Reference Öztekin, Brown and McKinley1994; Schiamberg et al. Reference Schiamberg, Shereda, Hu and Larson2006). In the present study, we analyse the same instability for viscoelastic fluids with pronounced shear thinning, which inherently drives a transition of the instability mechanism from purely elastic to elasto-inertial as we gradually decrease the polymer concentration. Finally, we use the following conventions for the governing rate-dependent dimensionless parameters: (1) a subscript $_0$ denotes dimensionless parameters determined using material function values in the zero-shear-rate limit, i.e. without incorporating rate-dependent viscosity effects, e.g. $Re_0=\rho \varOmega R^2/\eta _0, Wi_0 = \tau _s \dot {\gamma }$; (2) the absence of any subscript denotes dimensionless parameters determined with the incorporation of rate-dependent viscosity effects, e.g. $Re = \rho \varOmega R^2/\eta (\dot {\gamma }), Wi = \tau _s(\dot {\gamma }) \dot {\gamma }$; and (3) a subscript $_c$ denotes the critical condition $\dot {\gamma } \to \dot {\gamma }_c$, e.g. $Re_c=\rho \varOmega _c R^2/\eta (\dot {\gamma }_c), Wi_c = \tau _s(\dot {\gamma _c}) \dot {\gamma _c}$.
3.1. Time-dependent instability, power spectra and unsteady flow visualization
The experimentally measured time-evolving shear stress $\sigma (t,\dot {\gamma })$ and first normal stress difference $N_1(t,\dot {\gamma })$ for a constant shear rate experiment in a 40 mm 2$^\circ$ cone-and-plate geometry are shown in figure 3 for two different PIB fluids. The transient responses of $\sigma (t)$ and $N_1(t)$ are typical of viscoelastic fluids; at short times, a rapid stress increase with an overshoot (and initial quadratic growth in $N_1$) is observed. At longer times, constant steady-state values are observed only when the imposed shear rate is less than a critical value $\dot {\gamma }_c$. Beyond this point, both $\sigma (t,\dot {\gamma })$ and $N_1(t,\dot {\gamma })$ increase rapidly to a new fluctuating state similar to the earlier observations for Boger fluids (McKinley et al. Reference McKinley, Byars, Brown and Armstrong1991). Calculations of the apparent viscosity and first normal stress difference from averaging these enhanced time-dependent measurements result in the sudden jump observed in figure 1. However, the time dependence of the unstable flow of the shear-thinning viscoelastic fluid has notable differences compared with a Boger fluid.
As we reduce the PIB concentration, the viscosity of the solutions decreases and inertial effects become increasingly important. The competition between the effects of shear thinning and the effects of inertia results in the onset of instability at a lower $Wi_c$ than the elastic shear-thinning case of the 3 wt.% PIB solution. In addition, we observe a gradual decay in the mean amplitude of the stress fluctuations for the less concentrated PIB fluids due to slow irreversible sample ejection from the geometry edge as well as possible viscous heating effects (Calado, White & Muller Reference Calado, White and Muller2005). As a result, the eventual steady-state value reached after the fluctuations in the $\sigma (t)$ and $N_1(t)$ response die out can be lower than the quasi-steady-state value achieved after the initial transient is completed (at times $t\sim O(10)s$). Thus, reducing the PIB concentration gradually shifts the instability mechanism from purely elastic to elasto-inertial, and shear thinning plays a crucial role in determining the critical conditions for the onset of instability. This transition in the instability mechanism becomes clear from the power spectra of the shear stress $\sigma$ (denoted $E_{\sigma \sigma }$) and the first normal stress difference $N_1$ (denoted $E_{N_1{N_1}}$) obtained by taking Fourier transforms of $\sigma (t)$ and $N_1(t)$ after the onset of time-dependent fluctuations. The corresponding power spectra for the fluid response (at $\dot {\gamma }>\dot {\gamma }_c)$ are presented in figures 3(c) and 3(d), respectively. The slope of the power spectra progressively changes from $-$3 to $-$2 for $E_{\sigma \sigma }$ but remains unchanged for $E_{N_1N_1}$ as we gradually reduce the PIB concentration below 3 wt.%. Furthermore, when plotted as a function of an appropriate dimensionless frequency $\tau _s(\dot {\gamma }) f$, the measured spectra fall onto two distinct curves: one for the more concentrated solutions (3 and 2 wt.%) and a separate curve for the less concentrated solutions (1 and 0.3 wt %). This change in $E_{\sigma \sigma }$ appears to be a distinctive feature of a transition in the underlying instability mechanism (Steinberg Reference Steinberg2022).
Figure 4 shows instantaneous visualizations of the unsteady flow after the onset of instability in a 40 mm $2^\circ$ cone geometry, and the corresponding space–time diagrams (kymographs) depict the evolution of the unsteady flow with time. The instability originates at the edge of the conical fixture and spirals towards the centre (as can be seen from the fine dark streak spiralling inwards). As the perturbations propagate towards the centre of the cone, they dissipate, which can be clearly seen in the kymographs; the undulating fluctuations are most intense near the edge but become dimmer towards the centre of the cone. The progressive ejection of material from the gap due to time-dependent fluctuations is evident at long times.
3.2. Hysteresis and determination of critical conditions for onset of instability
Hysteresis has been observed for Boger fluids in the past following a sudden jump down in shear rate from a value larger than the critical value (after the induction of a time-dependent flow state) to a shear rate lower than the critical value that had initially been determined in a step ramp-up protocol resulting in the onset of time-dependent flow (McKinley et al. Reference McKinley, Byars, Brown and Armstrong1991). These observations are consistent with a subcritical Hopf bifurcation; however, directly accessing the hysteretic state and assessing its extent is easier if the control variable is switched to be the imposed stress. In the present work, we perform both controlled stress sweeps as well as stepped ramps up in imposed shear rate measurements.
The critical values we report for the onset of instability, i.e. $\dot {\gamma }_c$, are determined by loading a fresh fluid sample and then performing a series of step increases in the shear rate. At each new shear rate, the evolution in the shear stress and the normal stress difference is followed for a long time (${\approx }30$ min) to observe the possible onset of a time-dependent unstable flow at a constant applied shear rate (see figure 3a). The lowest shear rate at which such a transition was observed is then determined to be the critical shear rate. These critical values are indicated by the arrows in figure 5. We call this a ‘stepped shear rate ramp-up’ protocol. However, the transition might be hysteretic in the sense that if one has to ‘step down’ from an unstable flow rate, the flow might only become stable at a shear rate lower than the critical value obtained in a ‘stepped ramp-up’ experiment. To explore the presence and extent of flow hysteresis, we perform a continuous slow ramp-up and subsequent ramp-down in the shear stress (i.e. a saw-tooth stress profile) using the same geometry on a stress-controlled rheometer. In this protocol, we continuously increase the shear stress imposed on the sample (linearly in time) until the flow becomes unstable, and then slowly and continuously reduce the applied stress. By ‘slow’, we mean that the rate of increasing the applied stress is much slower than the stress relaxation time for the fluids; here, we use a ramp rate of $1/5\tau _s$. In the absence of flow instability, such a protocol should trace out the steady flow curve predicted by the fit of the Cross model for each fluid (shown by the dashed lines in figure 5). Hysteresis in the apparent flow curve of shear stress versus shear rate can thus be directly confirmed by departures from this flow curve, as shown in figure 5.
It is evident from figure 5 that the presence of strong viscoelastic shear-thinning reduces the magnitude of the flow hysteresis observed in the purely elastic case (for the 2 wt.%, 3 wt.% fluids); however, careful measurements using the stepped shear rate ramp-up and ramp-down protocols, and subtraction of the steady state flow stress $\sigma _s$, reveal that a small amount of hysteresis is still present (e.g. see the inset of figure 5 for the 2 wt.% fluid). In the case of the inertio-elastic instability observed for the 0.3 and 1 wt.% fluids, it is clear that the extent of flow hysteresis is very pronounced (with ${\rm \Delta} \dot {\gamma }_c = \dot {\gamma }_c^{up} - \dot {\gamma }_c^{down} \approx 100\ {\rm s}^{-1}$).
3.3. Mapping the instability transition in the $Wi$–$Re$ plane
The changes in the magnitude of the stress fluctuations and the power law slope of the power spectra suggest a transition in the underlying mechanism behind the onset of instability as we change the polymer concentration. Experiments with a range of cone-and-plate geometries ranging from $\{R,\theta \}=\{2~\textrm {cm}, 1^\circ \}$ to $\{1.25~\textrm {cm}, 5^\circ \}$ also show that the critical conditions vary systematically with changes in the geometric parameter $1/\theta$. These changes can be understood more quantitatively by projecting the critical conditions for the onset of instability onto the $Wi$–$Re$ plane (Datta et al. Reference Datta2022) as illustrated in figure 6(a). The transition from purely elastic to elasto-inertial can be clearly demonstrated by considering the evolution in the elasticity number $El=Wi/Re$ with shear-thinning effects (solid lines). The elasticity number quantifies the relative magnitude of viscoelastic and inertial effects in the $Wi$–$Re$ state diagram. The trajectory followed by highly elastic fluids with a constant viscosity $\eta _0$ (i.e. Boger fluids) is represented in figure 6(a) by dashed lines with constant slopes ($El_0=Wi_0/Re_0$). By including the effect of shear thinning through the rate-dependent viscosity of the PIB solutions, we can conclude that shear thinning shifts the critical Reynolds number $Re_c$ for the onset of instability to a substantially higher value (at a given Weissenberg number). The value of the elasticity number at the onset of instability $El_c$ reduces by four orders of magnitude (from $\approx$800 to 0.02) as we reduce the concentration of the PIB solution by one order of magnitude (from 3 wt.% to 0.3 wt.%), signifying the rapid rise in inertial effects compared with viscoelastic effects.
At low $Re$, shear thinning increases the stability of the torsional shear flow of viscoelastic fluids against the onset of purely elastic instabilities. For example, the Boger fluid used by Schiamberg et al. (Reference Schiamberg, Shereda, Hu and Larson2006) exhibited purely elastic instability beyond $Wi_c \approx O(1)$ compared with the values of $Wi_c \sim O(10)$ required for an elastic but strongly shear-thinning PIB solution. However, reducing the concentration of the PIB solutions also decreases the extent of shear thinning, i.e. the magnitude of changes in $\beta _P$ are reduced, and the onset of instability shifts to $Wi_c \sim O(1)$. Thus, reducing the extent of shear-thinning (or maintaining $\beta _P$ closer to unity) drives the earlier onset of elastic instability when represented in terms of $De$ or $Wi$.
The gradual increase in inertial effects is another crucial factor contributing to shifts in the onset conditions of the instability with a reduction in the PIB concentration. It is observed that the presence of elasticity makes the flow unstable at much lower Reynolds numbers ($Re_c \lessapprox 100$) compared with the critical value of $\rho \varOmega _c R^2/\eta$ required for the onset of purely inertial turbulence in cone–plate flows of Newtonian fluids. Finally, the shaded region, which is included solely as a guide for the eye, hints at common trends and the possible collapse of all the critical conditions onto a single master curve. In the following section, we seek a unified critical condition spanning purely elastic to elasto-inertial regimes for a range of polymer concentrations and flow geometries.
3.4. Bridging elasticity and inertia: a unified instability criterion for the onset of instability in torsional flows of shear-thinning viscoelastic fluids
The governing role of shear thinning in stabilizing the flow of the PIB solutions can be seen from the magnitudes of the decrease in $\beta _P(\dot {\gamma })$ shown in figure 6(b) with shear rate and/or decreasing PIB concentration. In addition, (2.3) shows that the shear relaxation time is also a shear-thinning function of the shear rate and this plays an important role in determining the onset of instability. These nonlinear effects of shear thinning can be incorporated heuristically in the critical conditions for the onset of a purely elastic instability by incorporating the rate-dependence of the relaxation time $\tau _s(\dot {\gamma })$ from (2.3) into (1.2):
In a cone and plate geometry with cone angle $\theta$ and radius $\mathcal {R}=R$ rotating with a rate $\varOmega$, we get the characteristic velocity $U=\varOmega R$ and shear rate $\dot {\gamma }=\varOmega /\theta$. Using these characteristic values and the general definition of the Weissenberg number $Wi = N_1/2\sigma$, a further simplification of (3.1) can be obtained for a cone-and-plate geometry:
Thus, (3.2) gives a critical condition for the onset of purely elastic instability incorporating the effects of shear thinning, viscoelasticity and changes in flow geometry. We note that the critical condition of (1.2) for the onset of purely elastic instabilities in constant viscosity Boger fluids (McKinley et al. Reference McKinley, Pakdel and Öztekin1996) can be recovered from (3.2) by substituting the well-known Oldroyd-B results for $N_1 = 2\eta _P\tau _s\dot {\gamma }^2$ and $\sigma =\eta _0\dot {\gamma }$. Curves of neutral stability consistent with (3.2) can also be drawn in the $Wi - \beta _P(\dot {\gamma })/\theta$ plane, where the effects of shear thinning can be incorporated by shifting the geometric parameter $1/\theta$ by an amount depending on the shear-thinning parameter $\beta _P(\dot {\gamma })$ (Öztekin et al. Reference Öztekin, Brown and McKinley1994). These curves of neutral stability (dotted lines) incorporating shear-thinning and geometry effects in the purely elastic critical instability criterion along with the experimentally measured critical conditions are presented in figure 7(a) using the values of $\beta _P({\dot {\gamma _c}})$ shown in figure 6(b). We can immediately conclude that the critical conditions shown in figure 6(a) would only be shifted vertically by this scaling. This modified critical stability condition still does not incorporate the effects of fluid inertia.
A unified critical condition must ultimately be deduced by performing a rigorous linear stability analysis of the elasto-inertial flow problem. In the absence of any such existing analysis, our empirical measurements can be harnessed to guide the formulation of a unified critical condition. Purely elastic instability arises due to nonlinearities in the fluid constitutive equations, while purely inertial Newtonian turbulence arises from the nonlinearities in the advective term of the equation of motion. The distinct origins of the two sources of nonlinearity suggest that from the viewpoint of infinitesimal perturbations, the two destabilizing terms may be coupled together so that instability ensues when some combination of elastic effects and inertial effects becomes larger than a threshold or a critical value.
The effects of nonlinear elasticity and shear thinning on the steady base torsional shearing flows are already included in (3.1); this purely elastic criterion just needs to be augmented by the appropriate nonlinear term to incorporate the effects of fluid inertia in determining the onset of instability. In addition, this unified criterion should: (1) recover the critical instability condition for the onset of purely elastic instability in the absence of inertia; (2) recover the critical instability condition for the onset of purely inertial instability in the absence of elasticity; and (3) predict a smooth transition between the two asymptotic regimes, as observed in figure 6(a).
Detailed linear stability analyses (Joo & Shaqfeh Reference Joo and Shaqfeh1994; Öztekin et al. Reference Öztekin, Brown and McKinley1994) consider infinitesimal perturbations to the base flow velocity and stress variables $(\boldsymbol {v}^0, \boldsymbol {\sigma }^0)$ as well as the corresponding velocity and stress gradients $(\boldsymbol {\nabla } \boldsymbol {v}^0, \boldsymbol {\nabla }\boldsymbol {\sigma }^0)$ present in the specific flow field of interest. In the inertialess limit, the perturbed velocity and stress fields can be written in the dimensionless form $\boldsymbol {v} = \boldsymbol {v}^0 + Wi \boldsymbol {v}' + O(Wi^2)$ and $\boldsymbol {\sigma } = \boldsymbol {\sigma }^0 + Wi \boldsymbol {\sigma }' + O(Wi^2)$ (where the stress tensor has been non-dimensionalized with the characteristic viscous stress $\sim \eta _0 U/\mathcal {R}$). Substituting these forms into the governing equations and collecting terms at the first order results in a complex eigenvalue problem involving terms of the form $Wi(\boldsymbol {\nabla } \boldsymbol {v}^{0} \boldsymbol {\cdot } \boldsymbol {\sigma }'), Wi (\boldsymbol {v}^0 \boldsymbol {\cdot } \boldsymbol {\nabla }\boldsymbol {\sigma }')$, etc. (see Joo & Shaqfeh Reference Joo and Shaqfeh1994, Shaqfeh Reference Shaqfeh1996 for details). From a scaling viewpoint, the resulting stability criterion is always in the form of (1.1), or equivalently,
Similarly, for the onset of Newtonian curved streamline instabilities, the perturbations can be written in the form $\boldsymbol {v} = \boldsymbol {v}^0 + Re \boldsymbol {v}' + O(Re^2)$ and the eigenvalue problem arises from coupling between the advective terms of the form $Re \boldsymbol {v}^0 \boldsymbol {\cdot } \boldsymbol {\nabla } \boldsymbol {v}', Re \boldsymbol {v}' \boldsymbol {\cdot } \boldsymbol {\nabla } \boldsymbol {v}^0$, etc. From a scaling viewpoint, the resulting relevant dimensionless group that arises is the Görtler number $G = Re (\delta /\mathcal {R})^{3/2}$ (Saric Reference Saric1994), which couples streamline curvature parametrized by a characteristic radius of curvature $\mathcal {R}$, inertial effects parametrized by a Reynolds number and a small parameter $\epsilon = \delta /\mathcal {R}$ parametrized by the ratio of boundary layer thickness $\delta$ to the radius of curvature $\mathcal {R}$. In a Taylor–Couette flow, however, the relevant dimensionless group that arises is the Taylor number $Ta=Re^2\epsilon ^3$, which also couples inertial effects parametrized by Reynolds number and a small parameter $\epsilon = h/\mathcal {R}$, where $h$ is the gap between two concentric cylinders. The dimensionless grouping on the left-hand side of (3.3) can be viewed as a viscoelastic Görtler number (Pakdel & McKinley Reference Pakdel and McKinley1996). These arguments all suggest that inertial effects can be parametrized by combining the Reynolds number and a flow geometry-dependent small parameter. So, to quantify inertial nonlinearities, we use the dimensionless ratio of the dynamic pressure $\rho {U}^2$ compared with the shear stress $\sigma$, i.e.
where $C$ is some critical value of this ratio. We note that this ratio is of the same functional form as the dimensionless stress ratio $N_1/\sigma$, which determines the purely elastic instability. Furthermore, as Renardy (Reference Renardy2000) has noted in strongly nonlinear elastic flows streamline tension (which scales with $N_1$), gives rise to a modified form of the inviscid Euler or Bernoulli equation. This scaling in (3.4) automatically gives rise to a combination of Reynolds number and a small parameter to quantify the inertial effects, similar to Görtler and Taylor numbers, as will be evident in the discussion below.
A full inertioelastic linear stability analysis must, ultimately, consider perturbations arising from both elastic effects and inertial effects as well as cross-coupling terms that scale as $Re Wi \boldsymbol {v}'$. The resulting eigenvalue problem will be extremely complex, but from a dimensional viewpoint, it must embody the two limits, purely elastic and purely inertial, discussed above, as well as relevant cross-terms. The simplest possible functional form that captures these limits as well as cross-coupling between the elastic stresses and inertial perturbations is a linear combination of the two dimensionless groupings in (3.3) and (3.4). This linear combination can be written in the generic form
where $\alpha$ is a dimensionless weighting parameter and $\tilde {M}_c$ is the critical disturbance magnitude (as modified by the inertial effects). Expanding the square in the second equality gives the locus of the critical conditions. Further simplification of (3.5) in a cone-and-plate geometry can be obtained by using (3.2) and substituting the characteristic velocity $U=\varOmega {R}$ and $\dot {\gamma } = \varOmega /\theta$ so that at the critical condition $\dot {\gamma } \to \dot {\gamma }_c$, we require
where a rate-dependent Reynolds number $Re=\rho \varOmega {R^2}/\eta (\dot {\gamma })$ and also the (small) geometry parameter $\theta$ appears naturally in the inertial term. This new generalized criterion in (3.6) asymptotically recovers the purely elastic instability condition ((3.2)) when inertia is negligible, i.e. $Re \to 0$, and also predicts the onset of secondary motion due to inertial effects in Newtonian fluids in a cone-and-plate geometry when the product $Re \theta$ exceeds a critical value consistent with earlier perturbation analyses and experiments (Fewell & Hellums Reference Fewell and Hellums1977; Sdougos et al. Reference Sdougos, Bussolari and Dewey1984).
In lieu of formal linear stability analysis, this unified critical condition can be validated empirically with experimental measurements. We use the result from (2.1), which indicates that $N_1$ varies quadratically with $\sigma$ to calculate the critical Weissenberg number ${Wi_c = A\sigma (\dot {\gamma }_c)/2}$ with $A=0.29$ obtained from the regression of the experimental measurements presented in figure 2. If the condition suggested by (3.6) is true, then plotting data for the onset of instability in the $Wi^2 \theta /\beta _P(\dot {\gamma })$ versus $Re \theta$ phase space should lie on a single curve. We find that this is indeed true for a suitable value of the dimensionless weighting parameter $\alpha =2.29$ obtained using a two-parameter best fit of the experimental data to the first expression in (3.6) and the result is shown by the solid line figure 7(b). The collapse of the data over three orders of magnitude in the elasticity number for four different viscoelastic fluids with varying shear-thinning strengths and for three different conical geometries strongly corroborates the functional form of the unified critical condition for the elasto-inertial instability given in (3.6). The least square curve fitting of (3.6) to experimental measurements of the critical conditions gives $\tilde {M}_c \approx 5.46$.
A further check of the unified critical condition (3.6) can be performed by checking if $\tilde {M}_c \approx 5.46$ obtained from our experimental measurements is consistent with critical conditions in the limiting cases of purely elastic instability of a Boger fluid and a purely inertial instability of a Newtonian fluid in a cone-and-plate geometry. We therefore project the dimensionless criterion represented by (3.6) (with $\alpha \simeq 2.29$ and $\tilde {M}_c \simeq 5.46$) in the modified $Wi$–$Re$ plane to incorporate the effects of shear thinning and flow geometry as shown in figure 7(b). Our heuristic blending rule provides a good collapse of all the experimental data shown in figure 6(a). A typical non-shear-thinning Boger fluid has $\beta _P({\dot {\gamma }}) \equiv \beta \simeq 0.4$ and $Re \to 0$, which recovers $Wi_c \approx 60$ consistent with previous experimental observations (McKinley et al. Reference McKinley, Byars, Brown and Armstrong1991). However, for a Newtonian fluid, (3.6) predicts the onset of the unsteady secondary flow in a cone-and-plate geometry at $Re \theta \lesssim 5.68$, which is consistent with, but smaller than, the value deduced from early visualization experiments of Sdougos et al. (Reference Sdougos, Bussolari and Dewey1984). This scaling form can, in principle, be readily extended to other geometries by the appropriate determination of the radius of curvature of the streamlines, the characteristic flow velocity and the identification of the appropriate characteristic shear rate (McKinley et al. Reference McKinley, Pakdel and Öztekin1996). For example, in a parallel plate geometry, the small geometric factor would be $H/R$, where $H$ is the gap height between the plates.
In the spirit of Occam's razor, we have proposed the simplest possible physically reasonable criterion (additively combining the destabilizing effects of elasticity and inertia). We hope that our results will serve as a motivator for a more detailed stability analysis that fully considers the cross-coupling of perturbation terms between the equation of motion and the nonlinear viscoelastic constitutive equation.
4. Conclusions
We have investigated the onset of time-dependent instabilities in the torsional shearing flow between a cone and a plate with $1^{\circ } \leq \theta \leq 5^{\circ }$ for a range of viscoelastic shear-thinning fluids. The onset of instability depends on the coupling between the fluid elasticity, inertia, shear thinning and flow geometry. Specifically, shear thinning in the fluid rheology stabilizes the unidirectional torsional base flow, while inertia destabilizes it. Furthermore, the systematic changes in the critical conditions and the change in the slope of the power spectra $E_{\sigma \sigma }$ suggests a transition in the instability mechanism from purely elastic to elasto-inertial as we reduce the polymer concentration (and the extent of shear thinning). We show how to represent this multidimensional flow instability transition by constructing a state diagram of critical conditions in the $Wi$–$Re$ plane. Finally, we propose a unified criterion (3.5) to predict the onset of elasto-inertial instability by augmenting the existing criterion for purely elastic instability to include shear thinning and coupling to perturbations arising from inertial effects in the Cauchy momentum equation. We validate the form of this unified criterion empirically with our experimental measurements of the critical conditions for a range of viscoelastic fluids and different conical geometries. The data collapse over a wide range of conditions corroborates the final functional form of (3.6) at least for a cone-and-plate geometry. Extending this unified criterion to other torsional flow geometries (like the parallel plate geometry and Taylor–Couette flow) using an appropriate Görtler number would provide further verification of this representation.
Supplementary material
Supplementary material are available at https://doi.org/10.1017/jfm.2024.254.
Acknowledgements
The authors would like to thank Lubrizol Inc. for financial support and for providing the components required to formulate the polymeric fluids used in this study.
Declaration of interests
The authors report no conflict of interest.