Hostname: page-component-74656f9884-j7cmm Total loading time: 0 Render date: 2026-11-08T05:13:42.778Z Has data issue: false hasContentIssue false
  • English
  • Français

Opposing roles of shear-thinning on successive instabilities in Taylor–Couette flow

Published online by Cambridge University Press:  03 November 2026

Fan Xu Open the ORCID record for Fan Xu [Opens in a new window] ,
Zifeng Li  and
Ning Yang Open the ORCID record for Ning Yang [Opens in a new window]
Fan Xu
Affiliation:
State Key Laboratory of Mesoscience and Process Engineering, Institute of Process Engineering, Chinese Academy of Sciences , Beijing 100190, PR China
Zifeng Li
Affiliation:
State Key Laboratory of Mesoscience and Process Engineering, Institute of Process Engineering, Chinese Academy of Sciences , Beijing 100190, PR China School of Chemical Engineering, University of Chinese Academy of Sciences, Beijing 100049, PR China
Ning Yang*
Affiliation:
State Key Laboratory of Mesoscience and Process Engineering, Institute of Process Engineering, Chinese Academy of Sciences , Beijing 100190, PR China School of Chemical Engineering, University of Chinese Academy of Sciences, Beijing 100049, PR China
*
Corresponding author: Ning Yang, nyang@ipe.ac.cn
Rights & Permissions [Opens in a new window]

Abstract

Content of image described in text.

The intrinsic effect of shear-thinning on successive instabilities in narrow-gap Taylor–Couette flow is investigated using direct numerical simulation (DNS) with axially periodic boundary conditions, thereby eliminating endwall effects. The Carreau model is employed with three shear-thinning indices ($n = 0.52$, $0.62$, $0.68$), matching the experimental parameters of Cagney and Balabani (Phys. Fluids, 2019, vol. 31, p. 053102). Shear-thinning exerts opposing influences on the two instabilities. For the primary transition from circular Couette flow to Taylor vortex flow, enhanced shear-thinning (decreasing $n$) is destabilising, with the critical Reynolds number $\textit{Re}_{c,1}$ decreasing from $101$ to $94$. In contrast, the secondary transition to wavy vortex flow is strongly stabilised, with $\textit{Re}_{c,2}$ increasing from $115$ to $142$ over the same range of $n$. This opposing behaviour arises from a competition between two consequences of radial viscosity variation: the reduced viscosity near the inner wall promotes centrifugal disturbances (destabilising), while the elevated viscosity in the bulk gap damps inter-vortex oscillations (stabilising). Comparison with the non-monotonic experimental measurements of Cagney and Balabani (Phys. Fluids, 2019, vol. 31, p. 053102) shows that the observed discrepancies originate from endwall-induced wavelength constraints and vortex drift rather than intrinsic rheological effects. These results establish the intrinsic bifurcation structure of shear-thinning Taylor–Couette flow in the narrow-gap configuration and clarify the distinct roles of rheology and confinement in successive flow transitions.

JFM classification

Non-Newtonian Flows: Rheology Instability: Shear-flow instability

Information

Type
JFM Papers
Information
Journal of Fluid Mechanics , Volume 1046 , 10 November 2026 , A37
Check for updates
Copyright
© The Author(s), 2026. Published by Cambridge University Press

1. Introduction

Taylor–Couette (TC) flow, the motion of fluid between two coaxial cylinders with one or both rotating, has served for over a century as a canonical system for investigating hydrodynamic stability and transition to turbulence (Couette Reference Couette1890; Taylor Reference Taylor1923). For Newtonian fluids with a fixed outer cylinder, the flow is governed by the inner Reynolds number italic Re equals rho omega Subscript i Baseline r Subscript i Baseline d divided by mu $\textit{Re}=\rho \omega _{i}r_{i}d/\mu$ . At low italic Re $\textit{Re}$ , the flow is purely azimuthal circular Couette flow (CCF) (Couette Reference Couette1890). Beyond a critical value italic Re Subscript c comma 1 $\textit{Re}_{c,1}$ , centrifugal instability gives rise to steady, axisymmetric toroidal Taylor vortex flow (TVF) (Taylor Reference Taylor1923), which in turn becomes unstable to time-dependent wavy vortex flow (WVF) at a secondary critical Reynolds number italic Re Subscript c comma 2 $\textit{Re}_{c,2}$ (Andereck et al. Reference Andereck, Liu and Swinney1986). This rich bifurcation cascade has motivated extensive studies on the influence of non-Newtonian rheology, particularly shear-thinning, a nonlinear decrease in viscosity with shear rate, mu tilde ModifyingAbove gamma With dot Superscript n minus 1 $\mu \sim \dot \gamma ^{n-1}$ ( n less than or equals 1 $n\le 1$ ), which is ubiquitous in polymer solutions, suspensions and biological fluids (Mueller et al. Reference Mueller, Llewellin and Mader2010; Larson & Desai Reference Larson and Desai2015; Majji et al. Reference Majji, Banerjee and Morris2018). In TC flow, the spatial variation of viscosity necessitates a reference scale; while various choices exist (Masuda et al. Reference Masuda, Horie, Hubacz, Ohta and Ohmura2017), the nominal shear rate ModifyingAbove gamma With dot Subscript n o m Baseline equals omega Subscript i Baseline r Subscript i Baseline divided by d $\dot \gamma _{nom}=\omega _{i}r_{i}/d$ is most common (Cagney & Balabani Reference Cagney and Balabani2019; Cagney et al. Reference Cagney, Lacassagne and Balabani2020; Elçiçek & Güzel Reference Elçiçek and Güzel2020). For the primary instability (CCF right arrow $\rightarrow$ TVF) in narrow gaps ( eta equals r Subscript i Baseline divided by r Subscript o Baseline greater than or equals 0.8 $\eta =r_i/r_o\ge 0.8$ ), linear stability analyses and numerical studies consistently predict that enhanced shear-thinning (decreasing n $n$ ) reduces italic Re Subscript c comma 1 $\textit{Re}_{c,1}$ , provided the Reynolds number is defined using the zero-shear-rate viscosity (Lockett et al. Reference Lockett, Richardson and Worraker1992; Coronado-Matutti et al. Reference Coronado-Matutti, Souza Mendes and Carvalho2004; Ashrafi Reference Ashrafi2011; Alibenyahia et al. Reference Alibenyahia, Lemaitre, Nouar and Ait-Messaoudene2012). This destabilisation is attributed to reduce near-wall viscosity, which amplifies local centrifugal effects, and is accompanied by a profound modification of the TVF structure, including vorticity concentration at outflow boundaries and a strong asymmetry between inward and outward radial jets (Topayev et al. Reference Topayev, Nouar, Bernardin, Neveu and Bahrani2019; Cagney & Balabani Reference Cagney and Balabani2019).

In contrast, the influence of shear-thinning on the secondary instability (TVF right arrow $\rightarrow$ WVF) remains poorly understood and experimentally contentious. In a narrow-gap configuration ( eta equals 0.83 $\eta =0.83$ ) with aspect ratio upper Gamma equals 12.97 $\varGamma =12.97$ , Cagney & Balabani (Reference Cagney and Balabani2019) reported a non-monotonic dependence of italic Re Subscript c comma 2 $\textit{Re}_{c,2}$ on n $n$ ( italic Re Subscript c comma 2 Baseline equals 664 $\textit{Re}_{c,2} = 664$ for n equals 0.52 $n = 0.52$ , italic Re Subscript c comma 2 Baseline greater than 874 $\textit{Re}_{c,2} \gt 874$ for n equals 0.62 $n = 0.62$ and italic Re Subscript c comma 2 Baseline equals 826 $\textit{Re}_{c,2} = 826$ for n equals 0.68 $n = 0.68$ ), with the WVF transition vanishing for intermediate n $n$ , alongside slow axial vortex drift and spatially varying wave amplitudes. In another finite-length experiment with a similar radius ratio ( eta equals 0.8 $\eta =0.8$ ) but a lower aspect ratio ( upper Gamma equals 6.0 $\varGamma =6.0$ ), Kawaguchi et al. (Reference Kawaguchi, Tano and Saito2022) studied two shear-thinning fluids ( n equals 0.7 $n=0.7$ and 0.52 $0.52$ ) and found that italic Re Subscript c comma 2 $\textit{Re}_{c,2}$ dropped significantly from 1330 (Newtonian) to 608.7 for n equals 0.7 $n=0.7$ and 93.1 for n equals 0.52 $n=0.52$ , exhibiting a monotonic dependence on n $n$ , yet they observed pronounced axial elongation of end vortices and drift of vortex centres, features absent in their Newtonian counterpart. Meanwhile, Elçiçek & Güzel (Reference Elçiçek and Güzel2020) observed a direct transition from TVF to modulated wavy flow in xanthan gum solutions, bypassing the classic WVF regime entirely. These contradictory experimental outcomes strongly implicate endwall-induced finite-length effects, which constrain the axial wavelength and suppress the amplitude of the azimuthal waves. These contradictory experimental outcomes strongly implicate endwall-induced finite-length effects, which constrain the axial wavelength and suppress the amplitude of the azimuthal waves. Previous theoretical and numerical works by Cole (Reference Cole1976), Jones (Reference Jones1985) and Topayev et al. (Reference Topayev, Nouar and Dusek2022) have demonstrated that such confinement effects can significantly delay or alter the secondary bifurcation, obscuring the intrinsic rheological stability of the shear-thinning fluid.

