Shear-banding and Taylor-Couette instability in thixotropic yield stress fluids

In the present work, we study the flow of thixotropic yield stress fluids between two concentric cylinders. In order to take into account the thixotropy, the constitutive relation uses a structural parameter which is driven by a kinetic equation. Here, the Houska's model is considered. Depending on the breakdown rate of the structural parameter, localization or shear-banding are observed. We show that for fragile structures, a shear-banding flow may be observed although for stronger structures, only localisation of the flow is observed such as in Bingham fluids. Physical explanations of the shear-banding discussed by several authors in the literature highlight that the shear-banding may be associated with a discontinuity into the structure of the material and a non-monotonic evolution of the stress according to the constitutive relation with the strain rate. Solving numerically the flow, we show that such a rheological model based on the existence of a structural parameter is able to predict shear-banding. Moreover, the consequences of the thixotropy on the linear stability of the azimuthal flow is studied in a large range of parameters. Although the thixotropy allows shear banding for the base flow, it does not modify fundamentaly the stability of the Couette flow compared to a simple yield stress fluid. The apparent shear-thinning behaviour 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.


I. INTRODUCTION
Shear banding occurs in many complex materials which manifest non-linear rheological behavior such as micellar solutions, granular paste or colloidal gels [1]. It has been shown that the coexistence of static and flowing regions can be associated with the existence of a yield stress which is often related to the existence of a rigid network between the elements of the fluid which has to be broken for the system starts to flow. Shear banding occurs, in many cases, in systems which exhibit a competition between at least two mechanisms: a breakdown process due to an external applied shear and a spontaneous restructuring of the fluid due to the non-linear interaction between its elements [2]. This feature suggests that the onset of shear banding is related to a coupling between the constitutive law of the material and the evolution of its internal structure. However questions remain on the quantitative effects of such a structure parameter on the behavior of the flow and its stability. The structure parameter represents some physical properties of the instantaneous structure state of the fluid. The constitutive law depends on this structural parameter whose evolution is driven by a kinetic equation. In this picture, thixotropy models [3][4][5] such as Houska's model are relevant to predict flow instabilities and heterogeneities in complex structured fluids. Ovarlez et al. 2009 [6] and Coussot & Ovarlez 2010 [7] recently discuss the physical origin of shear localization and shear banding in complex fluids. Authors highlight the existence of a discontinuity of the shear rate profile during shear banding contrary to shear localization where the shear rates goes to zero continuously as one approaches the static region. In the last case, the shear rate is zero on both side of the static and flowing regions. "Shear localization" refers to a stress heterogeneity inherent to the geometry for a yield stress fluid whereas "shear banding" is directly related to the non-linear rheology of the material. The onset of shear banding or shear localization highly depends on the stress distribution. Analysing in depth these phenomenon needs a control of that distribution. In that context, Couette flows are particularly relevant to study these phenomena since the shear stress distribution is heterogeneous but well controlled contrary to cone or plate geometry. Taylor-Couette flow is often used in rheology as a reference shear flow. Moreover, since the historical work of Taylor [8], it is a paradigm for studies of stabilities and transition to turbulence. With Newtonian fluids, the transitional regimes observed in Taylor-Couette flows have been widely studied [9] when the inner or the outer cylinders rotates. In the present work, we only consider the case where the inner cylinder rotates at a given angular velocity ω i such as the velocity at the inner radius r i is v i = ω i r i ( fig. 1). 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 [8][9][10][11]. One can notice that an elastic instability of the flow may occurs even at very low velocity [12] for non-Newtonanian fluids. According to the studies of the hydrodynamic stability of shear-thinning [11,13] and Bingham fluids [11,14,15] in Couette flow, 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. For a Bingham fluid, the dependence of the wavelength of the Taylor vortices to the yield stress changes when a non-yielded zone appears close to the outer cylinder. Chen et al. [15] have shown that the effects of yield stress on the energy transient growth and flow structure of the optimal perturbations are quite different for the wide or narrow-gap cases. Considering axisymmetric optimal modes for corotating cylinders, the peak of the amplitude of optimal perturbation is found to be shifted toward the inner cylinder with increasing yield stress for the wide-gap case whereas the peak is shifted toward the outer cylinder for the narrow-gap case. Three dimensional perturbations were also investigated in [11] but it was found that the most unstable modes are axisymmetric, even for very strong shear-thinning fluids.
Certainly, Taylor vortices do occur in yield stress fluids. This is evidenced by the numerical experiments of Lockett et al. (1992) [16], Coronado-Matutti, Souza Mendes & Carvalho (2004) [17] and Jeng & Zhu (2010) [18], by phenomenological evidence from applications such as oil well drilling (i.e. observed changes in frictional pressure), and also by only few direct experimental studies. Nouar, Devienne & Lebouche (1987) [19] and Naimi, Devienne & Lebouche (1990) [20] have studied axial flow through an annulus with rotating inner cylinder, using CMC and Carbopol solutions, respectively. The former behaves as a power-law fluid and the latter has a yield stress, as well as shear-thinning behaviour. In both cases only axisymmetric Taylor vortices are reported, with their appearance retarded by the presence of an axial flow. Naimi et al. (1990) [20] reports that the yield stress appears to stabilize the flow. Fardin et al. [21] show experimental evidences of Taylor-like vortices in shear-banding flow of giant micelles. This observation suggests that the Taylor-vortices can commonly occur in a very large range of thixotropic or non-thixotropic non-Newtonian fluids.
In this article, we propose a numerical study of the stationary flow of a thixotropic yield stress fluid in a Couette configuration. For modelling steady state flows, we use the Houska's model where the constitutive law depends on a structural parameter driven by a kinetic equation.It is commonly admitted in the litterature [22] that a diffusion term is recquired to model shear banding. We show, in this article, that the stress selection at the interface by the structural diffusion term admits a limit value equals to the yield stress of the fully structured fluid when the structural diffusion coefficient tends to zero. We show that the inherent thixotropy of the model controls the transition from a solution where a static and a flowing regions coexists in which the shear rate profile is continuous (shear localization) to a solution where the two regions coexist but in which the shear rate profile is discontinuous (shear banding). We show that the onset of shear banding is controled by the competition between restructuring effects and breakdown effects due to the flow. Although the thixotropy initiates shear banding in the base flow, a linear analysis of the stability of flow solutions shows that the nature of the linear unstable mode which is steady and axisymmetric in a large range of explored parameters does not depend on the thixotropic character of the flow.

