# Taylor-Couette instability in thixotropic yield stress fluids

Mathieu Jenny, Sébastien Kiesgen de Richter, Nicolas Louvet, Salahedine

Skali-Lami, Yvan Dossmann

► To cite this version:

Mathieu Jenny, Sébastien Kiesgen de Richter, Nicolas Louvet, Salahedine Skali-Lami, Yvan Dossmann. Taylor-Couette instability in thixotropic yield stress fluids. Physical Review Fluids, 2017, 2, pp.023302. 10.1103/PhysRevFluids.2.023302. [hal-01491186]

HAL Id: hal-01491186

https://hal.science/hal-01491186

Submitted on 16 Mar 2017

**HAL** is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire **HAL**, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Taylor-Couette instability in thixotropic yield stress fluids

Mathieu Jenny,\* Sébastien Kiesgen de Richter, Nicolas Louvet, Salahedine Skali-Lami, and Yvan Dossmann

LEMTA, UMR No. 7563, Université de Lorraine, CNRS, F-54500 Vandoeuvre-lès-Nancy, France
(Received 30 August 2016; published xxxxxx)

We consider the flow of thixotropic yield stress fluids between two concentric cylinders. To account for the fluid thixotropy, we use Houška's model [Houška, Ph.D. thesis, Czech Technical University, Prague, 1981] with a single structural parameter driven by a kinetic equation. Because of the yield stress and the geometric inhomogeneity of the stress, only a part of the material in the gap may flow. Depending on the breakdown rate of the structural parameter, the constitutive relation can lead to a nonmonotonic flow curve. This nonmonotonic behavior is known to induce a discontinuity in the slope of the velocity profile within the flowing material, called shear banding. Thus, for fragile structures, a shear-banded flow characterized by a very sharp transition between the flowing and the static regions may be observed. For stronger structures, the discontinuity disappears and a smooth transition between the flowing and the static regions is observed. The consequences of the thixotropy on the linear stability of the azimuthal flow are studied in a large range of parameters. Although the thixotropy allows shear banding in the base flow, it does not modify fundamentally the linear stability of the Couette flow compared to a simple yield stress fluid. The apparent shear-thinning behavior depends on the thixotropic parameters of the fluid and the results about the onset of the Taylor vortices in shear-thinning fluids are retrieved. Nevertheless, the shear banding modifies the stratification of the viscosity in the flowing zone such that the critical conditions are mainly driven by the width of the flowing region.

DOI: 10.1103/PhysRevFluids.00.003300

## I. INTRODUCTION

Yield stress fluids, such as emulsions, foams, mud, and gels, are of industrial interest. Because of their high number of applications, they have been intensively studied over the few past decades. Many of them have an inner microstructure, responsible for the yield stress when an external load is applied, that resists large-scale rearrangement. The destruction of this microstructure by the flow is responsible for a complex phenomenon, named thixotropy. The competition between the internal reorganization and the macroscopic flow induces a complex time dependence of the rheological parameters, such as the apparent viscosity. The complex behavior of such fluids raises the question of their flow stability in industrial conditions (melting, mass transfers, etc.).

The Taylor-Couette flow is often considered as a paradigm to study the stability and the transition to turbulence of complex fluids [1–4]. While the Taylor-Couette flow of Newtonian fluids has been extensively studied since the historical work of Taylor [5], much attention has been paid to complex fluids during the past decade. According to the studies of the hydrodynamic stability of shear thinning [1,3,6] and Bingham fluids [2–4], it is observed that when the viscosity is scaled with the inner-wall shear viscosity, shear thinning has a stabilizing effect, i.e., the appearance of the Taylor vortices is delayed [3]. For simple yield stress fluids, two regions of the flow coexist, a yielded zone close to the rotating cylinder and a static region close to the fixed cylinder. Landry et al. [2] have shown that the vortices are localized in the yield zone. Naimi et al. [7] have reported that the yield stress appears to stabilize the flow. Few studies have focused on the influence of the thixotropy on the stability of Taylor-Couette flow in thixotropic shear-thinning fluids [8]. Questions remain about the consequences of a microstructural-dependent yield stress on the stability of the flow.

<footnote>\*mathieu.jenny@univ-lorraine.fr</footnote>
<figure><img src="image_2.png" /><figcaption>FIG. 1. Taylor-Couette geometry. Here $e_r, e_\theta$, and $e_\zeta$ are the unit vectors of the cylindrical coordinate system $(r, \theta, z)$.</figcaption></figure>

In that context, we study of the stability of a thixotropic yield stress fluid in a Taylor-Couette configuration. To model the base flow of such fluids, structural parameter models allow us to take into account the inner dynamic, coupled with the surrounding flow. The inner structure is then entirely described by the structural parameter. In this article we use Houška’s model [9,10]. It is built from the Hershel-Bulkley model, commonly used for nonelastic yield stress fluids, considering that the consistency $K$ and the yield stress $\tau_0$ depend linearly on the structural parameter $\lambda$. This robust thixotropic fluid model was originally developed to characterize liquid foods such as ketchup or yogurt [9–11]. The existence and unicity of a steady solution of this model was recently established in pipe flows [12]. Houška’s model has been successfully used by Wachs et al. [13] to model start-up pipe flows of waxy crude oils despite the limitations of the model. In particular, irreversible effects are not accounted for, as stressed by Mendes et al. [14]. Billingham and Ferguson [15] also investigated also steady pipe flow of bentonite mud using Houška’s model. The simplicity of the model allows us to conclude about the consequences of the thixotropy upon the linear stability of the flow. The linear stability analysis of the flow shows that the nature of the linear unstable mode is steady and axisymmetric in the large range of the explored parameters and does not depend on the thixotropic character of the flow.

## II. CONSTITUTIVE EQUATIONS AND NUMERICAL APPROACH

### A. Flow geometry

In the present work we consider the case of an inner cylinder rotating at a given angular velocity $\hat{\omega}_i$ so that the velocity at the inner radius $\hat{r}_i$ is $\hat{v}_i = \hat{\omega}_i \hat{r}_i$ (Fig. 1). Here the caret denotes a dimensional variable. The outer cylinder is static. In this configuration, when the velocity of the inner cylinder is sufficiently low, the purely azimuthal steady flow is stable for viscous fluids [3,5,6,16].

### B. Houška’s model

For a nonzero strain rate, i.e., $\hat{\gamma} \neq 0$, the constitutive law of Houška’s model [9,10] is expressed as

$$ \overline{\overline{\hat{\tau}}} = [(K + \Delta K \lambda)\hat{\gamma}^{n_c} + \tau_0 + \tau_1 \lambda] \frac{\overline{\overline{\hat{\gamma}}}}{\hat{\gamma}}, \quad (1) $$
where $\overline{\overline{\hat{\tau}}}$ denotes the stress tensor, the parameters $\Delta K$ and $\tau_1$ respectively determine the sensitivity of the consistency and the yield stress with the structural parameter $\lambda$, and $n_c$ is the shear-thinning index. The strain rate tensor is given by $\overline{\overline{\hat{\gamma}}} = \overline{\overline{\nabla}}\hat{\nu} + \overline{\overline{\nabla}}\hat{\nu}^T$, where $\overline{\overline{\nabla}}\hat{\nu}$ and the index $T$ denote the gradient tensor of the velocity vector $\hat{\nu}$ and the transposition, respectively. The second invariant of the strain rate tensor $\overline{\overline{\hat{\gamma}}}$ is defined by

$$ \hat{\gamma} = (\frac{1}{2}\hat{\gamma}_{ij}\hat{\gamma}_{ij})^{1/2}, \qquad (2) $$

using the Einstein summation convention where the elements of the strain tensor $\hat{\gamma}_{ij}$ are defined with the components of the fluid velocity $\hat{\nu}$. In Eq. (2), the indices $i$ and $j$ stand for the cylindrical coordinates $r, \theta$, and $z$ (Fig. 1). The structural parameter $\lambda$ is determined by the kinetic equation

$$ \frac{\partial \lambda}{\partial \hat{\tau}} + \hat{\nu} \cdot \nabla \lambda = a(1 - \lambda) - b\lambda \hat{\gamma}^m, \qquad (3) $$

where $a$ and $b$ are, respectively, the building and the breakdown parameters. The thixotropic breakdown index $m$ is taken to be equal to 1 in the following. The values of the structural parameter are within the range $0 \le \lambda \le 1$. The value $\lambda = 1$ means that the fluid is fully structured and $\lambda = 0$ means that it is fully unstructured. The kinetic equation (3) governs the evolution of the microstructure, which influences the flow by modifying the stress tensor.

### C. Nondimensional equations

To nondimensionalize the constitutive equations of the flow in a cylindrical Couette geometry, we choose the following references for the density, the velocity, and the length, respectively:

$$ \rho_{\text{ref}} = \rho, \quad v_{\text{ref}} = \hat{v}_i, \quad l_{\text{ref}} = \hat{d}, \qquad (4) $$