Therefore, two fundamental questions arise. (a) Does shear-thinning intrinsically stabilise or destabilise the secondary instability? (b) What physical mechanism underlies this effect? To isolate the fundamental role of shear-thinning on the two successive bifurcations in a narrow gap, we perform high-resolution direct numerical simulations (DNS) with axially periodic boundary conditions, thereby removing the influence of endwalls that are always present in experiments, and match the geometry of Cagney & Balabani (Reference Cagney and Balabani2019) ( eta equals 0.83 $\eta =0.83$ ). This configuration enables a direct assessment of the intrinsic stability characteristics free from extrinsic wavelength constraints. The present work addresses two specific questions: whether shear-thinning intrinsically stabilises or destabilises the secondary instability, and what the underlying physical mechanism is. The remainder of this paper is organised as follows. We first briefly describe the numerical method and setting in § 2. Section 3 presents and analyses the results for radial viscosity variation and reveals the mechanism for the opposing effects of shear-thinning on primary and secondary instabilities. Finally, the work is concluded in § 4.

Figure 1.

(a) Sketch of the TC system. (b–d) Verification of the Carreau model implementation at Re = 40, corresponding to the CCF regime. Symbols denote the viscosity computed from the DNS at various shear rates for n equals 0.52 $n=0.52$ , 0.62 $0.62$ , and 0.68 $0.68$ ; solid lines represent the theoretical Carreau model (2.5) with parameters from table 1.

A diagram and three line graphs depict the Taylor-Couette flow system and the verification of the Carreau model implementation.
Figure 1. Long description

Panel A: A diagram of the Taylor-Couette (TC) system. The diagram shows a cylindrical setup with an inner cylinder of radius ri, an outer cylinder of radius ro, and a length L. The inner cylinder rotates with an angular velocity ωi, and the distance between the cylinders is denoted as d. Panel B: A line graph showing the verification of the Carreau model implementation at Re = 40, corresponding to the CCF regime. The x-axis is labeled λcβ and the y-axis is labeled μ(β)/μ0. The graph includes orange circles representing the viscosity computed from the DNS at various shear rates for n = 0.52, and a solid line representing the theoretical Carreau model. Panel C: A line graph showing the verification of the Carreau model implementation for n = 0.62. The x-axis is labeled λcβ and the y-axis is labeled μ(β)/μ0. The graph includes blue circles representing the viscosity computed from the DNS at various shear rates, and a solid line representing the theoretical Carreau model. Panel D: A line graph showing the verification of the Carreau model implementation for n = 0.68. The x-axis is labeled λcβ and the y-axis is labeled μ(β)/μ0. The graph includes red circles representing the viscosity computed from the DNS at various shear rates, and a solid line representing the theoretical Carreau model.

2. Numerical method and setting

The TC system consists of two concentric cylinders, as shown in figure 1(a). In the present study, the outer cylinder is at rest, the inner cylinder is rotating at angular velocity omega Subscript i $\omega _i$ and thus driving the flow. Both inner and outer walls are smooth and set as the no-slip boundary condition. To eliminate endwall effects inherent in experimental setups and isolate the intrinsic stability characteristics of the flow, axially periodic boundary conditions are imposed in the axial direction. The radius ratio is fixed at eta equals r Subscript i Baseline divided by r Subscript o Baseline equals 0.83 $\eta = r_{i}/r_{o} = 0.83$ , matching the experimental configuration of Cagney & Balabani (Reference Cagney and Balabani2019) to facilitate direct comparison. In addition, this narrow-gap configuration minimises curvature effects and allows the radial viscosity variation to be isolated from curvature-induced complexities; the base flow velocity profile is approximately linear, which facilitates interpretation. For wider gaps ( eta less than or equals 0.7 $\eta \le 0.7$ ), however, curvature effects become significant, and shear-thinning may interact with them to qualitatively modify the instability mechanisms. The axial periodicity length is set to upper L equals 2 pi d $L = 2\pi d$ , corresponding to an aspect ratio upper Gamma equals upper L divided by d equals 2 pi $\varGamma = L/d = 2\pi$ . This choice accommodates exactly three pairs of Taylor vortices, based on the experimentally observed wavelength of approximately 2.16 d $2.16d$ for the Newtonian case (Cagney & Balabani Reference Cagney and Balabani2019).

The working fluid is assumed to be incompressible and purely shear-thinning, with no viscoelastic effects. The flow is governed by the continuity and Navier–Stokes equations:

(2.1) bold nabla bold dot bold italic u equals 0 comma \begin{equation} \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u}=0, \end{equation}
(2.2) StartFraction partial differential bold italic u Over partial differential t EndFraction plus bold nabla bold dot left parenthesis bold italic u bold italic u right parenthesis equals minus StartFraction 1 Over rho EndFraction bold nabla p plus StartFraction 1 Over rho EndFraction bold nabla bold dot bold italic tau comma \begin{equation} \frac {\partial \boldsymbol{u}}{\partial t}+\boldsymbol{\nabla }\boldsymbol{\cdot }(\boldsymbol{uu})=-\frac {1}{\rho }\boldsymbol{\nabla }p+\frac {1}{\rho }\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\tau }, \end{equation}

where

(2.3) bold italic tau equals mu left parenthesis ModifyingAbove gamma With dot right parenthesis left parenthesis bold nabla bold italic u plus left parenthesis bold nabla bold italic u right parenthesis Superscript upper T Baseline right parenthesis comma \begin{equation} \boldsymbol{\tau } = \mu (\dot \gamma )(\boldsymbol{\nabla }\boldsymbol{u} + (\boldsymbol{\nabla }\boldsymbol{u})^T), \end{equation}

where bold italic u $\boldsymbol{u}$ , t and p are the velocity, time and pressure of the flow; rho $\rho$ and mu left parenthesis ModifyingAbove gamma With dot right parenthesis $\mu (\dot \gamma )$ are the density and the dynamic viscosity of the fluid, respectively. It should be noted that bold nabla bold dot left bracket mu left parenthesis ModifyingAbove gamma With dot right parenthesis bold nabla bold italic u right bracket $\boldsymbol{\nabla }\boldsymbol{\cdot }[\mu (\dot \gamma )\boldsymbol{\nabla }\boldsymbol{u}]$ cannot be simplified to mu left parenthesis ModifyingAbove gamma With dot right parenthesis bold nabla squared bold italic u $\mu (\dot \gamma )\boldsymbol{\nabla} ^2\boldsymbol{u}$ unless the viscosity is spatially constant. In the present study, the simulations are performed using the standard finite-volume solver in OpenFOAM, where the viscous term is discretised directly in the conservative stress-divergence form based on (2.3), and the viscosity is updated at each cell face according to the local velocity gradients using the Carreau model. Thus, the spatial viscosity variations are rigorously accounted for. The equations are non-dimensionalised using the gap width d $d$ as the length scale, the inner cylinder tangential velocity u Subscript i Baseline equals r Subscript i Baseline omega Subscript i $u_i = r_i\omega _i$ as the velocity scale, the large-eddy turnover time t Superscript asterisk $t^{\ast }$ = d divided by u Subscript i $d/u_i$ , and the inertial pressure scale rho u Subscript i Superscript 2 $\rho u_i^2$ . We also define the non-dimensional radius r Superscript asterisk $r^{\ast }$ to be r Superscript asterisk $r^{\ast }$ = ( r minus r Subscript i $r- r_{i}$ )/d.

In TC flow, the Reynolds number is defined as

(2.4) upper R e equals StartFraction rho omega Subscript i Baseline r Subscript i Baseline d Over mu Subscript r e f Baseline EndFraction comma \begin{equation} Re=\frac {\rho \omega _{i}r_{i}d}{\mu _{ref}}, \end{equation}

where the reference viscosity mu Subscript r e f $\mu _{ref}$ is evaluated at the nominal shear rate ModifyingAbove gamma With dot Subscript n o m Baseline equals omega Subscript i Baseline r Subscript i Baseline divided by d $\dot \gamma _{nom} = \omega _i r_i / d$ , following the convention established for shear-thinning TC flow (Cagney & Balabani Reference Cagney and Balabani2019). The shear-rate-dependent viscosity is described by the Carreau model:

(2.5) mu left parenthesis ModifyingAbove gamma With dot right parenthesis equals mu Subscript normal infinity Baseline plus left parenthesis mu 0 minus mu Subscript normal infinity Baseline right parenthesis left bracket 1 plus left parenthesis lamda Subscript c Baseline ModifyingAbove gamma With dot right parenthesis squared right bracket Superscript left parenthesis n minus 1 right parenthesis divided by 2 Baseline comma \begin{equation} \mu (\dot \gamma )=\mu _{\infty }+(\mu _{0}-\mu _{\infty })[1+(\lambda _{c}\dot \gamma )^2]^{(n-1)/2}, \end{equation}