A. Modified Houska's model
The Houska's model [23,24] is built from the Hershel-Bulkley model, commonly used for non-elastic yield stress fluids, considering that the consistency K and the yield stress τ 0 depend linearly on the structural parameter λ. Thus, for a non-zero strain rate, i. e.γ = 0, the stress tensorτ is given by : The parameters ∆K and τ 1 determine the sensitivity of the consistency and the yield stress respectively with the structural parameter λ. n c is the shear thinning index. The second invariant of the strain rate tensorγ is defined bŷ using the Einstein summation convention where the elements of the strain tensorγ ij = ∂v i /∂x j + ∂v j /∂x i are defined with the components of the fluid velocityv. The indices i and j stand for the cylindrical coordinates r, θ and z ( fig.  1) and the hat. denotes a dimensional variable. According Olmsted et al. [22], a stress diffusion has to be added to keep the unicity of the steady solution in case of shear-banding by selecting the stress at the interface between the bands. Fardin et al. [25] show that the length based on this stress diffusion is at the order of magnitude of the molecular size in wormlike micelles. Thus, the structural parameter λ is determined by the kinetic equation : D is the structural diffusion coefficient. The diffusion term in eq. (3) has been added to the original version of Houska's model to take into account the stress diffusion, via the structural diffusion, in our model. a and b are respectively the building and the breakdown parameters. The thixotropic breakdown index m is taken equal to 1 in the following. The values of the structural parameter are within the range 0 ≤ λ ≤ 1. The value λ = 1 means that the fluid is fully structured and, at the opposite, λ = 0 means that it is fully unstructured. The evolving dynamics of the stress is governed by the equation (3).