where $\rho$ is the density of the fluid and $\hat{v}_i$ the inner cylinder velocity. One can build a characteristic strain rate $\hat{v}_i/\hat{d}$ using the latter reference dimensional parameters. For the non-Newtonian fluids, several choices can be made for the reference viscosity. We choose the plastic viscosity of the fluid at the characteristic strain rate $\hat{v}_i/\hat{d}$ as a reference viscosity

$$ \mu_{\text{ref}} = \mu_0(1 + \Delta K^* \lambda_{\text{ref}}). \qquad (5) $$

The parameter $\mu_0 = K(\hat{v}_i/\hat{d})^{n_c-1}$ is the standard reference viscosity of a power-law fluid. The structural parameter $\lambda_{\text{ref}}$ is given by Eq. (3) at equilibrium

$$ \lambda_{\text{ref}} = \frac{a}{a + b(\hat{v}_i/\hat{d})^m} = \frac{1}{1 + b^*/a^*}, \qquad (6) $$

where the nondimensional building and breakdown parameters are

$$ a^* = \frac{a\hat{d}}{\hat{v}_i}, \quad b^* = b\left(\frac{\hat{v}_i}{\hat{d}}\right)^{m-1}, \qquad (7) $$

respectively. Here $\Delta K^* = \Delta K/K$ is the reduced thixotropic consistency factor. This parameter characterizes the dependence of the plastic viscosity with the inner structure of the fluid in comparison with the intrinsic consistency $K$, which itself depends on the solvent. The reference viscosity $\mu_{\text{ref}}$, depending on the ratio of the breakdown parameter $b^*$ over the building parameter $a^*$, decreases when $b^*/a^*$ increases, i.e., when the inner structure of the fluid becomes more and more fragile.
Using the previous reference dimensions, the Navier-Stokes and mass conservation equations for flows of incompressible fluids are

$$ \frac{\partial \mathbf{v}}{\partial t} + (\mathbf{v} \cdot \nabla)\mathbf{v} = -\nabla p + \frac{1}{\mathrm{Re}} \nabla \cdot \bar{\bar{\tau}}, \qquad (8) $$

$$ \nabla \cdot \mathbf{v} = 0, \qquad (9) $$

where $\mathbf{v} = \hat{\mathbf{v}}/\hat{v}_i$ stands for the reduced velocity and $p = \hat{p}/\rho \hat{v}_i^2$ for the reduced pressure. One can notice that $\hat{p}$ is the modified pressure including the hydrostatic pressure. The Reynolds number is defined using the reference viscosity (5) by

$$ \mathrm{Re} = \frac{\mathrm{Re}_0}{1 + \Delta K^* \lambda_{\mathrm{ref}}}, \qquad (10) $$

where $\mathrm{Re}_0 = \frac{\rho \hat{v}_i \hat{d}}{\mu_0}$. Thus, the reduced stress tensor reads

$$ \bar{\bar{\tau}} = \left[ \left( \frac{1 + \Delta K^* \lambda}{1 + \Delta K^* \lambda_{\mathrm{ref}}} \right) \dot{\gamma}^{n_c} + \mathrm{Bn} \left( \frac{1 + \tau_1^* \lambda}{1 + \tau_1^* \lambda_{\mathrm{ref}}} \right) \right] \frac{\bar{\bar{\dot{\gamma}}}}{\dot{\gamma}}, \qquad (11) $$

where $\dot{\gamma}$ and $\bar{\bar{\dot{\gamma}}}$ are the nondimensional strain rate and strain tensor, respectively, and $\tau_1^* = \tau_1/\tau_0$ is the reduced thixotropic yield stress. Equation (11) involves the Bingham number, which is the ratio of the yield stress with the plastic viscous stress

$$ \mathrm{Bn} = \mathrm{Bn}_0 \frac{1 + \tau_1^* \lambda_{\mathrm{ref}}}{1 + \Delta K^* \lambda_{\mathrm{ref}}}, \qquad (12) $$

where $\mathrm{Bn}_0 = \frac{\tau_0}{K(\hat{v}_i/\hat{d})^{n_c}}$ is the standard Bingham number of a Hershel-Bulkley fluid. The Bingham number increases with the yield stress and localization stops occurring for sufficiently high values of the Bingham number. According to Eq. (11), the reduced yield stress equals

$$ \tau_y = \mathrm{Bn} \left( \frac{1 + \tau_1^* \lambda}{1 + \tau_1^* \lambda_{\mathrm{ref}}} \right). \qquad (13) $$

For the structural parameter, the nondimensional version of Eq. (3) is

$$ \frac{\partial \lambda}{\partial t} + \mathbf{v} \cdot \nabla \lambda = a^*(1 - \lambda) - b^* \lambda \dot{\gamma}^m. \qquad (14) $$

Equations (11) and (14) describe the coupling between the flow properties and the evolution of the microstructure.

## D. Boundary conditions for the flow

The inner and outer reduced radii of the Couette setup are defined by

$$ r_i = \frac{\eta}{1 - \eta}, \qquad (15) $$

$$ r_e = \frac{1}{1 - \eta}, \qquad (16) $$

with $\eta = \hat{r}_i/\hat{r}_e$ the radii ratio. The velocity vector $\mathbf{v}$ is written in the cylindrical basis as $\mathbf{v} = v_r \mathbf{e}_r + v_\theta \mathbf{e}_\theta + v_z \mathbf{e}_z$. The boundary conditions are as follows.

(i) At the inner radius $r = r_i$, the velocity components are $v_\theta = 1$ and $v_r = v_z = 0$.

(ii) At the outer radius of the flowing zone $r = r_o$, the velocity components are $v_r = v_\theta = v_z = 0$.

(iii) In our case, there is a material limit at $r = r_e$. Thus, the outer radius $r_o$ is given by the following criterion: If $\tau(r_e) \ge \tau_y$, $r_o = r_e$; otherwise $\tau(r_o) = \tau_y$. In the following, we assume that
the stress at the interface between the flowing and static regions is the yield stress. Other assumptions

would be beyond the framework of Houška's model and it would demand a model for the solid phase.

At $r = r_o$, the yield stress $\tau_y$ given by Eq. (13) becomes

$$ \tau_{yo} = \text{Bn} \left( \frac{1 + \tau_1^* \lambda_o}{1 + \tau_1^* \lambda_{\text{ref}}} \right), \qquad (17) $$

with $\lambda_o$ the structural parameter at the interface. Thus, the last boundary condition becomes $\tau(r_o) =$

$\tau_{yo}$ and the stress $\tau_{yo}$ at the interface is defined by Eq. (17). Nevertheless, the stress condition at

the interface between the fluid and solidlike zone is well defined only if the structural parameter

$\lambda$ is continuous across the interface. When shear banding occurs, we expect a discontinuity of the

structural parameter $\lambda$ at the interface when no diffusion of the structural parameter $\lambda$ is included

in Eq. (14). According to Olmsted et al. [17], a stress diffusion term must be added in that case to

conserve the unicity of the steady solution by selecting the stress at the interface between the bands

(see also Lu et al. [18]). It was also shown that a spatially local model, i.e., without any diffusive

gradient of the stress (diffusive term for the structural parameter in our case), will not correctly

predict a shear banded state. The steady state depends then on the flow or numerical noise history

by selecting arbitrarily a stress value at the interface. When adding a diffusive term, the continuity

of the yield stress or the structural parameter across the interface is ensured. This kind of diffusive

term was recently interpreted as a nonlocal effect at the molecular scale in the flow of micellar

suspensions [19]. As the value of the stress diffusion coefficient, similar to the structural diffusion

coefficient, is found to be very small [19,20], we will focus in the following on the cases where the

coefficient of diffusion is equal to zero. This assumption implies that the stress interface is fixed. For

Bingham-like fluids, i.e., when $n_c = 1$ and $m = 1$, Eqs. (11) and (14) for one-dimensional steady

flows give the following equation for $\lambda$ over the gap:

$$ (b^* \tilde{\tau}_1 - a^* \Delta \tilde{K}) \lambda^2 + [a^*(\Delta \tilde{K} - \tilde{K}) + b^*(\tilde{\tau}_0 - \tau)] \lambda + a^* \tilde{K} = 0, \qquad (18) $$

where

$$ \tilde{K} = \frac{1}{1 + \Delta K^* \lambda_{\text{ref}}}, \qquad \Delta \tilde{K} = \frac{\Delta K^*}{1 + \Delta K^* \lambda_{\text{ref}}}, \qquad (19) $$

$$ \tilde{\tau}_0 = \frac{\text{Bn}}{1 + \tau_1^* \lambda_{\text{ref}}}, \qquad \tilde{\tau}_1 = \frac{\text{Bn} \tau_1^*}{1 + \tau_1^* \lambda_{\text{ref}}}. \qquad (20) $$

Assuming that the stress at the interface is the yield stress $\tau_{yo}$, the second-order polynomial (18) is

rewritten at $r = r_o$ setting $\tau = \tau_{yo}$ and it becomes

$$ -a^* \Delta \tilde{K} \lambda_o^2 + a^*(\Delta \tilde{K} - \tilde{K}) \lambda_o + a^* \tilde{K} = 0. \qquad (21) $$