where mu 0 $\mu _{0}$ and mu Subscript normal infinity $\mu _{\infty }$ are the zero- and infinite-shear-rate viscosities, lamda Subscript c $\lambda _{c}$ is the characteristic relaxation time and n $n$ is the shear-thinning index ( n equals 1 $n=1$ for a Newtonian fluid).

Equations (2.1) and (2.2) are solved using the open-source finite-volume platform OpenFOAM. A second-order-accurate, collocated grid arrangement is employed. Temporal discretisation uses the second-order backward scheme, while the convective terms are discretised with a second-order total variation diminishing (Van Leer) scheme. The computational domain is discretised using a structured mesh with upper N Subscript r Baseline equals 64 $N_r = 64$ , upper N Subscript theta Baseline equals 256 $N_\theta = 256$ and upper N Subscript z Baseline equals 128 $N_z = 128$ cells in the radial, azimuthal and axial directions, respectively. The mesh is refined in the near-wall regions to adequately resolve the steep velocity gradients and the thin boundary layers characteristic of shear-thinning TC flow. The minimum radial grid spacing near the walls is approximately 0.008 d $0.008d$ , while the maximum spacing in the gap centre is 0.03 d $0.03d$ . The azimuthal and axial directions employ uniform meshes with spacings of approximately 0.13 d $0.13d$ and 0.05 d $0.05d$ , respectively. A fixed time step is chosen such that the Courant–Friedrichs–Lewy (CFL) number remains below 0.7 throughout all simulations. As a global measure of numerical convergence, we monitor the torque conservation for every Reynolds number. In all cases, the relative difference between the torque acting on the inner cylinder and that acting on the outer cylinder remains below 1 % after the flow reaches a statistically steady state. The simulations are advanced until the time-averaged torques on the inner and outer cylinders differ by less than 1 %, and statistics are collected over at least 100 large-eddy turnover times ( d divided by u Subscript i $d/u_i$ ) after this quasi-steady state is reached (Dong Reference Dong2007; Brauckmann & Eckhardt Reference Brauckmann and Eckhardt2013). Velocity and pressure fields computed in Cartesian coordinates are transformed to cylindrical coordinates for post-processing.

The solver has been validated for Newtonian TC flow in our previous work (Xu et al. Reference Xu, Su, Lan, Zhao, He, Sun and Wang2022, Reference Xu, Zhao, Sun, He and Wang2023) and further verified for shear-thinning fluids by comparing DNS-based effective viscosities with the Carreau model at low italic Re $\textit{Re}$ (Re = 40), where the flow is purely azimuthal. The agreement confirms correct implementation of the rheology (figure 1 b–d). To determine the critical Reynolds numbers italic Re Subscript c comma 1 $\textit{Re}_{c,1}$ and italic Re Subscript c comma 2 $\textit{Re}_{c,2}$ , the simulations are initiated from a low Reynolds number italic Re equals 40 $\textit{Re}=40$ , where the flow remains in the CCF state. The Reynolds number was increased gradually in coarse steps to identify the approximate onset of the TVF state. Once the TVF state was established, the Reynolds number was reduced back into the stable CCF regime to obtain a robust initial condition. Subsequently, italic Re $\textit{Re}$ was increased again with a much finer increment of normal upper Delta upper R e equals 1 $\Delta Re = 1$ to accurately determine the critical Reynolds number for the primary instability italic Re Subscript c comma 1 $\textit{Re}_{c,1}$ . The same procedure was applied for the secondary transition from TVF to WVF. The onset of TVF is identified by steady axisymmetric vortex structures, while WVF is detected by the emergence of travelling azimuthal waves, as verified by the time series of the axial velocity at multiple monitor points. For each fluid, simulations are conducted up to italic Re almost equals 1000 $\textit{Re} \approx 1000$ , sufficient to capture both instabilities. Results at higher italic Re $\textit{Re}$ are provided in the supplementary material. Three shear-thinning indices ( n equals 0.52 comma 0.62 comma 0.68 $n=0.52, 0.62, 0.68$ ) are considered using parameters listed in table 1, taken from Cagney & Balabani (Reference Cagney and Balabani2019). The critical Reynolds numbers obtained from the present DNS for each fluid are also included in table 1.

Table 1.

Carreau model parameters for the three shear-thinning fluids investigated (reproduced from Cagney & Balabani (Reference Cagney and Balabani2019)) together with the critical Reynolds numbers italic Re Subscript c comma 1 $\textit{Re}_{c,1}$ and italic Re Subscript c comma 2 $\textit{Re}_{c,2}$ obtained from the present DNS. Here mu 0 $\mu _0$ and mu Subscript normal infinity $\mu _{\infty }$ are the zero- and infinite-shear-rate viscosities, lamda Subscript c $\lambda _c$ is the characteristic relaxation time and n $n$ is the shear-thinning index.

A table comparing parameters of three shear-thinning fluids.
Table 1. Long description

A table comparing parameters of three shear-thinning fluids. The table has three rows and six columns. Column headers are Fluid, μ0 (mPa s), μ∞ (mPa s), λc (s), n, Rec,1, and Rec,2. Row 1: Fluid, 1; μ0 (mPa s), 85.32; μ∞ (mPa s), 0.935; λc (s), 1.29; n, 0.52; Rec,1, 94; Rec,2, 142. Row 2: Fluid, 2; μ0 (mPa s), 23.11; μ∞ (mPa s), 0.976; λc (s), 0.89; n, 0.62; Rec,1, 97; Rec,2, 134. Row 3: Fluid, 3; μ0 (mPa s), 13.64; μ∞ (mPa s), 1.732; λc (s), 0.50; n, 0.68; Rec,1, 101; Rec,2, 115.

3. Results and discussion

Figure 2(a) summarises the critical Reynolds numbers for the onset of TVF ( italic Re Subscript c comma 1 $\textit{Re}_{c,1}$ ) and WVF ( italic Re Subscript c comma 2 $\textit{Re}_{c,2}$ ) as functions of the shear-thinning index n $n$ , obtained from the present DNS and compared with the experimental measurements of Cagney & Balabani (Reference Cagney and Balabani2019). The phase diagram in figure 2(a) delineates the stability boundaries separating the CCF, TVF and WVF regimes (Andereck et al. Reference Andereck, Liu and Swinney1986). Figure 2(b) displays representative contours of the azimuthal vorticity in the meridional ( r $r$ – z $z$ ) plane for the TVF and WVF regimes, illustrating the structural changes in the flow associated with each instability.

Figure 2.

(a) Critical Reynolds numbers for the onset of TVF ( italic Re Subscript c comma 1 $\textit{Re}_{c,1}$ ) and WVF ( italic Re Subscript c comma 2 $\textit{Re}_{c,2}$ ) as functions of the shear-thinning index n $n$ . The DNS results are shown as open symbols connected by lines, while the open symbols denote the experimental data of Cagney & Balabani (Reference Cagney and Balabani2019). The regimes CCF, TVF and WVF are denoted by text labels in the diagram. (b) Instantaneous contours of the azimuthal vorticity (normalised by the rotating angular velocity of inner cylinder) in the meridional ( r minus z $r{-}z$ ) plane at two azimuthal positions ( theta equals 0 $\theta =0$ and theta equals pi $\theta =\pi$ ) for TVF ( italic Re Subscript c comma 1 Baseline equals 94 $\textit{Re}_{c,1}=94$ ) and WVF ( italic Re Subscript c comma 2 Baseline equals 142 $\textit{Re}_{c,2}=142$ ) at n equals 0.52 $n=0.52$ . The steady axisymmetric TVF structure and the symmetry breaking in WVF are clearly visible. Three triangular markers at z equals upper L divided by 4 $z = L/4$ , upper L divided by 2 $L/2$ and 3 upper L divided by 4 $3L/4$ are superimposed on the WVF contour at theta equals 0 $\theta =0$ , indicating the locations of the three monitor points used in (c); all three probes are at the mid-gap r Superscript asterisk Baseline equals 0.5 $r^*=0.5$ . (c) Time evolution of the axial velocity u Subscript z Baseline divided by u Subscript i $u_z/u_i$ at the mid-gap ( r Superscript asterisk Baseline equals 0.5 $r^*=0.5$ ) for n equals 0.52 $n=0.52$ , recorded at three axial monitor points ( z equals upper L divided by 4 comma upper L divided by 2 comma 3 upper L divided by 4 $z=L/4, L/2, 3L/4$ ). Dashed lines denote the steady state at the primary instability threshold ( italic Re Subscript c comma 1 Baseline equals 94 $\textit{Re}_{c,1}=94$ ), while solid lines illustrate the quasi-periodic oscillations at the secondary instability threshold ( italic Re Subscript c comma 2 Baseline equals 142 $\textit{Re}_{c,2}=142$ ). The overlapping traces confirm the existence of an axially coherent, global travelling wave in the WVF regime. (d) Circumferential profile of u Subscript z Baseline divided by u Subscript i $u_z/u_i$ at the mid-gap ( r Superscript asterisk Baseline equals 0.5 $r^*=0.5$ ) for n equals 0.52 $n=0.52$ (the corresponding profiles for other n $n$ values are provided in the Appendix C). The profile exhibits three complete wave periods, corresponding to a dominant azimuthal mode m = 3.

Multiple graphs depict critical Reynolds numbers, vorticity contours, axial velocity, and circumferential profiles in Taylor-Couette flow.