B. Nondimensional equations
To non-dimensionalize 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: where ρ is the density of the fluid. For the non-Newtonian fluids, several choices may be done for the reference viscosity. The choice of the reference viscosity will be discussed in the section V. The reference viscosity is the plastic viscosity of the fluid at the given reference strain rate v i /d. When the strain rate is v i /d, the corresponding structural parameter λ ref is given by eq. (3) at equilibrium: where the nondimensional building and breakdown parameters are respectively. Thus, the reference viscosity built using the reference strain rate reads : where µ 0 = K(v i /d) nc−1 can be recognized as the standard reference viscosity of a power law fluid. ∆K = ∆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 depends on the solvent. The reference viscosity µ ref depends on the ratio of the breakdown parameter b over the building parameter a . The reference viscosity decreases when the ratio 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 incompressible fluids are : where v =v/v i stands for the reduced velocity and p =p/(ρv 2 i ) for the reduced pressure. One can notice thatp is the modified pressure including the hydrostatic pressure. The Reynolds number is defined using the reference viscosity (7) by: Thus, the reduced stress tensor reads: whereγ andγ are the nondimensional strain rate and strain tensor, τ 1 = τ 1 /τ 0 is the reduced thixotropic yield stress. The equation (11) involves the Bingham number which is the ratio between the yield stress and the plastic viscous stress: is the standard Bingham number of a Hershel-Bulkley fluid. For a strong yield stress behavior, the Bingham number is high and localized flows are expected. According the equation (11), the reduced yield stress is: For the structural parameter, the non-dimensional version of the equation (3) is: with the non-dimensional structural diffusion coefficient C. Boundary conditions for the flow The inner and outer radii of the Couette setup are now defined by with η = r i /r e the radii ratio. The velocity vector v is written in the cylindrical basis as v = v r e r + v θ e θ + v z e z . We only consider the case where the inner cylinder rotates and the outer cylinder is fixed. Thus, the boundary conditions are: • At the inner radius r = r i , the velocity components are v θ = 1 and v r = v z = 0.
• At the outer radius of the flowing zone r = r o , the velocity components are v r = v θ = v z = 0.
• 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 τ (r e ) ≥ τ y , r o = r e , else, τ (r o ) = τ y .
At r = r o , the yield stress τ y given by eq. (13) becomes with λ o the structural parameter at the interface. Thus, the last boundary condition becomes τ (r o ) = τ yo and the stress τ yo at the interface is defined by eq. (18). Nevertheless, the stress condition at the interface between the fluid and solid-like zone is well defined only if the structural parameter λ is continuous across the interface. When the shear banding appears, the structural parameter λ is sharply discontinuous across the interface when no diffusive gradient term is added (D = 0) in the governing equation of λ and thus the yield stress can not be well defined (Olmsted et al. [22] and Lu et al. [26] in 2000). Authors showed 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/numerical noise history because it will select arbitrary 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 has been recently interpret as a non-local effect at the molecular scale in flow of micellar suspensions in [25]. As the value of the stress diffusion coefficient, similar to the structural diffusion coefficient, is found to be very small [25,27], we will focus on the cases where D tends to zero. In this case, we found a limit value for the stress at the interface in case of shear-banding and steady flows. The limit is the yield stress of the fully structured material τ ys (Fig. 2). Thus, we can solve the steady equations with D = 0, setting the stress at the interface : D. NUMERICAL METHOD FOR TRANSIENT FLOWS

Arbitrary Lagrangian Eulerian method
To perform the numerical solution, we use a finite difference method for the spatial discretisation. The mesh is split in two regions: the first one is the fluid zone for radii r i ≤ r ≤ r o and the second one is the solid zone for radii r o ≤ r ≤ r e . The mesh points are regularly spread between the inner radius r i and the outer radius r o in the first region and between r o and the external radius r e in the second region. The number of points are the same in the two regions. In order to solve the transient flow, we use an Arbitrary Lagrangian-Eulerian (ALE) method [28,29] for the interface tracking. The method consists in solving the governing equations in the gap considering a mesh where one node (1D problem) sticks to the interface. To keep a constant radial length between two successive nodes of the mesh at each time step, the arbitrary velocity of nodes is given by: where e r is the radial unit vector. u o is the radial velocity of the interface. The set of equations (8, 9 and 14) has to be written in the moving domain [see 28,29]: In the last equation (23), the convective gradient term, u∇λ ensures the continuity of the structural parameter and the yield stress across the interface as long as the interface velocity is not zero, even if there is no diffusive term in the initial equation (14), i.e. when D tends to zero, as shown in the figure 3. Thus, in the case Re = 0 and u = 0 (finite time), the stress at the interface given by eq. (18) is defined even if D tends to zero in shear banded flows because of the smoothing due to the motion of the interface. The steady state of the shear banded flow, i. e. with a sharp discontinuity, can only be considered as an asymptotic state when the velocity tends to zero at t → +∞. One can remark that the figure 3-b shows that the stress value at the interface tends to the yield stress of the fully structured material τ ys as found previously with D tending to zero in the steady flows.
To compute the flow in the fluid region r i ≤ r ≤ r o , one solves the set of equations (21, 22 and 23) with the boundary conditions given in the section II C. The stress at the interface is equal to the local yield stress. This latter condition provides the interface velocity u o . The structural parameter is also computed in the solid-like region consideringγ = 0. One can notice that when r o = r e , there is no interface between a fluid and a solid-like region. Thus, the stress has just to be above the local yield stress, i. e.γ ≥ 0 everywhere.