The only positive root of (21) is $\lambda_o = 1$. In the framework of the considered model without any

diffusive term, the stress at the fluid-solid interface is always the yield stress of the fully structured

material $\tau_{ys}$ in steady flows. Thus, we can either solve the steady equations with or without shear

banding, setting the stress at the interface

$$ \tau(r_o) = \tau_{ys} = \text{Bn} \left( \frac{1 + \tau_1^*}{1 + \tau_1^* \lambda_{\text{ref}}} \right). \qquad (22) $$

Note that adding a diffusive term in Eq. (14) would result in slightly different stress values. It might

be necessary to compare to experimental results, but it would not modify the conclusions of the

present paper.

## E. Numerical methods for steady flows

Only the flowing region needs to be considered to solve the steady flow. To perform the numerical

resolution, we use a finite-difference method for the spatial discretization. The mesh points are
regularly spaced between the inner radius $r_i$ and the outer radius $r_o$. The stress points are taken between two successive velocity points to ensure the numerical accuracy of the scheme for the velocity. For the derivative operations, the standard second-order centered scheme is used. The numerical method used for the spatial discretization is quite well established and is similar to the ones used, for instance, in Refs. [8,21,22]. A validation and a convergence test are performed in Sec. IV B.

To calculate the base flow, we consider the steady axisymmetric solution of Eqs. (8), (9) and (14), i.e., $v_b = V_b(r)e_\theta$ and $\lambda = \lambda_b(r)$. In the flowing region, i.e., $r_i \le r \le r_o$, the only nonzero element of the strain rate tensor is $\dot{\gamma}_{r\theta}$. The strain rate is always nonzero, negative in the flowing region. The only nonzero element of the stress tensor $\bar{\tau}$ is then $\tau_{r\theta}$. Considering the previous assumptions for the flow, the well known result for steady Couette flow applies:

$$ \tau_{r\theta,b} = -\frac{C}{r^2}, \qquad (23) $$

where the positive constant $C$ is related to the torque imposed by the inner rotating cylinder. The radius $r_o$ can be obtained from the stress condition on the interface between the yielded and unyielded regions:

$$ r_o = \sqrt{\frac{C}{\tau_{yo}}}. \qquad (24) $$

If $r_o \ge r_e$ according to Eq. (24), all the material in the gap flows and $r_o = r_e$. In the next section, the flow curves show that the minimal value $\tau_{min}$ of the stress may be below $\tau_{yo}$. Thus, for $\tau_{min}r_e^2 \le C < \tau_{yo}r_e^2$, an alternative to Eq. (24) is to set $r_o = r_e$. In practice, this means that if there is no interface at the initial state, the fluid region fits the whole gap for $\tau_{min}r_e^2 \le C < \tau_{yo}r_e^2$ and if there is a solidlike region in the initial state, the flowing region is confined between $r_i$ and $\dot{r}_o < r_e$ according Eq. (24).

To compute the flow velocity in the yielded region, we calculate the strain rate by solving the regular setup of equations at each point of the mesh:

$$ \lambda_b = \frac{1}{1 + (b^*/a^*)\dot{\gamma}_b^m}, \qquad (25) $$

$$ \dot{\gamma}_b = \left( \frac{(C/r^2 - \tau_{yb})(1 + \Delta K^* \lambda_{ref})}{1 + \Delta K^* \lambda_b} \right)^{1/n_c}, \qquad (26) $$

with $\tau_{yb}$ the yield stress given by Eq. (13) replacing $\lambda$ by $\lambda_b$. Once we obtain the strain rate $\dot{\gamma}_b$ and the structural parameter $\lambda_b$ for a given constant $C$ by solving the setup of Eqs. (25) and (26), the fluid velocity is calculated by the integration of the strain rate with $V_b(r_o) = 0$ as the boundary condition. Finally, one has to find the value of $C$ such that $V_b(r_i) = 1$ using the algorithm available in MATLAB to calculate the zero of a real nonlinear function. If needed, the pressure $P_b$ of the base flow can be obtained by integrating the equation

$$ \frac{\partial P_b}{\partial r} = \frac{V_b^2}{r} \qquad (27) $$

and setting the inner pressure $P_b(r_i) = 0$ for instance.

### III. EFFECT OF THE THIXOTROPY ON THE BASE FLOW

In the following sections, we set the thixotropic index breakdown to $m = 1$. It seems reasonable to argue that the structural parameter $\lambda$ modifies the viscous term and the yield stress with the same order of magnitude. Thus, we set $\Delta K^* = \tau_1^*$ and so $Bn = Bn_0$. As we focus our study on the cases where both flowing and solidlike regions exist, i.e., $r_o < r_e$, the Bingham number is set equal to
<figure><img src="image_3.png" /><figcaption>FIG. 2. Composite flow curves, stress $\tau$ vs the strain rate $\dot{\gamma}$, become nonmonotonic depending on the breakdown parameter $b^*$ (Bn = 2, $n_c$ = 1, $\Delta K^* = 1$, $\tau_1^* = 1$, and $a^* = 1$).</figcaption></figure>

Bn = 2. Finally, for steady-state flows only the ratio $b^*/a^*$ appears and we choose to set $a^* = 1$ without loss of generality.

A. Steady-state flow curves

The base flow is computed using 50 grid points in the flowing region of the gap to ensure good accuracy for the linear stability analysis as shown in Sec. IV B. We report in Fig. 2 the evolution of the composite flow curves under controlled shear rate for different values of the breakdown parameter $b^*$. The composite curves are obtained straightforwardly by replacing the structural parameter $\lambda$ by its relation to $\dot{\gamma}$ (25) in the constitutive law (11). As $b^*$ increases, the composite curve drops from a monotonic to a nonmonotonic behavior, which presents an unstable branch leading to shear banding [23–25]. This result indicates that, in the range of parameters studied here, shear banding may occur in the base flow. Shear banding is characterized by a discontinuous strain rate $\dot{\gamma}$ and structural parameter $\lambda$ across the fluid-solid interface [Figs. 3(b) and 3(d)], although in the simple localized flow, the transition between the flowing and static regions is smooth [Figs. 3(a) and 3(c)]. The mechanism underlying shear banding is different from those considering viscoelastic fluids as in the major classical works about the stability of shear-banded flows [26]. Structural discontinuity is here the underlying cause of shear banding. The localized flow, due to the yield stress, is also observed in simple yield stress fluids such as Bingham fluids.

Thus, the shear banding may appear if the sign of the derivative of the constitutive relation $\tau(\dot{\gamma})$ changes at a critical strain rate $\dot{\gamma}_0 > 0$. In other words, the necessary condition to allow the shear-banded flow is

$$ \exists \dot{\gamma}_0 \ge 0: \frac{\partial \tau}{\partial \dot{\gamma}} = 0. \qquad (28) $$

In the steady state, the derivative of the stress given by Houška's model is

$$ \frac{\partial \tau}{\partial \dot{\gamma}} = \dot{\gamma}^{n_c-1} \left( \frac{n_c(b^*/a^*)^2 \dot{\gamma}^{2m} + (b^*/a^*)[2n_c + (n_c - m)\Delta K^*]\dot{\gamma}^m + n_c(1 + \Delta K^*)}{(1 + \Delta K^*\lambda_{\text{ref}})[1 + (b^*/a^*)\dot{\gamma}^m]^2} - \frac{m(b^*/a^*)\text{Bn}\tau_1^*\dot{\gamma}^{m-n_c}}{(1 + \tau_1^*\lambda_{\text{ref}})[1 + (b^*/a^*)\dot{\gamma}^m]^2} \right). \qquad (29) $$
<figure><img src="image_4.png" /><figcaption>FIG. 3. (a) and (b) Strain rate $\dot{\gamma}$ in the gap and (c) and (d) stress and yield stress in the gap in (a) and (c) a localized flow at $b^* = 0.5$ and (b) and (d) a shear-banded flow at $b^* = 2$. Other parameters are the Bingham number $Bn = 2$, the shear-thinning index $n_c = 1$, the structural dependence of the consistency $\Delta K^* = 1$, the structural dependence of the yield stress $\tau_1^* = 1$, and the building parameter $a^* = 1$.</figcaption></figure>

213 For Bingham-like fluids where $m = n_c = 1$, Eq. (28), using the derivative of $\tau$ given by Eq. (29), 214 admits only one positive root $\dot{\gamma}_0$:

$$ \dot{\gamma}_0 = \frac{a^*}{b^*} \left[ \sqrt{Bn\tau_1^* \frac{b^*}{a^*} \left( \frac{1 + \Delta K^* \lambda_{\text{ref}}}{1 + \tau_1^* \lambda_{\text{ref}}} \right) - \Delta K^* - 1} \right] \quad (30) $$

215 if

$$ 1 + \Delta K^* - Bn\tau_1^* \frac{b^*}{a^*} \left( \frac{1 + \Delta K^* \lambda_{\text{ref}}}{1 + \tau_1^* \lambda_{\text{ref}}} \right) < 0. \quad (31) $$