For the primary instability (CCF right arrow $\rightarrow$ TVF), the DNS results reveal a clear monotonic trend: as the fluid becomes more shear-thinning (decreasing n $n$ ), the critical Reynolds number italic Re Subscript c comma 1 $\textit{Re}_{c,1}$ decreases from 101 $101$ at n equals 0.68 $n=0.68$ to 94 $94$ at n equals 0.52 $n=0.52$ . This destabilising effect is consistent with the physical picture that enhanced shear-thinning reduces the effective viscosity in the high-shear region adjacent to the inner cylinder, thereby amplifying the local centrifugal forces that drive the Taylor instability (Alibenyahia et al. Reference Alibenyahia, Lemaitre, Nouar and Ait-Messaoudene2012; Khali et al. Reference Khali, Nebbali and Bouhadef2013; Cagney & Balabani Reference Cagney and Balabani2019). In contrast, for the secondary instability (TVF right arrow $\rightarrow$ WVF), shear-thinning stabilises the wavy mode: italic Re Subscript c comma 2 $\textit{Re}_{c,2}$ increase from 115 $115$ at n equals 0.68 $n=0.68$ to 142 $142$ at n equals 0.52 $n=0.52$ . This opposing effect implies that the radial viscosity variation, while promoting centrifugal destabilisation, simultaneously suppresses the oscillatory shear-layer instability associated with inflectional velocity profiles between adjacent vortices. The experimental data of Cagney & Balabani (Reference Cagney and Balabani2019), however, exhibit a distinctly different, non-monotonic behaviour for both instabilities. It should be noted that the experimental values of italic Re Subscript c comma 2 $\textit{Re}_{c,2}$ are consistently far larger than their DNS counterparts. These discrepancies are direct manifestations of endwall effects inherent in finite-length experimental setups (Ji et al. Reference Ji, Burin, Schartman and Goodman2006), which are absent in the present axially periodic DNS (Cole Reference Cole1976; Czarny et al. Reference Czarny, Serre, Bontoux and Lueptow2003). While no-slip endwalls impose axial wavelength constraints, the axially periodic boundary conditions used here also restrict the admissible axial wavelengths. The key distinction is that the rigid endwalls mechanically suppress the amplitude of the travelling azimuthal waves that characterise the WVF regime (Cole Reference Cole1976; Jones Reference Jones1985), thereby artificially elevating the apparent italic Re Subscript c comma 2 $\textit{Re}_{c,2}$ . Indeed, the significant spatial variability of wave amplitudes observed in the experiments of Cagney & Balabani (Reference Cagney and Balabani2019) for shear-thinning fluids is itself a signature of these endwall-induced constraints. In contrast, the periodic boundary condition allows the travelling wave to propagate continuously across the axial boundaries without endwall-induced damping.

The vorticity contours in figure 2(b) provide structural evidence for the two distinct flow states. In the TVF regime, the azimuthal vorticity appears as three regularly spaced, alternating-sign vortex pairs along the axial direction, consistent with the imposed axial periodicity of upper Gamma equals 2 pi $\varGamma = 2\pi$ . The steady, axisymmetric nature of these Taylor vortices is confirmed by the constant-time signals at the three monitor points shown in figure 2(c). The spatial structure emerges spontaneously from the nonlinear evolution within the periodic domain rather than from any prescribed forcing. Even if the axial length upper L $L$ accommodates an integer number of vortex pairs, the vortex boundaries do not necessarily coincide with the periodic boundaries. In the present simulations, the observed alignment is a natural consequence of the specific mode selected during the nonlinear evolution, rather than a prescribed condition or an artefact of the chosen axial length. Upon transition to WVF, this symmetry is broken: the vorticity contours become distorted and lose their axial regularity, indicative of the travelling azimuthal waves that characterise the wavy instability (Wereley & Lueptow Reference Wereley and Lueptow1998). Overall, shear-thinning has a destabilising effect on the primary instability but a significantly stronger stabilising effect on the secondary instability.

The rheological parameters in table 1 are extracted from the experimental fluids, where mu 0 $\mu _0$ , mu Subscript normal infinity $\mu _\infty$ , lamda Subscript c $\lambda _c$ and n $n$ vary simultaneously. To isolate the individual influence of the shear-thinning index n $n$ from the viscosity levels themselves, we performed a controlled set of simulations where mu 0 $\mu _0$ , mu Subscript normal infinity $\mu _\infty$ and lamda Subscript c $\lambda _c$ were fixed while n $n$ was varied. The results confirm that the opposing trends are robust regardless of whether the other Carreau parameters are varied or fixed: a decreasing n $n$ consistently lowers italic Re Subscript c comma 1 $\textit{Re}_{c,1}$ and elevates italic Re Subscript c comma 2 $\textit{Re}_{c,2}$ . The quantitative differences observed between the fixed-parameter simulations and the original experimental fluids, such as the variation in italic Re Subscript c comma 2 $\textit{Re}_{c,2}$ for n equals 0.68 $n=0.68$ , arise primarily from the substantial differences in the absolute viscosity levels. The parameters and numerical results of these controlled simulations are provided in Appendix A. We also validated our numerical framework against the linear stability benchmarks of Alibenyahia et al. (Reference Alibenyahia, Lemaitre, Nouar and Ait-Messaoudene2012) for power-law fluids. Although their results were obtained for a slightly different radius ratio and for a power-law rather than Carreau rheology, both sets exhibit the same physical trend: italic Re Subscript c comma 1 $\textit{Re}_{c,1}$ decreases monotonically with decreasing n when the Reynolds number is defined using the nominal shear-rate viscosity. This consistency confirms the reliability of our DNS methodology for capturing the primary instability. The parameters and numerical results of this validation are provided in Appendix B.

Furthermore, the use of axially periodic boundary conditions imposes a constraint on the admissible axial wavelengths. To evaluate the sensitivity of our results to the specific choice of the axial periodicity length, we performed additional simulations for the strongest shear-thinning fluid ( n equals 0.52 $n=0.52$ ) using two alternative periodic lengths ( upper L equals 6 d $L=6d$ and upper L equals 6.6 d $L=6.6d$ ). The resulting critical Reynolds numbers are as follows: for upper L equals 6 d $L=6d$ , italic Re Subscript c comma 1 Baseline equals 94 $\textit{Re}_{c,1}=94$ , italic Re Subscript c comma 2 Baseline equals 142 $\textit{Re}_{c,2}=142$ ; for upper L equals 6.6 d $L=6.6d$ , italic Re Subscript c comma 1 Baseline equals 94 $\textit{Re}_{c,1}=94$ , italic Re Subscript c comma 2 Baseline equals 141 $\textit{Re}_{c,2}=141$ . These results demonstrate that the primary instability threshold ( italic Re Subscript c comma 1 $\textit{Re}_{c,1}$ ) is independent of the chosen periodic length, while the secondary instability threshold ( italic Re Subscript c comma 2 $\textit{Re}_{c,2}$ ) varies by only 1 Reynolds number unit. These findings confirm that the conclusions presented in this study are not artefacts of a specific axial wavelength selection, in line with the established understanding of wavelength constraints in periodic TC configurations (Cole Reference Cole1976; Jones Reference Jones1985; Topayev et al. Reference Topayev, Nouar and Dusek2022).

To further characterise the unsteady and three-dimensional features of the WVF regime, we present quantitative analyses in figures 2(c) and 2(d). Figure 2(c) shows the time evolution of the axial velocity u Subscript z $u_z$ at the mid-gap ( r Superscript asterisk Baseline equals 0.5 $r^*=0.5$ ) for n equals 0.52 $n=0.52$ , recorded at three distinct axial monitor points ( z equals upper L divided by 4 comma upper L divided by 2 comma 3 upper L divided by 4 $z=L/4, L/2, 3L/4$ ). In the TVF regime, u Subscript z $u_z$ remains almost constant over time at all three positions, consistent with the steady, axisymmetric nature of Taylor vortices. As the flow transitions to WVF, the axial velocity at all three monitor points exhibits clear quasi-periodic oscillations. These quasi-periodic oscillations are captured at italic Re equals italic Re Subscript c comma 2 Baseline $\textit{Re}=\textit{Re}_{c,2}$ , the onset of the secondary instability. At this critical condition, the WVF perturbation has already emerged with a finite amplitude, as is typical for supercritical bifurcations in TC flow. Notably, the temporal signals from the three spatially separated axial points almost perfectly overlap, the maximum absolute deviation among the three signals is approximately 5 % of the signal amplitude, which is far smaller than the dynamic range of the oscillations and confirms that the three signals are indeed highly coherent. Although the three axial positions correspond to different local phases of the steady Taylor-vortex structure, the WVF perturbation itself is a global travelling wave that occupies the entire periodic domain (Wereley & Lueptow Reference Wereley and Lueptow1998). Because our simulations employ axially periodic boundary conditions, which eliminate the endwall-induced spatial damping that is inevitably present in experimental setups, the travelling wave is free to propagate without axial attenuation. Consequently, its axial envelope remains approximately uniform over the periodic length, yielding nearly identical amplitudes at all three monitoring points regardless of their different positions within the underlying vortex structure.