Time discretisation
The time scheme used for the velocity is a semi-implicit time scheme at the order one. Here, we focus on the 1D equations where considering v = V (t, r)e θ and u = U (t, r)e r , U being defined by eq. (20): The index n denotes the time step t n and V n , for instance, stands for V (t n , r n ). The radii r n and r n+1 are related by r n+1 = r n + U n ∆t with the time step defined as ∆t = t n+1 − t n . The stress τ = −τ rθ is given by the equation (11). To implicit the time scheme, one writes the developpment for the stress at the order one in time: whereγ = −γ rθ = V /r − ∂V /∂r. For numerical stability reasons, the time discretization scheme used λ n (r n+1 ) is obtained by the interpolated value of the former field λ n at the current node position r n+1 . Thus, the equation (23) is rewritten as: To test that the final steady state does not depend on the starting condition, as it would be expect when D = 0 [22], we define two cases. The first one is a start at rest with the material fully structured, i. e. V 0 = 0 and λ 0 = 1 everywhere. As the thickness of the fluid region collapses at the initial time, we take an arbitrary profile for the first non-zero velocity at t 1 = ∆t. The velocity V 1 is defined with a second order polynomial function of the radius. The velocity field V 1 verifies all the boundary conditions of the velocity and the strain rate is zero at the interface r o,1 . The stress τ 1 in the initial fluid region is calculated in order to balance the pulse momentum provided to accelerate the flow from 0 to the arbitrary velocity V 1 (eq. 24). Thus, the stress depends on U 0 , i. e. on u o,1 = (r o,1 − r i )/∆t. The position r o,1 of the interface is calculated to equalize the stress at the interface and the yield stress. From this time, the velocity profile is computed solving numerically the equations (24), (25) and (27). In the second scenario, the calculation starts from a fully unstructured and flowing material in the whole gap. The figures 4 show that both transient stage tends to the same steady state even if the asymptotic steady state is a shear banded flow as it will be discussed in the following. The asymptotic case corresponds to the steady state where the interfacial stress is the yield stress of the fully structured material. Setting D = 0, in the right-hand-side of the equation (14) would select an intermediary value for the stress at the interface in a shear banded flow. This stress would be smaller than the yield stress of the fully structured material and depend on the scale length defined by the diffusive coefficient. Indeed, it would smooth the variation of λ between its values in the flowing region near the interface (i. e. < 1) and 1 in the far distance, at the scale of the diffusive length, in the solid-like region.
In the following, we will focus on the case where the diffusive coefficient D is 0. Knowing that the asymptotic steady state of the transient stage corresponds to the steady state where the stress at the interface τ o = τ (r o ) = τ ys defined by the equation (19), i. e. eq. (18) with λ o = 1, whatever shear banding or not, one can solve directly the steady equations.

E. NUMERICAL METHODS FOR STEADY FLOWS
Only the fluid region needs to be considered to solve the steady flow. One uses the same finite difference method for the spatial discretisation as for the unsteady problem. The stress points are taken between two successive velocity points to insure the numerical accuracy of the scheme for the velocity. For the derivative operations, the standard second order centred scheme is used. The numerical method used for the spatial discretisation is quite well established and is similar to the ones used, for instance, in [30][31][32]. Moreover, a validation and a convergence test are performed in the subsection IV B.
To calculate the base flow, we consider the steady axisymetric solution of the equations (8), (9) and (14), i. e. v b = V b (r)e θ and λ = λ b (r). In the fluid domain, i. e. r i ≤ r ≤ r o , the only non-zero element of the strain rate tensor isγ rθ . The strain rate does not reach zero in the flowing region and its sign is always negative. Thus, one can write: The only non-zero element of the stress tensor τ is then τ rθ . Considering the previous assumptions for the flow, the azimuthal component of the equation (8) becomes: leading to the well known result for the steady Couette flow: 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: If r o ≥ r e according the equation (31), all the material in the gap flows and r o = r e . In the next section, the flow curves show that the minimal value τ min of the stress may be below τ yo . Thus, for τ min re 2 ≤ C < τ yo r 2 e , an alternative to the equation (31) is to set r o = r e . In practice, it means that if there is no interface at the initial state, the fluid region fits the whole gap for τ min re 2 ≤ C < τ yo r 2 e and if there is a solid-like region in the initial state, the flowing region is confined between r i and r o < r e according Eq. (31).
To compute the flow velocity in the yielded region, we calculate the strain rate in the yielded region by solving the regular setup of equations at each point of the mesh: with τ yb the yield stress given by eq. (13) replacing λ by λ b . Once we obtain the strain rateγ b and the structural parameter λ b for a given constant C by solving the setup of equations (32)(33), the linear equation (28) is solved to calculate the fluid velocity with V b (r o ) = 0 as boundary condition. Finally, one has to find the value of C such as V b (r i ) = 1 using the algorithm available in Matlab to calculate the zero of a real non-linear function. If needed, the pressure P b of the base flow can be obtained by integrating the equation: and setting the inner pressure P b (r i ) = 0 for instance.