216 The derivative of $\tau$ [Eq. (29)] is negative for $\dot{\gamma} \in [0, \dot{\gamma}_0[$. As $\dot{\gamma} = 0$ lies in the forbidden range of 217 strain rates, the simple localization is no longer stable. The strain rate at the structural discontinuity 218 is the nonzero strain rate

$$ \dot{\gamma}_c = Bn\tau_1^* \left( \frac{1 + \Delta K^* \lambda_{\text{ref}}}{1 + \tau_1^* \lambda_{\text{ref}}} \right) - \frac{1 + \Delta K^*}{b^*/a^*}, \quad (32) $$

219 which ensures the stress condition (22). Figure 2 and Eqs. (30) and (32) show that $\dot{\gamma}_c \ge \dot{\gamma}_0$.
<figure><img src="image_5.png" /><figcaption>FIG. 4. (a) Strain rate $\dot{\gamma}_b$, (b) viscosity $\mu_b$, (c) azimuthal velocity $V_b$, and (d) structural parameter $\lambda_b$ of the base flow vs the reduced gap position $y = (r - r_i)/(r_e - r_i)$ with $Bn = 2$, $n_c = 1$, $\Delta K^* = 1$, $\tau_1^* = 1$, and $a^* = 1$ for a large gap $\eta = 0.5$.</figcaption></figure>

Note that if $n_c < 1$ and $m \ge 1$,

$$ \frac{\partial \tau}{\partial \dot{\gamma}_{\dot{\gamma} \to 0^+}} \to +\infty \qquad (33) $$

and thus there is a range of strain rate values close to zero where $\tau(\dot{\gamma})$ is a growing function of $\dot{\gamma}$.

The range of positive strain rates where the stress $\tau$ decreases cannot start at a zero value. This case would be similar to the flow curve of a semidilute wormlike micelle solution with a yield stress like in Fig. 1(b) of [27]. As there is no diffusion term in our set of equations, sharp discontinuities of the strain rate and the structural parameter can appear within the fluid region at a radius $r_i < r < r_o$ when $n_c < 1$. Our numerical method does not allow such discontinuous fields in the fluid domain except at the interface between the fluid and solidlike region, i.e., at $r = r_o$. Thus, in the following, we limit our parametric study to shear-banded flows with an interface between the flowing and static regions only [flow curve corresponding to Fig. 1(c) of [27]], i.e., with $n_c = 1$ and shear localization with $n_c \le 1$.

## B. Velocity profiles and structural parameter

The strain rate, the viscosity, the velocity, and the structural parameter profiles are shown in Figs. 4(a)–4(d). We see that for all values of $b^*$, a flowing and a static region coexist with the considered value of Bingham number $Bn = 2$. However, the discontinuity of the strain rate profile
depends on the latter parameter and is related to a discontinuity of the structural parameter. For lower values of $b^*$, the velocity profiles as the structure parameter is continuous as observed, for example, in carbopol gels, emulsions, and foams [28,29]. In that case, the shear localization is inherent to the existence of the yield stress and no shear banding is observed. By contrast, for $b^*$ greater than the critical value $b_c^* = 1$, the shear rate becomes nonzero at the outer boundary of the flowing region [Fig. 4(a)]. The underlying cause of this discontinuity of the strain rate is the discontinuity of the structural parameter $\lambda$ across the boundary of the flowing region [Fig. 4(d)]. Such a discontinuous strain rate profile between a static and a flowing region has been observed using magnetic resonance imaging measurements in cement pastes [30] and bentonite suspensions [29]. It corresponds to a steady-state shear-banded velocity profile where the shear rate is equal to a critical shear rate in the liquid region and is equal to zero in the solid region.

According to Eq. (31), increasing the parameter $\tau_1^*$ may produce the same effect as increasing $b^*/a^*$. The steady-state flow is controlled by the competition between the restructuring and the breakdown effects. The more the structure close to the interface is broken efficiently by the strain rate [Fig. 4(d)], the more the viscosity drops significantly and rapidly.

Finally, we explore the effect of the shear-thinning index $n_c$. When $n_c < 1$, the shear-banded flows are not observed because the constitutive relation of the material is always a growing function for $\dot{\gamma}$ sufficiently close to zero [Eq. (33)], allowing small values for $\dot{\gamma}$ in the flow. In that case the flow is always simply shear localized and smooth. This contrasts with the previous cases discussed above ($n_c = 1$), where small values of $\dot{\gamma}$ fall in the unstable branch of the flow curve and then lead to shear-banded flows. As would be expected from the velocity profiles obtained with Carreau fluids by Alibenyahia et al. [3], the flow is confined close to the inner cylinder when the shear-thinning index $n_c$ decreases. This confirms that the shear-thinning behavior confines the flow in the inner region of the gap where the viscosity is lower.

### C. Interface between the static and the flowing regions

Now we focus on the evolution of the width of the flowing region $y_o = r_o - r_i$ depending on the thixotropic parameters and the shear-thinning index $n_c$; $y_o$ can be obtained in Fig. 4(d) by reading the abscissa where $\lambda$ reaches 1. In Figs. 5(a)–5(c) the curve of $y_o$ separates the inner flowing region from the outer static region. The base flow evolves smoothly from shear localization to shear-banding regimes by increasing the value of $b^*$ or $\Delta K^*$ and $\tau_1^*$ [Figs. 5(a) and 5(b)]. Nevertheless, it can be observed in Fig. 5(a) that $y_o$ decreases faster in shear-localized flows than in shear-banded flows where $b^*$ or $\Delta K^*$ and $\tau_1^*$ increase. When the breakdown parameter $b^*$ tends to infinity, the size of the flowing region, characterized by the reduced position $y_o = r_o - r_i$ of the interface between the flowing and solidlike regions, decreases to a minimum size corresponding to the one of the fully unstructured equivalent Bingham fluid [Fig. 5(a)]. A quite similar remark can be made about $\Delta K^*$ and $\tau_1^*$. Indeed, increasing the parameters $b^*$, $\Delta K^*$, or $\tau_1^*$ makes the shear-thinning behavior stronger in the steady flowing region. Thus, it is not surprising that when the shear-thinning index $n_c$ decreases, the width of the flowing zone also decreases [Fig. 5(c)]. Nevertheless, the model with one structural parameter predicts shear banding when the thixotropic parameters grow above some critical values given by Eq. (31). The regular Bingham or Hershel-Bulkley laws only describe the shear localization.

Now we will study the linear stability of the base flow, whether shear banding is present or not.

## IV. LINEAR STABILITY ANALYSIS

### A. Equation setup

To perform a linear analysis of stability, the fluid velocity, the structural parameter, and the pressure are decomposed as

$$v = v_b + \tilde{v}(r)\exp(\sigma t + in\theta + ikz), \quad (34)$$

$$\lambda = \lambda_b + \tilde{\lambda}(r)\exp(\sigma t + in\theta + ikz), \quad (35)$$

$$p = P_b + \tilde{p}(r)\exp(\sigma t + in\theta + ikz), \quad (36)$$
<figure><img src="image_6.png" /><figcaption>FIG. 5. Reduced position of the interface between flowing and static regions $y_o = r_o - r_i$ for a large gap $\eta = 0.5$ with $B_n = 2$, $n_c = 1$, and $a^* = 1$. (a) $y_o$ vs $b^*$ with $\Delta K^* = 1$ and $\tau_1^* = 1$. The vertical dashed line stands for the critical value of $b^* = 1$ where the strain rate at the interface $\dot{\gamma}_o$ becomes nonzero. (b) $y_o$ vs $\Delta K^*$ or $\tau_1^*$ with $b^* = 1$. The vertical dashed line stands for the critical value of $\Delta K^* = 1$ (or $\tau_1^*$) where the strain rate at the interface $\dot{\gamma}_o$ becomes nonzero. (c) $y_o$ vs $n_c$ with $\Delta K^* = 0.5$, $\tau_1^* = 0.5$, and $b^* = 1$. Only shear-localization cases are considered when $n_c \neq 1$.</figcaption></figure>

$$ \sigma \tilde{\mathbf{v}} = -\overline{\tilde{\nabla}} \mathbf{v}_b \cdot \tilde{\mathbf{v}} - \overline{\tilde{\nabla}} \tilde{\mathbf{v}} \cdot \mathbf{v}_b + \frac{1}{\mathrm{Re}} \mathrm{div} \left( \left. \frac{\partial \overline{\tilde{\mathbf{\tau}}}}{\partial \dot{\gamma}_{ij}} \right|_b \dot{\gamma}_{ij}(\tilde{\mathbf{v}}) + \left. \frac{\partial \overline{\tilde{\mathbf{\tau}}}}{\partial \lambda} \right|_b \tilde{\lambda} \right) - \nabla \tilde{p}, \quad (37) $$

$$ \sigma \tilde{\lambda} = -\mathbf{v}_b \cdot \nabla \tilde{\lambda} - \tilde{\mathbf{v}} \cdot \nabla \lambda_b - (a + b \dot{\gamma}_b^m) \tilde{\lambda} - mb \lambda_b \dot{\gamma}_b^{m-1} \left. \frac{\partial \dot{\gamma}}{\partial \dot{\gamma}_{ij}} \right|_b \dot{\gamma}_{ij}(\tilde{\mathbf{v}}), \quad (38) $$

$$ 0 = \mathrm{div}(\tilde{\mathbf{v}}). \quad (39) $$

<footnote n="279">279 where $\tilde{\mathbf{v}}$, $\tilde{\lambda}$, and $\tilde{p}$ are the perturbation of the base flow considering the azimuthal mode $n$ and the</footnote>

<footnote n="280">280 axial wave number $k$. Injecting Eqs. (34)–(36) in the general setup of Eqs. (8), (14), and (9) and after</footnote>

<footnote n="281">281 withdrawing the nonlinear terms, the linear setup of equations for the perturbation of the base flow is</footnote>

<footnote n="282">282 The indices $ij$ stand for $r, \theta$, or $z$ and Einstein's convention for summation is used. As the stress is</footnote>

<footnote n="283">283 always continuous across $r_o$, we can use a method similar to those of Frigaard et al. [31] and Landry</footnote>

<footnote n="284">284 et al. [2]. It consists in writing the linear setup of Eqs. (37)–(39) in the flowing region of the base</footnote>

<footnote n="285">285 flow. The displacement of the yield stress boundary is fully driven by the perturbation of the flow</footnote>
<table><caption>TABLE I. Critical Reynolds number $Re_c$ and critical axial wave number $k_c$ for a Newtonian fluid vs the number of nodes $M$ in the gap with a radii ratio $\eta = 0.9$.</caption><thead><tr><th><i>M</i></th><th>Re<sub>c</sub></th><th><i>k</i><sub>c</sub></th></tr></thead><tbody><tr><td>20</td><td>132.492</td><td>3.1270</td></tr><tr><td>30</td><td>131.989</td><td>3.1280</td></tr><tr><td>40</td><td>131.822</td><td>3.1283</td></tr><tr><td>50</td><td>131.746</td><td>3.1285</td></tr><tr><td>60</td><td>131.705</td><td>3.1286</td></tr><tr><td>100</td><td>131.647</td><td>3.1287</td></tr></tbody></table>

[see Eqs. (B1) and (B2)] unlike in a viscoelastic context. We consider here rigid boundary conditions

for the velocity perturbation, i.e., $\tilde{v} = 0$ at $r = r_i$ and $r = r_o$. Thus, the generalized eigenvalue

problem given by the latter setup of Eqs. (37)–(39) can be straightforwardly written in matrix form

$$ \sigma \begin{bmatrix} I_v & 0 & 0 \\ 0 & I_\lambda & 0 \\ 0 & 0 & 0 \end{bmatrix} \begin{bmatrix} V \\ \Lambda \\ P \end{bmatrix} = \begin{bmatrix} L_{vv} & L_{v\lambda} & -G \\ L_{\lambda v} & L_{\lambda\lambda} & 0 \\ D & 0 & 0 \end{bmatrix} \begin{bmatrix} V \\ \Lambda \\ P \end{bmatrix}, \quad (40) $$

where $V$ is the vertical matrix of the values of the components of velocity $\tilde{v}$ at each inner point of the gap, $\Lambda$ is the vertical matrix of the values of $\tilde{\lambda}$ at each point of the mesh, including the inner and outer radii, and $P$ is the vertical matrix of the values of $\tilde{p}$ taken in the middle points of two successive nodes of the velocity and structural parameter mesh.

The linear problem (40) admits a number of infinite eigenvalues that is two times the number of degrees of freedom of the pressure. The infinite eigenvalues have to be removed because they correspond to nonzero divergence velocity fields.

## B. Convergence test and validation

In order to test the convergence of our numerical scheme and to validate our method, the critical Reynolds number $Re_c$ and the critical axial wave number $k_c$ are determined using different number of nodes $M$ in the gap. The results are given in Tables I and II. For Newtonian fluids, many works allow us to validate our results. In a recent work [8], a similar numerical method gave $Re_c = 131.66$ and $k_c = 3.130$ in Newtonian fluids with $\eta = 0.9$. Those values are in very good agreement with ours. In addition, for Bingham fluids, Alibenyahia et al. [3] found $Re_c = 127.74943$ and $k_c = 3.183706$ with a spectral method at $Bn = 1$ and $\eta = 0.5$. Once again, our results in Table II agree with these values within an error below 0.1%. According to Tables I and II and the results of [3,8], we can estimate the relative error for the critical values of the Reynolds number $Re_c$ and the axial wave number $k_c$ below 0.1% when $M \ge 50$. Thus, we use $M = 50$ in the following.

<table><caption>TABLE II. Critical Reynolds number $Re_c$ and critical axial wave number $k_c$ for a Bingham fluid vs the number of nodes $M$ in the gap with a radii ratio $\eta = 0.5$.</caption><thead><tr><th><i>M</i></th><th>Re<sub>c</sub></th><th><i>k</i><sub>c</sub></th></tr></thead><tbody><tr><td>20</td><td>128.472</td><td>3.1695</td></tr><tr><td>30</td><td>128.057</td><td>3.1776</td></tr><tr><td>40</td><td>127.919</td><td>3.1803</td></tr><tr><td>50</td><td>127.857</td><td>3.1816</td></tr><tr><td>60</td><td>127.823</td><td>3.1822</td></tr><tr><td>100</td><td>127.776</td><td>3.1832</td></tr></tbody></table>
<figure><img src="image_7.png" /><figcaption>FIG. 6. Critical Reynolds number (a) Re_c and (b) Re_c0 = y_o Re_c and critical axial wave number (c) k_c and (d) k_c0 = y_o k_c vs b*. Large gap η = 0.5, Bn = 2, n_c = 1, ΔK* = τ_1* = 1, and a* = 1. The vertical dashed line stands for the critical value of b* = 1 where the strain rate at the interface ˙γ_o becomes nonzero.</figcaption></figure>

### C. Stability analysis of Couette flow of thixotropic yield stress fluids

To determine the critical eigenmode of the linear setup of Eq. (40), the algorithm seeks the minimum of the critical value of the Reynolds number Re_c depending on the wave number k for a given azimuthal mode n. The critical Reynolds number is reached when the real part of the eigenvalue σ is zero. The minimal value of Re_c is reached at the critical wave number k_c. We have verified that the critical perturbation is always axisymmetric, i.e., n = 0, by computing the critical Reynolds number for the azimuthal modes n from 0 to 3 within the range of our parameters for the thixotropic yielded fluids. The result is that the Taylor vortices, steady and axisymmetric, correspond always to the most unstable eigenmode of (40). Thus, the critical mode is not oscillating or three dimensional, just as expected for nonthixotropic fluids such as those of Bingham [2]. This suggests that the structural parameter λ plays a passive role. Indeed, the operator L_{λλ} in Eq. (40) is

$$L_{\lambda\lambda} = -\mathrm{diag}\left(a^* + b^* \dot{\gamma}_b^m + \frac{inV_b}{r}\right), \quad (41)$$
<figure><img src="image_8.png" /><figcaption>FIG. 7. (a)–(d) Streamlines of the critical velocity field of the perturbation. (e)–(g) Contour plot of the perturbation of the structural parameter $\lambda$. The contour lines step is 5% of the normalized amplitude of the perturbation. Blue stands for negative values and red for positive values (counterclockwise and clockwise spin for Taylor vortices). The large gap $\eta = 0.5$, $B_n = 2$, $n_c = 1$, $\Delta K^* = \tau_1^* = 1$, $a^* = 1$, and (a) $b^* = 0$, (b) and (e) $b^* = 1$, (c) and (f) $b^* = 5$, and (d) and (g) $b^* = 10$.</figcaption></figure>

where diag(·) stands for the diagonal matrix generated by the value of the argument at each node of the mesh. Equation (41) shows that it generates eigenvalues $\sigma$ such that their real part is always negative. This means that the perturbations of $\lambda$ vanishes without any coupling terms with the perturbation of the velocity. This confirms the passive role of Eq. (38). Only the values of the critical wave number $k_c$ and the critical Reynolds number $Re_c$ are modified compared to the Newtonian or shear-thinning cases. One can notice this result because it means that the unsteady effects of the thixotropy in cylindrical Couette flow that might occur above the threshold of the primary instability are nonlinear. Nevertheless, shear-banded flows may occur with thixotropic yielded fluids. This is a real difference from simple yield stress fluids.

The ratio $b^*/a^*$ denotes the resistance against the strain rate $\dot{\gamma}$ of the structure described by $\lambda$. The higher $b^*/a^*$ is, the easier the inner structure of the fluid is broken down by the shear. As shown previously in Fig. 5, the flowing region decreases because the yield stress collapses with the breakdown of the structure. If the fluid would be a viscous Newtonian fluid, the critical Reynolds number would increase because the gap becomes small, stabilizing the flow. In Fig. 6(a) the variation of $Re_c$ with $b^*$ suggests that our choice for the reference viscosity $\mu_{\text{ref}}$ is representative of an equivalent Newtonian fluid and thus we retrieve the stabilizing effect of the reduction of the gap width. Moreover, as the viscosity of the fluid decreases when the inner structure is broken, the fluid is stronger with stronger shear thinning. Thus, it is not surprising that the growth of $b^*/a^*$ ends by stabilizing the flow [Fig. 6(a)], as observed experimentally for a large gap by Escudier et al. [6] and shown by Alibenhahia et al. [3] when the shear-thinning index $n_c < 0.6$ for $\eta = 0.5$. Nevertheless, it has been known that the yielded flow region width is the relevant length scale since the earliest studies on the linear stability of yield stress fluids, e.g., Frigaard et al. [31]. The Reynolds number $Re_o = y_o Re$ is calculated taking into account the gap width of the yielded region $y_o$. The collapse of the width $y_o$ of the flowing zone (Fig. 5) is stronger than the stabilizing effect of shear thinning. Thus, the critical Reynolds number $Re_{c0}$ decreases until $y_o$ is close to its minimum value corresponding to a fully unstructured fluid [Fig. 6(b)]. The critical wave number $k_c$ mainly follows the evolution of the fluid gap width and thus it increases with $b^*$ [Fig. 6(c)]. By recalculating the wave
TAYLOR-COUETTE INSTABILITY IN THIXOTROPIC ...

<figure><img src="image_9.png" /><figcaption>FIG. 8. Critical Reynolds number (a) $Re_c$ and (b) $Re_{co} = y_o Re_c$ and critical axial wave number (c) $k_c$ and (d) $k_{co} = y_o k_c$. The large gap $\eta = 0.5$, $B_n = 2$, $n_c = 1$, $a^* = 1$, and $b^* = 1$. The vertical dashed line stands for the critical value of $\Delta K^* = 1$ where the strain rate at the interface $\dot{\gamma}_o$ becomes nonzero.</figcaption></figure>

345 number considering the fluid gap as the effective gap, we show that the wavelength $L_z = 2\pi/k_{co}$ first increases ($k_{co}$ decreases) [Fig. 6(d)] because the vortices fill better the flowing region [Figs. 7(a)–7(d)] when $b^*$ increases. It can be characterized by the width between the contour line corresponding to the 346 5% level of the streamlines and the outer limit of the flowing region [Figs. 7(a)–7(d)]. This zone is the 347 so-called dead zone [Fig. 7(a)]. The vortices are squeezed toward the inner wall [Fig. 7(a)] because of 348 the stratification of the viscosity. When the shear-banding appears, the stratification of the viscosity 349 is weakened and the Taylor vortices thicken [Figs. 7(c) and 7(d)]. The wave number $k_{co}$ tends to an 350 asymptotic value that corresponds to the rolls of an equivalent nonstructured fluid ($\Delta K^* = \tau_1^* = 0$). 351 To compute the critical Reynolds number of an equivalent nonthixotropic fluid when the gap width 352 corresponds to $y_o$ for $b^* = 10$, we have to set $\eta = 0.7403$ and $B_n = 0.7018$. The critical Reynolds 353 number is $Re_c = 126.9870$ and the critical wave number $k_c = 2.9177$. These values have to be 354 compared to the $Re_{co} = 119.4021$ and $k_{co} = 2.8981$ found when $b^* = 10$. As expected, for high 355 values of $b^*$, this fits with an equivalent nonthixotropic (simple) fluid flowing in a smaller gap. 356 357

It is worth noticing that the perturbation of the structural parameter $\lambda$ corresponds to the convection 358 of the structure by the Taylor vortices. In Figs. 7(e)–7(g) the negative zone of the perturbation 359
<figure><img src="image_10.png" /><figcaption>FIG. 9. Critical Reynolds number (a) $Re_c$ and (b) $Re_{co} = y_o Re_c$ and critical axial wave number (c) $k_c$ and (d) $k_{co} = y_o k_c$. The large gap $\eta = 0.5$, $B_n = 2$, $\Delta K^* = \tau_1^* = 0.5$, $a^* = 1$, and $b^* = 1$.</figcaption></figure>

<footnote n="360">360 corresponds to the convection from the inner cylinder where the strain rate destroys the microstructure toward the outer cylinder. Thus, in our parameters range, the linear stability is driven by the flow that governs the perturbation of the microstructure.</footnote>

<footnote>361</footnote>

<footnote n="362">362
The parameters $\Delta K^*$ and $\tau_1^*$ stabilize the flow according to Fig. 8(a). Nevertheless, as previously, the material gap size is not the most relevant to define the critical Reynolds number. From this point of view, the Reynolds number $Re_{co}$ decreases and the flow is destabilized [Fig. 8(b)]. Indeed, the shear-thinning behavior is strengthened when $\Delta K^*$ and $\tau_1^*$ increase. The high-viscosity area near the interface between the flowing and solidlike regions, corresponding to the dead zone, collapses after the onset of the shear banding. This is responsible for the increase of $k_{co}$ in Fig. 8(d) with $\Delta K^*$ and $\tau_1^*$. 363
364
365
366
367
368</footnote>

<footnote n="369">369
As expected, the effect of $n_c$ seems to be either destabilizing or stabilizing depending on whether we track $Re_c$ [Fig. 9(a)] or $Re_{co}$ [Fig. 9(b)]. For the shear-localized flows, increasing the shear-thinning behavior with the parameters $b^*$, $\Delta K^*$, and $\tau_1^*$ or $n_c$ has similar effects.</footnote>

<footnote n="370">370
371
372
To conclude on the effect of the thixotropy on the linear stability of the Couette flow, the critical mode corresponds to the axisymmetric Taylor vortices that are also found for simple yield stress fluids, such as Bingham fluids [2]. The critical eigenvalue is real, as in simple fluids. The results found by Landry et al. [2] with Bingham fluids or Alibenyahia et al. [3] with shear-thinning fluids</footnote>

<footnote>373</footnote>

<footnote n="374">374
375</footnote>
<figure><img src="image_11.png" /><figcaption>FIG. 10. Critical Reynolds number (a) Re₀₇ = (1 + ΔK*λₐₐₐ)Re₇ and (b) Re₀₇₀ = y₀Re₀₇ vs b*. The large gap η = 0.5, Bn = Bn₀ = 2, ΔK* = τ₁* = 1, and a* = 1.</figcaption></figure>

376 are retrieved. The critical perturbation is driven by the inertial term, i.e., by the centrifugal force.

377 The yield stress and the shear-thinning behavior confine the rolls toward the inner cylinder in the

378 wide-gap case. In a Bingham fluid, Chen et al. [4] show that the optimal perturbation is also shifted

379 toward the inner cylinder in the wide-gap case. Thus, even during the transient growth preceding

380 the onset of the instability, only the inner zone of the gap is perturbed. The perturbation of the

381 structural parameter is driven by the convection of the material because of the Taylor vortices.

382 The stabilizing or destabilizing effect depends on the reference viscosity used for the definition of

383 the Reynolds number. Nevertheless, the key point is that increasing the shear-thinning behavior

384 reduces the width of the inner region where the viscosity is low. As would be the case if the material

385 gap size would be reduced, it stabilizes the flow. The shear thinning is driven not only by n₇ but also

386 by the thixotropic parameters, i.e., the ratio b*/a*, ΔK*, and τ₁*. Although the thixotropy does not

387 produce a qualitative modification of the linear stability of the flow in the case of shear localization,

388 it allows for shear banding.

389 V. REFERENCE VISCOSITY AND REFERENCE YIELD STRESS

390 Figures 6(a) and 10(a) show that the characteristic value chosen for the viscosity may dramatically

391 change the conclusion about the effect of the parameters on the critical value of the Reynolds

392 number. Nevertheless, the asymptotic behavior for large b* can also be retrieved from Fig. 10(b):

393 The critical Reynolds Re₀₇₀ = 130.2569 is close to the one of the equivalent case with Bingham fluid,

394 Re₇ = 126.9870. The viscosity µ₀ is also a good choice to interpret the results, but our reference

395 µₐₐₐ might be more relevant from a practical point of view. Moreover, it reproduces the stabilizing

396 effect of thinning the gap, which would be observed with Newtonian fluids. Finally, we defined the

397 wall Reynolds number as

$$
\begin{equation}
\mathrm{Re}_w = \mathrm{Re}/\mu_w, \tag{42}
\end{equation}
$$

398 where µw is the shear viscosity of the fluid on the inner cylinder. This viscosity is relevant because

399 it fixes the resistive torque on the rotating cylinder, which is measured in classical rheological

400 experiments. Moreover, the centrifugal instability at the origin of the onset of the Taylor vortices

401 is triggered in the low-viscosity region, close to the inner cylinder. Figure 11(a) shows that the

402 inner-wall shear viscosity µw decreases for b* = 0.6 just before the onset of the shear banding.

403 For the shear-banded flows, the inner-wall shear viscosity increases to reach a limit value when b*

404 becomes high, i.e., when the structure is broken down even when the strain rate is low. The critical
<figure><img src="image_12.png" /><figcaption>FIG. 11. (a) Inner wall shear-viscosity $\mu_w$ and (b) critical Reynolds number $Re_{cw} = Re_c/\mu_w$ scaled with $\mu_w$ vs $b^*$. The large gap $\eta = 0.5$, $Bn = Bn_0 = 2$, $\Delta K^* = \tau_1^* = 1$, and $a^* = 1$.</figcaption></figure>

Reynolds number calculated with the wall viscosity is strongly growing before the onset of the shear banding but it is slightly constant (~250) for the shear-banded flow [Fig. 11(b)]. This observation suggests that the inner-wall shear viscosity is more relevant for the onset of the Taylor vortices when the velocity profile of the base flow corresponds to the shear banding, i.e., when the strain rate does not approach zero and the viscosity values are finite and moderate.

Thus, the thixotropic yield stress fluids behave mainly as a viscous fluid when the structure is fragile and the shear banding appears. In this case, the inner-wall viscosity is a good reference to predict the onset of the Taylor vortices.

## VI. CONCLUSION

In this work we have studied the base flow and the linear stability in a Couette cell of a thixotropic yield stress material modeled by Houška's model. This model with a single structural parameter allows for nonmonotonic composite flow curves depending on the ratio between the building and the breakdown parameters $b^*/a^*$ (Fig. 2). Nonmonotonic composite curves are known to trigger shear banding [25]. In shear banding, the structural parameter $\lambda$ jumps abruptly from a value below 1 to 1 across the interface between the fluid and the solidlike zones. The shear rate $\dot{\gamma}$ exhibits a discontinuity across the interface and the width of the high-viscosity zone, called the dead zone, collapses when shear banding occurs. In shear-banded flows, the selected stress at the interface is the yield stress of the solid material where $\lambda = 1$ if there is no diffusive term for the structural parameter. For the flows of simple yield stress fluids, only a smooth transition between the flowing and static regions is possible and the so-called dead zone is thicker than for thixotropic yield stress fluids.

The primary instability of the Couette flow is studied for a large gap ($\eta = 0.5$) when the Bingham number $Bn$ is sufficiently high to have a solidlike region in the gap. The thixotropy does not modify the kind of linear unstable mode. It remains steady and axisymmetric in the large range of parameters explored in comparison with simple yield stress fluids [2]. The choice of the reference viscosity in the definition of the Reynolds number plays a critical role in determining the stabilizing or the destabilizing effect of the thixotropic and the shear-thinning behaviors. By taking our reference viscosity $\mu_{ref}$ or the inner-wall shear-viscosity $\mu_w$, the thixotropy stabilizes the Couette flow in both cases because it increases the stratification of the viscosity over the gap. This stabilizing effect was also found in Ref. [3] with shear-thinning fluids. Nevertheless, rescaling the Reynolds number with the flowing gap size $y_o$ shows that there is a competition between the reduction of the stratification of the viscosity that destabilizes the flow and the reduction of the effective gap size because of the
breakdown of the fluid structure. The shear banding dramatically reduces the stratification of the viscosity, which remains finite and moderate in the fluid region. The flow is squeezed near the inner wall, which becomes similar to the one of a circular Couette setup with a smaller gap and without any dead zone. Thus, the y_o-scaled Reynolds number Re_{co} is only slightly growing and tends towards a constant value for high values of the breakdown parameter b*. The unsteady effects of the thixotropy may thus appear at the secondary instability of the Taylor vortices or it could be due to nonlinear effects, which should be studied in future works.

## 443 APPENDIX A: DETAILS OF THE LINEAR OPERATORS FOR THE VELOCITY PERTURBATION

In cylindrical coordinates, the perturbation vector of the velocity is written

$$ \bar{\mathbf{v}} = u\mathbf{e}_r + v\mathbf{e}_\theta + w\mathbf{e}_z. \quad (A1) $$

The expressions of the linear operators that appear in Eqs. (37)–(39) in the cylindrical coordinates system are

$$ -\overline{\nabla} \mathbf{v}_b \cdot \bar{\mathbf{v}} - \overline{\nabla} \bar{\mathbf{v}} \cdot \mathbf{v}_b = \frac{V_b}{r} \left( \frac{2}{r} v - i n u \right) \mathbf{e}_r - \left[ \left( \frac{V_b}{r} + \frac{\partial V_b}{\partial r} \right) u + \frac{i n V_b}{r} v \right] \mathbf{e}_\theta - \frac{i n V_b}{r} w \mathbf{e}_z \quad (A2) $$

for the inertial term and

$$ -\mathbf{v}_b \cdot \nabla \tilde{\lambda} - \bar{\mathbf{v}} \cdot \nabla \lambda_b = - \frac{\partial \lambda_b}{\partial r} u - \frac{i n V_b}{r} \tilde{\lambda} \quad (A3) $$

for the convective term of the structural parameter. The stress terms are

$$ \left. \frac{\partial \overline{\mathbf{\tau}}}{\partial \dot{\mathbf{\gamma}}_{ij}} \right|_{b} \dot{\mathbf{\gamma}}_{ij}(\bar{\mathbf{v}}) = \left[ \frac{1}{2}(\mu_1 - \mu_b)(\delta_{ir}\delta_{j\theta} + \delta_{i\theta}\delta_{jr})(\mathbf{e}_r \otimes \mathbf{e}_\theta + \mathbf{e}_\theta \otimes \mathbf{e}_r) + \mu_b(\mathbf{e}_i \otimes \mathbf{e}_j) \right] \dot{\mathbf{\gamma}}_{ij}(\bar{\mathbf{v}}), \quad (A4) $$

$$ \left. \frac{\partial \overline{\mathbf{\tau}}}{\partial \lambda} \right|_{b} \tilde{\lambda} = \tau_2 \lambda (\mathbf{e}_r \otimes \mathbf{e}_\theta + \mathbf{e}_\theta \otimes \mathbf{e}_r), \quad (A5) $$

with

$$ \mu_b = \left( \frac{1 + \Delta K^* \lambda_b}{1 + \Delta K^* \lambda_{\text{ref}}} \right) \dot{\gamma}_b^{n_c - 1} + \frac{\text{Bn}}{\dot{\gamma}_b} \left( \frac{1 + \tau_1^* \lambda_b}{1 + \tau_1^* \lambda_{\text{ref}}} \right), \quad (A6) $$

$$ \mu_1 = n_c \left( \frac{1 + \Delta K^* \lambda_b}{1 + \Delta K^* \lambda_{\text{ref}}} \right) \dot{\gamma}_b^{n_c - 1}, \quad (A7) $$

$$ \tau_2 = \left( \frac{\Delta K^*}{1 + \Delta K^* \lambda_{\text{ref}}} \dot{\gamma}_b^{n_c - 1} + \frac{\text{Bn}}{\dot{\gamma}_b} \frac{\tau_1^*}{1 + \tau_1^* \lambda_{\text{ref}}} \right) \dot{\gamma}_{r\theta,b}. \quad (A8) $$

The divergence of the stress tensor in cylindrical coordinates is

$$ \operatorname{div} \left( \left. \frac{\partial \overline{\mathbf{\tau}}}{\partial \dot{\mathbf{\gamma}}_{ij}} \right|_{b} \dot{\mathbf{\gamma}}_{ij}(\bar{\mathbf{v}}) \right) = \left[ 2\mu_b \frac{\partial^2 u}{\partial r^2} + 2\left(\frac{\mu_b}{r} + \frac{\partial \mu_b}{\partial r}\right) \frac{\partial u}{\partial r} - \left(n^2 \frac{\mu_1}{r^2} + k^2 \mu_b + \frac{2\mu_b}{r^2}\right) u \right. \\ + i n \frac{\mu_1}{r} \frac{\partial v}{\partial r} - i n \left(\frac{2\mu_b + \mu_1}{r^2}\right) v + i k \mu_b \frac{\partial w}{\partial r} \left] \mathbf{e}_r \right. \\ + \left[ i n \frac{\mu_1}{r} \frac{\partial u}{\partial r} + i n \left(\frac{2\mu_b + \mu_1}{r^2} + \frac{1}{r} \frac{\partial \mu_1}{\partial r}\right) u \right. \\ + \mu_1 \frac{\partial^2 v}{\partial r^2} + \left(\frac{\mu_1}{r} + \frac{\partial \mu_1}{\partial r}\right) \frac{\partial v}{\partial r} \right] $$
$$ 
\begin{aligned}
& - \left( \frac{2n^2\mu_b + \mu_1}{r^2} + k^2\mu_b + \frac{1}{r}\frac{\partial\mu_1}{\partial r} \right) v - nk\frac{\mu_b}{r}w \right] \mathbf{e}_{\theta} \\
& + \left[ ik\mu_b\frac{\partial u}{\partial r} + ik\left(\frac{\mu_b}{r} + \frac{\partial\mu_b}{\partial r}\right)u - nk\frac{\mu_b}{r}v \right. \\
& \left. + \mu_b\frac{\partial^2 w}{\partial r^2} + \left(\frac{\mu_b}{r} + \frac{\partial\mu_b}{\partial r}\right)\frac{\partial w}{\partial r} - \left(n^2\frac{\mu_b}{r^2} + 2k^2\mu_b\right)w \right] \mathbf{e}_z
\end{aligned} 
\quad (A9) $$

451 and

$$ \operatorname{div}\left(\frac{\partial \bar{\tau}}{\partial \lambda}\Big|_{b} \tilde{\lambda}\right) = in \frac{\tau_2}{r} \tilde{\lambda} \mathbf{e}_r + \left[ \tau_2 \frac{\partial \tilde{\lambda}}{\partial r} + \left( \frac{2\tau_2}{r} + \frac{\partial \tau_2}{\partial r} \right) \tilde{\lambda} \right] \mathbf{e}_{\theta}. \quad (A10) $$

452 Finally, the coupling term with the velocity in Eq. (38) is

$$ \frac{\partial \dot{\gamma}}{\partial \dot{\gamma}_{ij}}\Big|_{b} \dot{\gamma}_{ij}(\bar{\mathbf{v}}) = \frac{in}{r}u + \frac{\partial v}{\partial r} - \frac{v}{r}. \quad (A11) $$

453

## APPENDIX B: DEFORMATION OF THE YIELD STRESS BOUNDARY

For the Bingham-like problem ($n_c = 1$), the displacement of the interface is fully driven by the perturbation of the strain rate by an expression similar to Eq. (3.11) of Landry et al. [2]. The details of the mathematical proofs are given in Refs. [2,31]. They are based on the continuity of the stress across the interface of the plug. The displacement of the interface $\tilde{r}_o$ for a perturbation of the stress $\tilde{r}_o$ at the interface is

459 If $n_c = m = 1$, the perturbation of the stress at $r = r_o$ is linked to the perturbation of the strain rate

460 where $\tilde{K}$, $\Delta\tilde{K}$, and $\tilde{\tau}_1$ are defined by Eqs. (19) and (20). We recall that $\sigma$ is an eigenvalue of the linear setup of Eq. (40). Near the critical point, the real part $|\text{Re}(\sigma)| \ll 1$. As $\sigma$ is found to be a real number, we have $\sigma = 0$ at the critical point.

463 For $n_c < 1$, the linear perturbation of the strain rate $\delta\dot{\gamma}$ on the interface is negligible because $|\delta\dot{\gamma}| \propto |\tilde{\tau}_o|^{1/n_c}$ at $r = r_o$ [we only consider the case $\dot{\gamma}_b(r_o) = 0$]. This leads to no linear displacement of the yield stress boundary. This is compliant with the so-called dead zone, where the rolls vanish near the yield stress boundary.

[1] Z. Li and R. E. Khayat, A non-linear dynamical system approach to finite amplitude Taylor-Vortex flow of shear-thinning fluids, Int. J. Numer. Methods Fluids 45, 321 (2004).

[2] M. P. Landry, I. A. Frigaard, and D. M. Martinez, Stability and instability of Taylor-Couette flows of a Bingham fluid, J. Fluid Mech. 560, 321 (2006).

[3] B. Alibenyahia, C. Lemaitre, C. Nouar, and N. Ait-Messaoudene, Revisiting the stability of circular Couette flow of shear-thinning fluids, J. Non-Newtonian Fluid Mech. 183, 37 (2012).

[4] C. Chen, Z.-H. Wan, and W.-G. Zhang, Transient growth in Taylor-Couette flow of a Bingham fluid, Phys. Rev. E 91, 043202 (2015).

[5] G. I. Taylor, Stability of a viscous liquid contained between two rotating cylinders, Philos. Trans. R. Soc. London Ser. A 223, 289 (1923).

[6] M. P. Escudier, I. W. Gouldson, and D. M. Jones, Taylor vortices in Newtonian and shear-thinning liquids, Proc. R. Soc. London A 449, 1935 (1995).
[7] M. Naimi, R. Devienne, and M. Lebouche, Etude dynamique et thermique de l'écoulement de Couette-Taylor-Poiseuille; cas d'un fluide présentant un seuil d'écoulement, Int. J. Heat Mass Transfer **33**, 381 (1990).

[8] M. Pourjafar, E. Chaparian, and K. Sadeghy, Taylor-Couette instability of thixotropic fluids, Meccanica **50**, 1451 (2015).

[9] M. Houška, Engineering aspects of the rheology of thixotropic liquids, Ph.D. thesis, Czech Technical University, Prague, 1981.

[10] J. Šesták, R. Žitný, and M. Houška, Simple rheological models of food liquids for process design and quality assessment, J. Food Eng. **2**, 35 (1983).

[11] R. Žitný, A. A. Landfeld, J. Skočilas, J. Štancl, V. Flegl, M. Zachariášová, M. Jírů, and M. Houška, Hydraulic characteristic of collagen, Czech J. Food Sci. **34**, 479 (2015).

[12] A. Ahmadpour and K. Sadeghy, An exact solution for laminar, unidirectional flow of Houska thixotropic fluids in a circular pipe, J. Non-Newtonian Fluid Mech. **194**, 23 (2013).

[13] A. Wachs, G. Vinay, and I. A. Frigaard, A 1.5 D numerical model for the start up of weakly compressible flow of a viscoplastic and thixotropic fluid in pipelines, J. Non-Newtonian Fluid Mech. **159**, 81 (2009).

[14] R. Mendes, G. Vinay, G. Ovarlez, and P. Coussot, Modeling the rheological behavior of waxy crude oils as a function of flow and temperature history, J. Rheol. **59**, 703 (2015).

[15] J. Billingham and J. W. J. Ferguson, Laminar, unidirectional flow of a thixotropic fluid in a circular pipe, J. Non-Newtonian Fluid Mech. **47**, 21 (1993).

[16] C. D. Andereck, S. S. Liu, and H. L. Swinney, Flow regimes in a circular Couette system with independently rotating cylinders, J. Fluid Mech. **164**, 155 (1986).

[17] P. D. Olmsted, O. Radulescu, and C. Y. D. Lu, Johnson-Segalman model with a diffusion term in cylindrical Couette flow, J. Rheol. **44**, 257 (2000).

[18] C. Y. D. Lu, P. D. Olmsted, and R. C. Ball, Effects of Nonlocal Stress on the Determination of Shear Banding Flow, Phys. Rev. Lett. **84**, 642 (2000).

[19] M. A. Fardin, O. Radulescu, A. Morozov, O. Cardoso, J. Browaeys, and S. Lerouge, Stress diffusion in shear banding wormlike micelles, J. Rheol. **59**, 1335 (2015).

[20] O. Radulescu, P. D. Olmsted, J. P. Decruppe, S. Lerouge, J. F. Berret, and G. Porte, Time scales in shear banding of wormlike micelles, Europhys. Lett. **62**, 230 (2003).

[21] M. Jenny and B. Nsom, Primary instability of a Taylor-Couette flow with a radial stratification and radial buoyancy, Phys. Fluids **19**, 108104 (2007).

[22] A. M. Philippe, C. Baravian, M. Jenny, F. Meneau, and L. J. Michot, Taylor-Couette Instability in Anisotropic Clay Suspensions Measured Using Small-Angle X-Ray Scattering, Phys. Rev. Lett. **108**, 254501 (2012).

[23] P. D. Olmsted, Perspectives on shear banding in complex fluids, Rheol. Acta **47**, 283 (2008).

[24] P. C. F. Møller, S. Rodts, M. A. J. Michels, and D. Bonn, Shear banding and yield stress in soft glassy materials, Phys. Rev. E **77**, 041507 (2008).

[25] T. Divoux, M. A. Fardin, S. Manneville, and S. Lerouge, Shear banding of complex fluids, Ann. Rev. Fluid Mech. **48**, 81 (2016).

[26] H. J. Wilson and S. M. Fielding, Linear instability of planar shear banded flow of both diffusive and non-diffusive Johnson-Segalman fluids, J. Non-Newtonian Fluid Mech. **138**, 181 (2006).

[27] C. Perge, M. A. Fardin, and S. Manneville, Surfactant micelles: Model systems for flow instabilities of complex fluids, Eur. Phys. J. E **37**, 1 (2014).

[28] G. Ovarlez, S. Rodts, X. Chateau, and P. Coussot, Phenomenology and physical origin of shear localization and shear banding in complex fluids, Rheol. Acta **48**, 831 (2009).

[29] G. Ovarlez, S. Cohen-Addad, K. Krishan, J. Goyon, and P. Coussot, On the existence of a simple yield stress fluid behavior, J. Non-Newtonian Fluid Mech. **193**, 68 (2013).

[30] S. Jarny, N. Roussel, S. Rodts, F. Bertrand, R. Le Roy, and P. Coussot, Rheological behavior of cement pastes from MRI velocimetry, Cement Concrete Res. **35**, 1873 (2005).

[31] I. A. Frigaard, S. D. Howison, and I. J. Sobey, On the stability of Poiseuille flow of a Bingham fluid, J. Fluid Mech. **263**, 133 (1994).