This interpretation is further supported by the observation that the dominant azimuthal mode m = 3 remains unchanged with n (figure 2(d) and Appendix C), confirming that the secondary instability is governed by a large-scale global mode rather than a locally confined fluctuation. Regarding the azimuthal structure, figure 2(d) presents the circumferential profiles of u Subscript z $u_z$ for the n = 0.52 case at the mid-gap ( r Superscript asterisk Baseline equals 0.5 $r^*=0.5$ ) (the corresponding profiles for n = 0.62 and 0.68 are provided in Appendix C). The profiles consistently exhibit three complete wave periods for all the investigated n $n$ values, confirming a dominant azimuthal mode of m equals 3 $m=3$ and the three-dimensional travelling-wave character of the WVF regime (Jones Reference Jones1985; Andereck et al. Reference Andereck, Liu and Swinney1986). More importantly, the fact that m equals 3 $m=3$ remains constant across the different n $n$ values provides critical physical insight into the stabilising mechanism. It indicates that shear-thinning, while significantly postponing the onset of WVF (increasing italic Re Subscript c comma 2 $\textit{Re}_{c,2}$ ), does not alter the preferred azimuthal mode selection. This implies that the stabilisation of the secondary instability originates from the enhanced viscous damping in the bulk gap that suppresses the growth of the travelling waves, rather than from a change in the spatial wavelength or mode structure.

It is also worth clarifying that the apparent shrinking of the Taylor-vortex stability range with increasing shear-thinning reported by Topayev et al. (Reference Topayev, Nouar, Bernardin, Neveu and Bahrani2019) is not in conflict with our present findings. This discrepancy originates from the different reference viscosity used in the definition of the Reynolds number. Following the convention of experimental studies (Cagney & Balabani Reference Cagney and Balabani2019), we define the Reynolds number using the viscosity evaluated at the nominal shear rate. In contrast, Topayev et al. (Reference Topayev, Nouar, Bernardin, Neveu and Bahrani2019) employed the zero-shear-rate viscosity mu 0 $\mu _0$ as the reference. When our critical Reynolds numbers are recalculated using mu 0 $\mu _0$ as the reference scale, the resulting values reproduce the qualitative trend reported by Topayev et al. (Reference Topayev, Nouar, Bernardin, Neveu and Bahrani2019), namely a shrinking TVF stability interval with decreasing n $n$ . This demonstrates that the apparent contradiction is solely a consequence of the choice of the viscosity scale, rather than a conflict in the underlying physical mechanisms.

To clarify the physical mechanism underlying the contrasting effects of shear-thinning on italic Re Subscript c comma 1 $\textit{Re}_{c,1}$ and italic Re Subscript c comma 2 $\textit{Re}_{c,2}$ , we examine the radial profiles of the time- and spatially-averaged dimensionless viscosity and its gradient at the respective critical Reynolds numbers. Figure 3 presents these profiles for all three values of n, evaluated at italic Re equals italic Re Subscript c comma 1 Baseline $\textit{Re} = \textit{Re}_{c,1}$ (a,b) and italic Re equals italic Re Subscript c comma 2 Baseline $\textit{Re} = \textit{Re}_{c,2}$ (c,d). The viscosity is normalised by the viscosity mu Subscript 0 $\mu {_0}$ at the zero shear rate, and the radial coordinate is r Superscript asterisk $r^{\ast }$ = ( r minus r Subscript i $r-r_{i}$ )/d, with r Superscript asterisk $r^{\ast }$ = 0 and r Superscript asterisk $r^{\ast }$ = 1 corresponding to the inner and outer cylinder walls, respectively.

Figure 3.

Radial profiles of the time-, azimuthally- and axially-averaged dimensionless viscosity mu left parenthesis ModifyingAbove gamma With dot right parenthesis divided by mu 0 $\mu (\dot \gamma )/\mu _0$ (a,c) and its radial gradient d mu left parenthesis ModifyingAbove gamma With dot right parenthesis slash d r $\text{d}\mu (\dot \gamma )/\text{d}r$ (b,d) for three shear-thinning fluids with n equals 0.52 $n = 0.52$ , 0.62 $0.62$ and 0.68 $0.68$ . Panels (a,b) correspond to the primary instability threshold ( italic Re equals italic Re Subscript c comma 1 Baseline $\textit{Re} = \textit{Re}_{c,1}$ ); panels (c,d) correspond to the secondary instability threshold ( italic Re equals italic Re Subscript c comma 2 Baseline $\textit{Re} = \textit{Re}_{c,2}$ ). The dashed lines in (a,c) represent the average values of dimensionless viscosity at different n, and the dashed lines in (b,d) represent the radial gradient where d mu left parenthesis ModifyingAbove gamma With dot right parenthesis slash d r equals 0 $\text{d}\mu (\dot \gamma )/\text{d}r = 0$ .

Multiple graphs depict the radial profiles of dimensionless viscosity and its radial gradient for shear-thinning fluids at different instability thresholds.
Figure 3. Long description

The image contains four line graphs. Panel A shows the radial profiles of the time-, azimuthally- and axially-averaged dimensionless viscosity for three shear-thinning fluids with different values of n at the primary instability threshold. The x-axis is labeled r* and the y-axis is labeled μ(γ)/μ0. Panel B shows the radial gradient of the dimensionless viscosity for the same fluids and threshold, with the x-axis labeled r* and the y-axis labeled dμ(γ)/dr. Panel C presents similar data to Panel A but at the secondary instability threshold. Panel D presents similar data to Panel B but at the secondary instability threshold. Each graph includes dashed lines representing average values or gradients.

At the primary instability threshold ( italic Re equals italic Re Subscript c comma 1 Baseline $\textit{Re} = \textit{Re}_{c,1}$ , figure 3 a,b), the viscosity profiles exhibit a pronounced variation across the gap for all shear-thinning cases. The vertical axes in figure 3 represent the normalised viscosity mu divided by mu 0 $\mu /\mu _0$ . As n $n$ decreases, this normalised near-wall viscosity is substantially reduced, reflecting the stronger local shear-thinning effect described by the Carreau model, although the absolute viscosity mu 0 $\mu _0$ for smaller n $n$ remains higher (see table 1). This reduction is accompanied by a sharpening of the viscosity gradient d mu slash d r $\text{d} \mu / \text{d} r$ in the inner region ( r Superscript asterisk Baseline less than or equals 0.3 $r^{\ast } \le 0.3$ ), which reaches a maximum at approximately r Superscript asterisk Baseline almost equals 0.25 $r^{\ast } \approx 0.25$ in the inner portion of the gap and then decays towards the gap centre. The gradient values at the wall ( r Superscript asterisk Baseline equals 0 $r^*=0$ ) and at approximately r Superscript asterisk Baseline almost equals 0.45 $r^* \approx 0.45$ are comparable, indicating a non-monotonic gradient distribution. The enhanced positive gradient near the rotating inner cylinder effectively lubricates the flow, diminishing the local viscous resistance and thereby facilitating the centrifugal instability that triggers the formation of Taylor vortices (Ashrafi Reference Ashrafi2011; Alibenyahia et al. Reference Alibenyahia, Lemaitre, Nouar and Ait-Messaoudene2012).The subsequent negative gradient, although weaker, may further influence the radial distribution of angular momentum transport. Nevertheless, the overall effect accounts for the monotonic decrease of italic Re Subscript c comma 1 $\textit{Re}_{c,1}$ with decreasing n $n$ observed in figure 2(a).

At the secondary instability threshold ( italic Re equals italic Re Subscript c comma 2 Baseline $\textit{Re} = \textit{Re}_{c,2}$ , figure 3 c,d), the overall viscosity levels are lower than at italic Re Subscript c comma 1 $\textit{Re}_{c,1}$ due to the higher rotation rates, yet the qualitative trends with n $n$ persist. Notably, the viscosity gradient profiles at italic Re Subscript c comma 2 $\textit{Re}_{c,2}$ display a more pronounced peak near the inner wall and a steeper decay in the outer region compared to those at italic Re Subscript c comma 1 $\textit{Re}_{c,1}$ . This intensified near-wall viscosity gradient has a dual effect. While it continues to promote centrifugal destabilisation, it also acts to confine the high-shear region to a thin boundary layer adjacent to the inner cylinder, leaving the bulk of the gap with a relatively elevated and more uniform viscosity. The wavy vortex instability, by contrast, is primarily associated with the destabilisation of the azimuthal velocity shear layers, partial differential u Subscript theta divided by partial differential z $\partial u_\theta / \partial z$ , generated by the axial modulation of the Taylor vortex structure (Martinand et al. Reference Martinand, Serre and Lueptow2014; Dessup et al. Reference Dessup, Tuckerman, Wesfreid, Barkley and Willis2018). Shear-thinning stabilises the WVF by increasing the viscous damping of these shear layers, while the accompanying weakening of the radial jets further reduces the nonlinear coupling between adjacent vortices. Consequently, italic Re Subscript c comma 2 $\textit{Re}_{c,2}$ increases as n $n$ decreases, yielding the opposing trend evident in figure 2(a).