III. EFFECT OF THE THIXOTROPY ON THE BASE FLOW
In the next sections, we set the thixotropic index breakdown to m = 1. It seems reasonable to argue that the structural parameter λ modifies the viscous term and the yield stress with the same order of magnitude. Thus, we fix ∆K = τ 1 and so Bn = Bn 0 . As we focus our study on the cases where both a flowing and a solid-like regions exist, i. e. r o < r e , the Bingham number is fixed to Bn = 2. Finally, for steady state flow 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 nodes in the flowing region of the gap to insure a good accuracy for the linear stability analysis as shown in the validation subsection IV B.
We report in Fig.5 the evolution of the composite flow curves under controlled shear rate for different values of the breakdown parameter b . The composite curves are obtained straightfully by replacing the structural parameter λ by its relation toγ (32) in the constitutive law (11). As b increases, the composite curve drops from a monotonic to a non-monotonic behavior which presents an unstable branch leading to shear banding [33]. This result indicates that, in the range of parameters studied here, Houska's model is able to predict shear banding.
The existence of a critical shear rateγ c in thixotropic yield stress fluids is often explained in terms of an underlying decreasing branch of the flow curve at low shear rates [2]. In such a scenario, the constitutive relation of the material is a decreasing function of the strain rate between 0 andγ 0 [see 1, for instance]. In the case of the shear localization, the constitutive law is a growing function ofγ and its derivative is strictly positive. Thus, the shear-banding may appear if the sign of the derivative changes at a critical strain rateγ 0 > 0. In other words, the necessary condition to allow the shear-banded flow is: ∃γ 0 ≥ 0, such as ∂τ ∂γ = 0 Thus, the range of strain rate [0,γ 0 ], where ∂τ /∂γ < 0, is unstable and the shear-banding allows to avoid the forbidden values ofγ. In the steady-state, the derivative of the stress given by the Houska's model is: When eq. (35) admits at least one solution with the function (36), the shear-banding appears if the stress values reach levels that correspond to multiple possible strain rates according the steady constitutive function of the material. For Bingham-like behavior where m = n c = 1, eq. (35) using the derivative of τ given by eq. (36) admits only one positive rootγ 0 :γ The derivative of τ (36) is negative forγ ∈ [0,γ 0 [. In case of the localisation of the flow in r ∈ [r i , r o ], r o < r e , the boundary condition at r = r o for the stress is: The latter equation (39) . 8). Under controlled torque, the flow may freeze in a Couette with a small gap when the radii ratio η > τ min /τ yo andγ o <γ 0 . This result is consistent with the multiple yield stresses observed by Ovarlez et al. [6]. One can notice that if n c < 1 and m ≥ 1, and thus there is a range of strain rate values close to zero where τ (γ) is a growing function ofγ even if there is a range of positive strain rates where the stress τ decreases whenγ grows. This case would be similar to the flow curve of a semidilute wormlike micelle solution with a yield stress like in the Fig. 1b of [34]. 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 solid-like region, i. e. at r = r o . Thus, in the following, we limit our parametric study to shear banded flows with a banding interface between the flowing and static regions only (flow curve like in Fig. 1c of [34]), i. e. with n c = 1 and shear localization with n c ≤ 1.

B. Velocity profiles and structure parameter
The strain rate, the viscosity the velocity and the structure parameter profiles are shown in Fig. 9. We see that whatever the value of b , a flowing and a static region coexists 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 structure 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 [6,35]. In that case, the shear localization is inherent to the existence of the yield stress and no shear banding is observed. By contrast, for larger values of b , the shear rate becomes non-zero at the outer boundary of the flowing region and the structural parameter λ does not reach 1 as it stands in the solid-like region when the breakdown parameter is greater than a critical value b c = 1 (figs. 9). Discontinuous strain rate profiles between a static and a flowing region have been observed using MRI measurements in cement pastes [36] and bentonite suspensions [35]. 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 eq. (38) and figs. 10, increasing the parameter τ 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. 9d and 10d), the more the viscosity drops significantly and rapidly. A band in terms of structure takes place next to the unyielded region ( fig. 9b and  10b). Then, the transition from shear localization to shear banding is feasible.
At last, we explore the effect of the shear-thinning index n c . When n c < 1, the simple shear-banded flows as observed in figs. 9 or 10 are not observed because the constitutive relation of the material is always a growing function forγ sufficiently close to zero (eq. 41), allowing small values forγ in the flow . In that case the flow is always shear-localized. This contrast with the previous cases discussed above (n c = 1) where small values ofγ fall in the unstable branch of the flow curve and then lead to shear banded flows. In order to keep only one flowing zone in the explored range of n c , we set ∆K = τ 1 = 0.5. Staying on the continuous solution for the structural parameter λ, the inner strain rate increases with the shear-thinning index n c and the structural parameter λ reaches smoothly the fully structured state of the solid-like zone because of the growing viscosity near y o (figs. 11). As it would be expected from the velocity profiles obtained with Carreau fluids by Alibenyahia et al. [11], the flow is confined close to the inner cylinder when the shear-thinning index n c decreases. It confirms that the shear-thinning behaviour confined Detail of the black box   fig. 12-a). Quite similar remark can be done when ∆K and τ 1 increase. Indeed, increasing b , ∆K or τ 1 makes the shear-thinning behaviour 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. 12-c). Nevertheless, the model with one structural parameter allows to predict the shear banding which can not be described by the Bingham law or the Hershel-Bulkley law. When the thixotropic parameters are increasing above some critical values given by eq. (38), the shear-banding appears smoothly from the shear localization where the strain rate is continuous between the flowing and static regions.
Now, we will study the linear stability of the base flow, whatever we have shear-banding or not.

A. Equations setup
To perform a linear analysis of stability, the fluid velocity, the structural parameter and the pressure are decomposed such as: v,λ andp describe the perturbation of the base flow considering the azimuthal mode n and the axial wave number k. Injecting eqs. (42 -44) in the general setup of equations (8,14,9) and after withdrawing the nonlinear terms, the linear setup of equations for the perturbation of the base flow is: The indices ij stand for r, θ or z and Einstein's convention for summation is used. Thus, the generalized eigenvalues problem given by the latter setup of equations (45-47) can be straightfully written in matrix form: V is the vertical matrix of the values of the components of velocityṽ at each inner point of the gap. Λ is the vertical matrix of the values ofλ at each point of the mesh, including the inner and outer radii. P is the the vertical matrix of the values ofp taken in the middle points of two successive nodes of the velocity and structural parameter mesh. The linear problem (48) admits a number of infinite eigenvalues which is two times the number of degree of freedom of the pressure. The infinite eigenvalues have to be eliminate because they correspond to non-zero divergence velocity fields. Nevertheless, it is not harmful because Matlab's algorithm generates true infinite eigenvalues avoiding any ambiguity.

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 the tables I and II. For Newtonian fluids, many works allows us to validate our results. In a recent work [32], a similar numerical method gives Re c = 131.66 and k c = 3.130 in Newtonian fluids with η = 0.9. Their values are in very good agreement with ours. More, for Bingham fluids, Alibenyahia et al. [11] found Re c = 127.74943 and k c = 3.183706 with a spectral method at Bn = 1 and η = 0.5. Once again, our results in the table II agree this values within an error below 0.1%.
According, the tables I and II and the results of [11] and [32], 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 ≥ 50. Thus, we use M = 50 in the following.