The viscosity gradient profiles also reveal a subtle difference between the three fluids. For the most strongly shear-thinning case ( n equals 0.52 $n=0.52$ ), the gradient at italic Re Subscript c comma 2 $\textit{Re}_{c,2}$ reaches its maximum closer to the inner wall and decays more abruptly than for n equals 0.62 $n=0.62$ and 0.68 $0.68$ . This correlates with the particularly pronounced stabilisation of the wavy mode observed for n equals 0.52 $n=0.52$ in figure 2(a). The radial confinement of the low-viscosity region thus appears to be a key factor in the suppression of the secondary instability.

In summary, the radial viscosity variation induced by shear-thinning exerts a dual influence: it promotes the primary centrifugal instability by reducing the effective viscosity near the driving wall, while it stabilises the secondary wavy instability by increasing the viscous damping in the bulk of the gap and restricting the low-viscosity region to a thin boundary layer. This differential effect constitutes the fundamental origin of the opposing effects reported in the present study. In the following section, we further explore this mechanism by analysing the full two-dimensional flow fields for the n equals 0.52 $n=0.52$ case at both critical Reynolds numbers.

To gain deeper insight into the mechanisms underlying the opposing stability trends, we now examine the full two-dimensional flow fields in the meridional ( r minus z $r{-}z$ ) plane for the most strongly shear-thinning fluid ( n equals 0.52 $n=0.52$ ). Figure 4 presents contours of the instantaneous dimensionless shear rate lamda Subscript c Baseline ModifyingAbove gamma With dot $\lambda {_c}\dot \gamma$ , viscosity mu 0 $\mu _0$ , and azimuthal velocity u Subscript theta Baseline divided by u Subscript i $u_\theta /u_i$ superposed velocity vectors (the radial and axial components) at the critical Reynolds numbers for the primary instability (figure 4 a, italic Re Subscript c comma 1 Baseline equals 94 $\textit{Re}_{c,1} = 94$ ) and the secondary instability (figure 4 b, italic Re Subscript c comma 2 Baseline equals 142 $\textit{Re}_{c,2} = 142$ ). The local shear rate ModifyingAbove gamma With dot $\dot {\gamma }$ is calculated from the second invariant of the strain-rate tensor, ModifyingAbove gamma With dot equals StartRoot 2 bold italic upper D colon bold italic upper D EndRoot $\dot {\gamma } = \sqrt {2\boldsymbol{D}:\boldsymbol{D}}$ , where bold italic upper D equals left parenthesis bold nabla bold italic u plus left parenthesis bold nabla bold italic u right parenthesis Superscript upper T Baseline right parenthesis divided by 2 $\boldsymbol{D} = (\boldsymbol{\nabla }\boldsymbol{u} + (\boldsymbol{\nabla }\boldsymbol{u})^T)/2$ is the strain-rate tensor.

Figure 4.

Contours of the instantaneous (left) shear rate lamda Subscript c Baseline ModifyingAbove gamma With dot $\lambda _c\dot \gamma$ , (middle) dimensionless viscosity mu divided by mu 0 $\mu /\mu _0$ and (right) azimuthal velocity u Subscript theta Baseline divided by u Subscript i $u_\theta /u_i$ superposed velocity vectors (the radial and axial components) in the meridional ( r $r$ – z $z$ ) plane for the most shear-thinning fluid ( n equals 0.52 $n = 0.52$ ): (a) at the primary instability threshold ( italic Re Subscript c comma 1 Baseline equals 94 $\textit{Re}_{c,1} = 94$ ) and (b) at the secondary instability threshold ( italic Re Subscript c comma 2 Baseline equals 142 $\textit{Re}_{c,2} = 142$ ).

Heat maps showing fluid dynamics data with velocity vectors.

At italic Re Subscript c comma 1 $\textit{Re}_{c,1}$ (figure 4 a), the flow exhibits the characteristic structure of steady, axisymmetric TVF. Three counter-rotating vortex pairs are clearly discernible, regularly spaced along the axial direction with a wavelength corresponding to the imposed periodicity upper Gamma equals 2 pi $\varGamma = 2\pi$ . The shear-rate field reveals concentrated bands of high ModifyingAbove gamma With dot $\dot \gamma$ adjacent to the inner cylinder wall and along the radial outflow jets separating adjacent vortices, where fluid is ejected from the inner boundary towards the outer cylinder. Correspondingly, the viscosity contours display pronounced spatial variations: the viscosity is lowest in the high-shear near-wall region and along the outflow jets, while significantly higher values occupy the vortex cores and the region near the stationary outer cylinder, where the local shear rate is diminished (Topayev et al. Reference Topayev, Nouar, Bernardin, Neveu and Bahrani2019). This spatial heterogeneity in viscosity is a direct consequence of the shear-thinning rheology and underpins the radial variation analysed in figure 3. The azimuthal velocity field shows the expected monotonic decay from the inner to the outer cylinder, with weak undulations induced by the toroidal vortices. The velocity vectors reveal a mild asymmetry: the outflow jet accelerates slightly as it moves outwards and the inflow jet decelerates gently inwards.

At italic Re Subscript c comma 2 $\textit{Re}_{c,2}$ (figure 4 b), the flow has transitioned to WVF, and the structural symmetry evident at italic Re Subscript c comma 1 $\textit{Re}_{c,1}$ is broken. The vortex boundaries are no longer straight axial lines; instead, they exhibit undulations characteristic of the travelling azimuthal waves that define the WVF regime (Andereck et al. Reference Andereck, Liu and Swinney1986; Wereley & Lueptow Reference Wereley and Lueptow1998). The shear-rate contours become distorted, with the outflow boundaries becoming inclined and tilted relative to the axial direction. The low-viscosity regions exhibit a similar tilted orientation to the outflow boundaries, and their shapes are similar across the flow field. Importantly, the azimuthal velocity field reveals a weakening of the inward jets relative to the outward jets, a feature that has been identified as a hallmark of shear-thinning TC flow (Escudier et al. Reference Escudier, Gouldson and Jones1995; Cagney & Balabani Reference Cagney and Balabani2019; Topayev et al. Reference Topayev, Nouar, Bernardin, Neveu and Bahrani2019). The velocity vectors in figure 4(b) show a much stronger asymmetry: the outflow jet accelerates rapidly as it moves outwards, while the inflow jet decelerates sharply when approaching the inner wall. This progression is consistent with the steeper viscosity gradient at the higher Reynolds number (figure 3 d): the low-viscosity zone adjacent to the inner wall facilitates outward acceleration, but the same zone fails to sustain the inward return flow, thereby weakening the inter-vortex coupling required for wavy instability. The reduced strength of the inward jets diminishes the axial modulation of the Taylor vortex structure, thereby weakening the azimuthal velocity shear layers, partial differential u Subscript theta divided by partial differential z $\partial u_\theta / \partial z$ , which are the primary drivers of the wavy instability (Martinand et al. Reference Martinand, Serre and Lueptow2014; Dessup et al. Reference Dessup, Tuckerman, Wesfreid, Barkley and Willis2018). In other words, the inward jets are essential for sustaining the strong axial gradients of the azimuthal velocity that characterise the WVF regime (Jones Reference Jones1985). The diminished nonlinear coupling between adjacent vortices, combined with the elevated bulk viscosity noted in figure 3(c), thus provides a consistent explanation for the delayed onset of WVF with enhanced shear-thinning.

The contrasting flow structures at italic Re Subscript c comma 1 $\textit{Re}_{c,1}$ and italic Re Subscript c comma 2 $\textit{Re}_{c,2}$ thus corroborate the differential effect of shear-thinning on the two instabilities. At the primary threshold, the near-wall viscosity reduction dominates, promoting the centrifugal overturning that generates Taylor vortices. At the secondary threshold, however, the stabilising influences, elevated viscosity in the vortex cores, weakened inward jets and confinement of the low-viscosity region to a thin boundary layer, collectively outweigh any residual destabilising effects, leading to the observed increase in italic Re Subscript c comma 2 $\textit{Re}_{c,2}$ with decreasing n $n$ .

4. Conclusions

DNS with axially periodic boundary conditions reveal that shear-thinning exerts opposing effects on successive instabilities in narrow-gap TC flow: enhanced shear-thinning (decreasing n $n$ ) destabilises the primary instability ( italic Re Subscript c comma 1 $\textit{Re}_{c,1}$ decreases) but strongly stabilises the secondary instability ( italic Re Subscript c comma 2 $\textit{Re}_{c,2}$ increases). This behaviour originates from a competition between two opposing consequences of radial viscosity variation: the reduced viscosity near the inner wall promotes centrifugal instability (destabilising), while the elevated viscosity in the bulk gap damps inter-vortex oscillations (stabilising). Comparison with experiments has demonstrated that previously reported non-monotonic trends arise from endwall-induced extrinsic constraints, such as wavelength confinement and amplitude suppression of the azimuthal waves, rather than intrinsic rheological effects. These findings establish the intrinsic bifurcation structure of shear-thinning TC flow in the narrow-gap configuration ( eta equals 0.83 $\eta =0.83$ ) and clarify the distinct roles of rheology and conffnement in successive transitions.

Supplementary material

Supplementary material are available at https://doi.org/10.1017/jfm.2026.12226.

Funding

This study is financially supported by National Key R&D Program of China (grant number 2023YFE0106600) and National Natural Science Foundation of China (grant numbers 22421003 and 22508389).

Declaration of interests