C. Stability analysis of Couette flow of thixotropic yield stress fluids
To determine the critical eigenmode of the linear setup of equations (48), the algorithm seeks the minimum of the critical value of the Reynolds number Re c depending on the wave number k for a given azimutal 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 axisymetric, 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 results is that the Taylor vortices, steady and axisymmetric, corresponds always to the most unstable eigenmode of (48). One can observed that the equation (46) does not generate oscillating or three-dimensional modes and only the values of the critical wave number k c and the critical Reynolds number Re c are modified compared to the Newtonien or shear-thinning cases. One can notice this result because it means that the unsteady effects of the thixotropy in cylindrical Couette flow which 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 with simple yielded fluids. The ratio b /a denotes the resistance against the strain rateγ of the structure described by λ. The higher b /a is, the easier the inner structure of the fluid is broken down by the shear. As shown previously in fig. 12, 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 increases because of the gap becomes small. In the fig. 13-a, the variation of Re c with b suggests that our choice for the reference viscosity µ 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 decrease when the inner structure is broken, the fluid is stronger and stronger shear-thinning. Thus, it is not surprising that the grow of b /a ends by stabilizing the flow ( fig. 13-a) as it is observed experimentally for large gap by Escudier et al. [10] and shown by Alibenhahia et al. [11] when the shear-thinning index n c < 0.6 for η = 0.5. Nevertheless, for yielded fluids, one can argue that the effective gap corresponds to the fluid zone. The Reynolds number Re o = y o Re is calculated with the gap width of the yielded region y o . In this case, the collapse of the width y o of the fluid zone ( fig. 12) because of the breakdown of the inner micro-structure is stronger than the stabilizing effect of shear-thinning and the critical Reynolds number Re co decreases until y o is close to its minimum value corresponding to a fully destructured fluid ( fig. 13-b). The critical wave number k c mainly follows the evolution of the fluid gap width and thus it increases with b ( fig. 13-c). Of course, by recalculating the wave number considering the fluid gap as the effective gap, show that the wave length L z = 2π/k co first increase because the high viscosity near the end of It is worthy to notice that the perturbation of the structural parameter λ corresponds to the convection of the structure by the Taylor vortices. In the figs. 14-(e-g), the negative zone of the perturbation corresponds to the convection from the inner cylinder where the structural parameter is low toward the outer cylinder. Thus, in our parameter range, the linear stability is driven by the Navier-Stokes equation and the perturbation of the structural parameter is a passive mode. The parameters ∆K and τ 1 stabilise the flow according the figure 15-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. 15-b). Indeed, the decrease of the viscosity with λ is greater when ∆K and τ 1 increase. The collapse of the dead zone due to high viscosity near the transition between the fluid and solid-like zones after the onset of the shear-banding is responsible of the increasing of k co in fig. 15-d.
Finally, the effect of n c seems to be either destabilizing or stabilizing depending if we track Re c ( fig. 16-a) or Re co ( fig. 16-b). The figures 17 show that the Taylor-vortices are shifted toward the inner cylinder as the shear-thinning behavior increases. The effect is similar as increasing the shear-thinning behavior with the parameters b , ∆K and τ 1 ) for the shear-localized flows. In shear-banded flows (figs 14c-d), the width of the dead zone collapses and thus, the Taylor-vortices expand in the whole fluid region.
To conclude on the effect of the thixotropy on the linear stability of the Couette flow, first the critical mode corresponds to the axisymmetric Taylor-vortices which are also found for simple yield stress fluids, such as Bingham fluids [14]. The critical eigenvalue is real as in simple fluids. The results found by Landry et al. [14] with Bingham fluids or Alibenyahia et al. [11] with shear-thinning fluids are retrieved. The critical perturbation is driven by the inertial term, i. e. by the centrifugal force. The yield stress and the shear thinning behavior confine the rolls toward the inner cylinder in the wide-gap case. In a Bingham fluid, Chen et al. [15] show that the optimal perturbation is also shifted toward the inner cylinder in the wide-gap case. Thus, even during the transient growth preceding the onset of the instability, only the inner zone of the gap is perturbed. The perturbation of the structural parameter is driven by the convection of the material because of the radial and axial velocity of the Taylor vortices. The stabilizing or destabilizing effect depends on the reference viscosity used for the definition of the Reynolds number. Nevertheless, the key point is that increasing the shear-thinning behaviour reduces the width of the inner region where the viscosity is low. As it would be the case if the material gap size would be reduced, it stabilizes the flow. The shear-thinning is not only driven by n c but also by the thixotropic parameters, i. e. the ratio b /a , ∆K and τ 1 . Although the thixotropy does not produced a qualitative modification of the linear stability of the flow in case of shear localization, it can make appear a second kind of base flow. In the shear localization case, the strain rate drops smoothly to zero at the interface between the fluid and solid-like regions. Thus, there is a so-called dead zone where the viscosity becomes very strong near the interface. The width of the dead zone collapses in the shear banded flows (Figs. 9b, 10b). Thus, the Taylor vortices can occupied the fluid gap as it would be the case with shear-thinning fluids in an equivalent ratio of radii η eq = r o /r i , i. e. in a narrow-gap case. In every cases, the width of the fluid gap y o is the most relevant reference length to describe the onset of the linear instability, as it would be the case for simple yield stress fluids, such as Bingham's fluids.