The authors report no conflict of interest.

Appendix A. Details of the controlled simulations

The controlled simulations were performed with mu 0 equals 85.32 $\mu _0=85.32$ mPa s, mu Subscript normal infinity Baseline equals 0.935 $\mu _\infty =0.935$ mPa s and lamda Subscript c Baseline equals 1.29 $\lambda _c=1.29$ s fixed at the values of the original n equals 0.52 $n=0.52$ fluid, while n was varied to 0.62 and 0.68. The resulting critical Reynolds numbers are summarised in table 2. These results confirm that the opposing trends are robust regardless of variations in the other Carreau parameters.

Table 2.

Critical Reynolds numbers obtained from controlled simulations, where only the shear-thinning index n $n$ is varied, while mu 0 comma mu Subscript normal infinity Baseline $\mu _0, \mu _\infty$ and lamda Subscript c $\lambda _c$ are fixed at the values of the original n equals 0.52 $n=0.52$ fluid.

A table summarizing critical Reynolds numbers from controlled simulations with varied shear-thinning index.
Table 2. Long description

A table with three rows and three columns. The columns are labeled n, Re_c,1, and Re_c,2. The row labels and values are as follows: Row 1: n, 0.52; Re_c,1, 94; Re_c,2, 142. Row 2: n, 0.62; Re_c,1, 96; Re_c,2, 132. Row 3: n, 0.68; Re_c,1, 99; Re_c,2, 127.

Appendix B. Validation against linear stability benchmarks

To demonstrate the numerical accuracy of our DNS methodology, we compare our results with the linear stability analysis of Alibenyahia et al. (Reference Alibenyahia, Lemaitre, Nouar and Ait-Messaoudene2012) for power-law fluids, as presented in table 3. The results of Alibenyahia et al. (Reference Alibenyahia, Lemaitre, Nouar and Ait-Messaoudene2012) (shown in their figure 8a) were obtained for a power-law fluid, whereas the present study employs the Carreau model. Despite these differences, both sets exhibit the same physical trend: italic Re Subscript c comma 1 $\textit{Re}_{c,1}$ decreases with decreasing n when defined using the nominal shear-rate viscosity, confirming the reliability of our DNS methodology.

Table 3.

Comparison of the primary instability threshold italic Re Subscript c comma 1 $\textit{Re}_{c,1}$ between the linear stability results of Alibenyahia et al. (Reference Alibenyahia, Lemaitre, Nouar and Ait-Messaoudene2012) using the power-law model and the present DNS using the Carreau model. Both sets exhibit the same physical trend, where decreasing n $n$ leads to a lower italic Re Subscript c comma 1 $\textit{Re}_{c,1}$ , confirming the reliability of our numerical framework.

A table comparing primary instability thresholds between different studies and models.
Table 3. Long description

A table with six rows and four columns comparing the primary instability threshold between the linear stability results of Alibenyahia et al. (2012) using the power-law model and the present DNS using the Carreau model. The columns are labeled Reference, η, n, and Re_c,1. Row 1: Reference Alibenyahia et al. (2012), η 0.8, n 0.7, Re_c,1 ~86. Row 2: Reference Alibenyahia et al. (2012), η 0.8, n 0.6, Re_c,1 ~83. Row 3: Reference Alibenyahia et al. (2012), η 0.8, n 0.5, Re_c,1 ~80. Row 4: Reference Present DNS, η 0.83, n 0.68, Re_c,1 101. Row 5: Reference Present DNS, η 0.83, n 0.62, Re_c,1 97. Row 6: Reference Present DNS, η 0.83, n 0.52, Re_c,1 94.

Figure 5.

Circumferential profiles of u Subscript z Baseline divided by u Subscript i $u_z/u_i$ at the mid-gap ( r Superscript asterisk Baseline equals 0.5 $r^* = 0.5$ ) for n equals 0.52 $n = 0.52$ , 0.62 $0.62$ and 0.68 $0.68$ , evaluated at their respective secondary instability thresholds ( italic Re equals italic Re Subscript c comma 2 Baseline $\textit{Re} = \textit{Re}_{c,2}$ ). All three profiles exhibit three complete wave periods, indicating a dominant azimuthal mode m equals 3 $m = 3$ that is independent of n $n$ .

A line graph showing circumferential profiles for different values of n.

Appendix C. Circumferential profiles and azimuthal mode selection

As shown in figure 5, the circumferential profiles of the axial velocity u Subscript z Baseline divided by u Subscript i $u_z/u_i$ at the mid-gap ( r Superscript asterisk Baseline equals 0.5 $r^* = 0.5$ ) for all three shear-thinning fluids exhibit three complete wave periods, corresponding to a dominant azimuthal mode m equals 3 $m = 3$ . The independence of m $m$ from n $n$ confirms that the stabilising effect originates from enhanced viscous damping in the bulk gap rather than from a change in the preferred azimuthal mode structure.

References

Alibenyahia, B., Lemaitre, C., Nouar, C. & Ait-Messaoudene, N. 2012 Revisiting the stability of circular Couette flow of shear-thinning fluids. J. Non-Newtonian Fluid Mech. 183-184, 37–51.10.1016/j.jnnfm.2012.06.002CrossRefGoogle Scholar
Andereck, C.D., Liu, S.S. & Swinney, H.L. 1986 Flow regimes in a circular Couette system with independently rotating cylinders. J. Fluid Mech. 164, 155–183.10.1017/S0022112086002513CrossRefGoogle Scholar
Ashrafi, N. 2011 Stability analysis of shear-thinning flow between rotating cylinders. Appl. Math. Model. 35, 4407–4423.10.1016/j.apm.2011.03.010CrossRefGoogle Scholar
Brauckmann, H.J. & Eckhardt, B. 2013 Direct numerical simulations of local and global torque in Taylor–Couette flow up to Re = 30 000. J. Fluid Mech. 718, 398–427.10.1017/jfm.2012.618CrossRefGoogle Scholar
Cagney, N. & Balabani, S. 2019 Taylor–Couette flow of shear-thinning fluids. Phys. Fluids 31, 053102.10.1063/1.5088143CrossRefGoogle Scholar
Cagney, N., Lacassagne, T. & Balabani, S. 2020 Taylor–Couette flow of polymer solutions with shear-thinning and viscoelastic rheology. J. Fluid Mech. 905, A28.10.1017/jfm.2020.701CrossRefGoogle Scholar
Cole, J.A. 1976 Taylor-vortex instability and annulus-length effects. J. Fluid Mech. 75, 1–15.10.1017/S0022112076000098CrossRefGoogle Scholar
Coronado-Matutti, O., Souza Mendes, P.R. & Carvalho, M.S. 2004 Instability of inelastic shear-thinning liquids in a Couette flow between concentric cylinders. J. Fluids Engng 126, 385–390.10.1115/1.1760537CrossRefGoogle Scholar
Couette, M. 1890 Études sur le frottement des liquides. Ann. Chim. Phys. 21 (5), 433.Google Scholar
Czarny, O., Serre, E., Bontoux, P. & Lueptow, R.M. 2003 Interaction between Ekman pumping and the centrifugal instability in Taylor–Couette flow. Phys. Fluids 15, 467–477.10.1063/1.1534108CrossRefGoogle Scholar
Dessup, T., Tuckerman, L.S., Wesfreid, J.E., Barkley, D. & Willis, A.P. 2018 Self-sustaining process in Taylor–Couette flow. Phys. Rev. Fluids 3, 123902.10.1103/PhysRevFluids.3.123902CrossRefGoogle Scholar
Dong, S. 2007 Direct numerical simulation of turbulent Taylor–Couette flow. J. Fluid Mech. 587, 373–393.10.1017/S0022112007007367CrossRefGoogle Scholar
Elçiçek, H. & Güzel, B. 2020 Effect of shear-thinning behaviour on flow regimes in Taylor–Couette flows. J. Non-Newtonian Fluid Mech. 279, 104277.10.1016/j.jnnfm.2020.104277CrossRefGoogle Scholar
Escudier, M.P., Gouldson, I.W. & Jones, D.M. 1995 Taylor vortices in Newtonian and shear-thinning liquids. Proc. R. Soc. Lond. A 449, 155–176.10.1098/rspa.1995.0037CrossRefGoogle Scholar
Ji, H., Burin, M., Schartman, E. & Goodman, J. 2006 Hydrodynamic turbulence cannot transport angular momentum effectively in astrophysical disks. Nature 444 (7117), 343–346.10.1038/nature05323CrossRefGoogle ScholarPubMed
Jones, C.A. 1985 The transition to wavy Taylor vortices. J. Fluid Mech. 157, 135–162.10.1017/S0022112085002336CrossRefGoogle Scholar
Kawaguchi, T., Tano, Y. & Saito, T. 2022 Experimental studies on the Taylor–Couette flow of shear-thinning fluids. Intl J. GEOMATE 23 (97), 31–37.10.21660/2022.97.1379CrossRefGoogle Scholar
Khali, S., Nebbali, R. & Bouhadef, K. 2013 Numerical investigation of non-Newtonian fluids flows between two rotating cylinders using lattice-Boltzmann method. Intl J. Mech. Mechatron. Engng 7 (10), 2000–2006.Google Scholar
Larson, R.G. & Desai, P.S. 2015 Modeling the rheology of polymer melts and solutions. Annu. Rev. Fluid Mech. 47, 47–65.10.1146/annurev-fluid-010814-014612CrossRefGoogle Scholar
Lockett, T.J., Richardson, S.M. & Worraker, W.J. 1992 The stability of inelastic non-Newtonian fluids in Couette flow between concentric cylinders: a finite-element study. J. Non-Newtonian Fluid Mech. 43, 165–177.10.1016/0377-0257(92)80023-QCrossRefGoogle Scholar
Majji, M.V., Banerjee, S. & Morris, J.F. 2018 Inertial flow transitions of a suspension in Taylor–Couette geometry. J. Fluid Mech. 835, 936–969.10.1017/jfm.2017.754CrossRefGoogle Scholar
Martinand, D., Serre, E. & Lueptow, R.M. 2014 Mechanisms for the transition to waviness for Taylor vortices. Phys. Fluids 26, 094102.10.1063/1.4895400CrossRefGoogle Scholar
Masuda, H., Horie, T., Hubacz, R., Ohta, M. & Ohmura, N. 2017 Prediction of onset of Taylor–Couette instability for shear-thinning fluids. Rheol. Acta 56, 73–84.10.1007/s00397-016-0987-7CrossRefGoogle Scholar
Mueller, S., Llewellin, E.W. & Mader, H.M. 2010 The rheology of suspensions of solid particles. Proc. R. Soc. A 466, 1201–1228.10.1098/rspa.2009.0445CrossRefGoogle Scholar
Taylor, G.I. 1923 Stability of a viscous liquid contained between two rotating cylinders. Phil. Trans. R. Soc. Lond. A 223, 289–343.10.1098/rsta.1923.0008CrossRefGoogle Scholar
Topayev, S., Nouar, C., Bernardin, D., Neveu, A. & Bahrani, S.A. 2019 Taylor-vortex flow in shear-thinning fluids. Phys. Rev. E 100, 023117.10.1103/PhysRevE.100.023117CrossRefGoogle ScholarPubMed
Topayev, S., Nouar, C. & Dusek, J. 2022 Secondary instabilities in Taylor–Couette flow of shear-thinning fluids. J. Fluid Mech. 933, A4.10.1017/jfm.2021.1036CrossRefGoogle Scholar
Wereley, S.T. & Lueptow, R.M. 1998 Spatio-temporal character of non-wavy and wavy Taylor–Couette flow. J. Fluid Mech. 364, 59–80.10.1017/S0022112098008969CrossRefGoogle Scholar
Xu, F., Su, J., Lan, B., Zhao, P., He, Y., Sun, C. & Wang, J. 2023 Direct numerical simulation of Taylor–Couette flow with vertical asymmetric rough walls. J. Fluid Mech. 975, A30.10.1017/jfm.2023.884CrossRefGoogle Scholar
Xu, F., Zhao, P., Sun, C., He, Y. & Wang, J. 2022 Direct numerical simulation of Taylor–Couette flow: regime-dependent role of axial walls. Chem. Engng Sci. 263, 118075.10.1016/j.ces.2022.118075CrossRefGoogle Scholar
Figure 0