V. REFERENCE VISCOSITY AND REFERENCE YIELD STRESS
The figure 19-a shows that the reference chosen for the viscosity may dramatically change the conclusion about the effect of the parameters on the critical value of the Reynolds number. Nevertheless, the comparison between figs. 13-b and 19-b shows that the relevant length is the width of the fluid gap y o . Moreover, the asymptotic behaviour for large b can also be retrieved in fig. 19-b: the critical Reynolds Re co = 130.2569 close to the one of the equivalent case with Bingham fluid, Re c = 126.9870.
The viscosity µ 0 is also a good choice to interpret the results but our reference µ ref might be more relevant in a  practical point of view. Moreover, it reproduces the stabilizing effect of thinning the gap which would be observed with Newtonian fluids. Finally, we defined the wall Reynolds number as where µ w is the shear-viscosity of the fluid on the inner cylinder. This viscosity is relevant because it determines the resistive torque on the rotating cylinder which is measured in rheological experiment. Moreover, the centrifugal instability at the origin of the onset of the Taylor vortices is triggered at the low viscosity region close to the inner cylinder. The fig. 20-a shows that the inner wall shear-viscosity µ w decreases for b = 0.6 just before the onset of the shear-banding. For the shear-banded flows, the inner wall shear-viscosity increases to reach a limit value when b becomes high, i. e. when the structure is broken down even when the strain rate is low. The critical 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. 20-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 fall down zero and the viscosity values are finite and moderate. Thus, the thixotropic yielded fluids behave mainly as a viscous fluid when the structure is fragile and the shearbanding 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 modelized by the Houska's model. This model with one structural parameter allows non-monotonic composite flow curves depending on the ratio between the building and breakdown parameters b /a (Fig. 5). In case of nonmonotonic composite curves, a shear banded flow occurs if only a part of the gap flows. In shear banding, the structural parameter λ jumps abruptly from a value < 1 to 1 across the interface between the fluid and the solid-like zones. The shear rateγ is sharply discontinuous across the interface and the width of the high viscosity zone, called the dead zone, collapses with the shear-banding. In shear banded flows, the transient flow shows that the selected stress at the interface is the yield stress of the solid material where λ = 1 if there is no diffusive term for the structural parameter. One can notice that adding a diffusive term would select a lower stress corresponding to the yield stress of the partially structured material. For the following linear stability analysis, we only consider the case where the diffusive parameter is neglectable. The primary instability of the Couette flow is studied for a large gap (η = 0.5) when the Bingham number Bn is sufficiently high to have a solid-like region in the gap. The thixotropy does not modify the kind of the linear unstable mode which is still steady and axisymmetric in the large range of parameters explored in comparison with simple yield stress fluids [14]. We retrieve that the choice of the viscosity to scaled the Reynolds number may modify the conclusion about the effect of the thixotropy as for shear-thinning fluids. Considering a Reynolds number based on our reference viscosity µ ref or the inner-wall shear-viscosity µ w , lead to the conclusion that the thixotropy stabilizes the Couette flow because it increases the stratification of the viscosity over the gap. This stabilizing effect was also found by [11] with shear-thinning fluids. Nevertheless, re-scaling 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 which 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 but it squeezes the fluid area near the inner wall. Thus, the y o -scaled Reynolds number Re co is only slightly growing and tends to a constant value for high values of the breakdown parameter b .
Finally, The thixotropy allows for shear banding flows but does not induce transition to unsteady and nonaxisymmetric flows according to linear perturbations. For shear localized flows, the transition is driven by both the stratification of the viscosity and the size of the flowing region as in Bingham's fluids for instance [14]. When shear banding occurs, the stratification of the viscosity vanishes and the transition is mainly controlled by the effective gap (y o ). The shear-banding modifies the linear stability of the flow which becomes similar to the one of a circular Couette setup with a smaller gap and without any dead zone. The unsteady effects of the thixotropy may thus appears ∂γ ∂γ ij bγ ij (ṽ) = in r u + ∂v ∂r − v r (A11)