Figure 1. Figure 1 long description.(a) Sketch of the TC system. (b–d) Verification of the Carreau model implementation at Re = 40, corresponding to the CCF regime. Symbols denote the viscosity computed from the DNS at various shear rates for n=0.52$n=0.52$, 0.62$0.62$, and 0.68$0.68$; solid lines represent the theoretical Carreau model (2.5) with parameters from table 1.

Figure 1

Table 1. Carreau model parameters for the three shear-thinning fluids investigated (reproduced from Cagney & Balabani (2019)) together with the critical Reynolds numbers Rec,1$\textit{Re}_{c,1}$ and Rec,2$\textit{Re}_{c,2}$ obtained from the present DNS. Here μ0$\mu _0$ and μ∞$\mu _{\infty }$ are the zero- and infinite-shear-rate viscosities, λc$\lambda _c$ is the characteristic relaxation time and n$n$ is the shear-thinning index.Table 1 long description.

Figure 2

Figure 2. (a) Critical Reynolds numbers for the onset of TVF (Rec,1$\textit{Re}_{c,1}$) and WVF (Rec,2$\textit{Re}_{c,2}$) as functions of the shear-thinning index n$n$. The DNS results are shown as open symbols connected by lines, while the open symbols denote the experimental data of Cagney & Balabani (2019). The regimes CCF, TVF and WVF are denoted by text labels in the diagram. (b) Instantaneous contours of the azimuthal vorticity (normalised by the rotating angular velocity of inner cylinder) in the meridional (r−z$r{-}z$) plane at two azimuthal positions (θ=0$\theta =0$ and θ=π$\theta =\pi$) for TVF (Rec,1=94$\textit{Re}_{c,1}=94$) and WVF (Rec,2=142$\textit{Re}_{c,2}=142$) at n=0.52$n=0.52$. The steady axisymmetric TVF structure and the symmetry breaking in WVF are clearly visible. Three triangular markers at z=L/4$z = L/4$, L/2$L/2$ and 3L/4$3L/4$ are superimposed on the WVF contour at θ=0$\theta =0$, indicating the locations of the three monitor points used in (c); all three probes are at the mid-gap r∗=0.5$r^*=0.5$. (c) Time evolution of the axial velocity uz/ui$u_z/u_i$ at the mid-gap (r∗=0.5$r^*=0.5$) for n=0.52$n=0.52$, recorded at three axial monitor points (z=L/4,L/2,3L/4$z=L/4, L/2, 3L/4$). Dashed lines denote the steady state at the primary instability threshold (Rec,1=94$\textit{Re}_{c,1}=94$), while solid lines illustrate the quasi-periodic oscillations at the secondary instability threshold (Rec,2=142$\textit{Re}_{c,2}=142$). The overlapping traces confirm the existence of an axially coherent, global travelling wave in the WVF regime. (d) Circumferential profile of uz/ui$u_z/u_i$ at the mid-gap (r∗=0.5$r^*=0.5$) for n=0.52$n=0.52$ (the corresponding profiles for other n$n$ values are provided in the Appendix C). The profile exhibits three complete wave periods, corresponding to a dominant azimuthal mode m = 3.

Figure 3

Figure 3. Figure 3 long description.Radial profiles of the time-, azimuthally- and axially-averaged dimensionless viscosity μ(γ˙)/μ0$\mu (\dot \gamma )/\mu _0$ (a,c) and its radial gradient dμ(γ˙)/dr$\text{d}\mu (\dot \gamma )/\text{d}r$ (b,d) for three shear-thinning fluids with n=0.52$n = 0.52$, 0.62$0.62$ and 0.68$0.68$. Panels (a,b) correspond to the primary instability threshold (Re=Rec,1$\textit{Re} = \textit{Re}_{c,1}$); panels (c,d) correspond to the secondary instability threshold (Re=Rec,2$\textit{Re} = \textit{Re}_{c,2}$). The dashed lines in (a,c) represent the average values of dimensionless viscosity at different n, and the dashed lines in (b,d) represent the radial gradient where dμ(γ˙)/dr=0$\text{d}\mu (\dot \gamma )/\text{d}r = 0$.

Figure 4

Figure 4. Contours of the instantaneous (left) shear rate λcγ˙$\lambda _c\dot \gamma$, (middle) dimensionless viscosity μ/μ0$\mu /\mu _0$ and (right) azimuthal velocity uθ/ui$u_\theta /u_i$ superposed velocity vectors (the radial and axial components) in the meridional (r$r$–z$z$) plane for the most shear-thinning fluid (n=0.52$n = 0.52$): (a) at the primary instability threshold (Rec,1=94$\textit{Re}_{c,1} = 94$) and (b) at the secondary instability threshold (Rec,2=142$\textit{Re}_{c,2} = 142$).

Figure 5

Table 2. Critical Reynolds numbers obtained from controlled simulations, where only the shear-thinning index n$n$ is varied, while μ0,μ∞$\mu _0, \mu _\infty$ and λc$\lambda _c$ are fixed at the values of the original n=0.52$n=0.52$ fluid.Table 2 long description.

Figure 6

Table 3. Comparison of the primary instability threshold Rec,1$\textit{Re}_{c,1}$ between the linear stability results of Alibenyahia et al. (2012) using the power-law model and the present DNS using the Carreau model. Both sets exhibit the same physical trend, where decreasing n$n$ leads to a lower Rec,1$\textit{Re}_{c,1}$, confirming the reliability of our numerical framework.Table 3 long description.

Figure 7

Figure 5. Circumferential profiles of uz/ui$u_z/u_i$ at the mid-gap (r∗=0.5$r^* = 0.5$) for n=0.52$n = 0.52$, 0.62$0.62$ and 0.68$0.68$, evaluated at their respective secondary instability thresholds (Re=Rec,2$\textit{Re} = \textit{Re}_{c,2}$). All three profiles exhibit three complete wave periods, indicating a dominant azimuthal mode m=3$m = 3$ that is independent of n$n$.

Supplementary material: File

Xu et al. supplementary material

Xu et al. supplementary material
Download Xu et al. supplementary material(File)
File 730.6 KB