Ibrahim Cheddadi
email: ibrahim.cheddadi@inria.fr

François Pierre Saramito

F Graner

Steady Couette flows of elastoviscoplastic fluids are non-unique

Keywords: elastoviscoplastic / viscoelastoplastic fluids, non-Newtonian fluids, Couette experiment, liquid foam, suspensions, mathematical modelling

Steady Couette flows of elastoviscoplastic fluids are non-unique

Introduction

Localisation is a phenomenon often observed in two-or three-dimensional shear flows of complex materials: [START_REF] Coussot | Coexistence of liquid and solid phases in flowing soft-glassy materials[END_REF] observed it for emulsions, Salmon et al. (2003a) for colloids, and [START_REF] Howell | Stress fluctuations in a 2D granular Couette experiment: a continuous transition[END_REF]; [START_REF] Mueth | Signatures of granular microstructure in dense shear flows[END_REF]; [START_REF] Losert | Particle dynamics in sheared granular matter[END_REF]; [START_REF] Huang | Flow of wet granular materials[END_REF] for wet granular materials. It consists of a coexistence between a region localised near a moving boundary, where the material flows like a liquid, and another region where the material behaves like a solid.

Since the pioneering experiment of [START_REF] Debrégeas | Deformation and flow of a two-dimensional foam under continuous shear[END_REF] (Fig. 1), liquid foams (gas bubbles dispersed within a continuous liquid phase, as explained by [START_REF] Weaire | The physics of foam[END_REF] and [START_REF] Cantat | Les mousses: structure et dynamique[END_REF]) have been widely used for experimental, theoretical and numerical studies of localisation (for reviews see e.g. [START_REF] Höhler | Rheology of liquid foam[END_REF]; Schall and van Hecke (2010); [START_REF] Barry | Nonlocal effects in the continuum theory of shear localisation in 2d foams[END_REF]). Their discrete units, the gas bubbles, are easy to observe (especially in two dimensions) and to manipulate. Moreover, they display simultaneous elastic, viscous, plastic behaviours (referred to as elastoviscoplastic, or EVP), thus covering a wide range of behaviours observed in many complex materials.

The aim of this paper is to show that including tensorial elasticity in the classical viscoplastic Herschel-Bulkley (VP) model leads to many improvements in the understanding of Couette flows of non-thixotropic EVP fluids such as emulsions, liquid foams, or carbopol gel. These materials exhibit normal stresses that arise from the local anisotropy (hence the necessity of a tensorial description) of the elasticity related to their micro-structure. Localisation can appear if the material yields, that is, if the material is plastic; in the regions below the yield strain, the normal stresses can remain finite even in a steady-state flow. If, in addition, viscous dissipation occurs during plastic events, elasticity is coupled to viscosity in the flowing region, so that the normal stresses are coupled to the velocity gradient, even in the Couette geometry. [START_REF] Cheddadi | Numerical modelling of foam Couette flows[END_REF][START_REF] Cheddadi | Erratum to numerical modelling of foam Couette flows[END_REF]Cheddadi et al. ( , 2011a) ) have previously explored this approach with the Saramito ( 2007) model (Bingham-like plastic dissipation), for cylindrical Couette flows of liquid foams and other EVP flows around an obstacle; they have successfully explained the observations of normal stresses components measured by [START_REF] Janiaud | Foam in a two-dimensional Couette shear: a local measurement of bubble deformation[END_REF] in the experimental data of [START_REF] Debrégeas | Deformation and flow of a two-dimensional foam under continuous shear[END_REF].

In the present work, the theoretical predictions of the [START_REF] Saramito | A new elastoviscoplastic model based on the Herschel-Bulkley viscoplasticity[END_REF] model are compared with experimental measurements, including shear and normal stresses when available. This model is an extension of the [START_REF] Saramito | A new constitutive equation for elastoviscoplastic fluid flows[END_REF] model, that includes a Herschel-Bulkley-like plastic dissipation.

We study the influence of the dimensionless rheological parameters, including the yield strain and the cylindrical Couette geometry curvature (introduced in section 3). For simplicity, we focus here on the low velocity regime corresponding to most published foam Couette flow experiments; we postpone, for future work, analysis of both the quasistatic regime where viscosity plays no visible role, studied in simulations [START_REF] Wyn | Simulations of two-dimensional foam rheology: localization in linear Couette flow and the interaction of settling discs[END_REF][START_REF] Raufaste | Discrete rearranging disordered patterns: prediction of elastic and plastic behavior, and application to two-dimensional foams[END_REF], and the high velocity regime where viscous and friction effects are dominant (see e.g. [START_REF] Katgert | Rate dependence and role of disorder in linearly sheared two-dimensional foams[END_REF][START_REF] Katgert | Flow in linearly sheared twodimensional foams: From bubble to bulk scale[END_REF][START_REF] Katgert | Couette flow of twodimensional foams[END_REF]). The inclusion of elasticity and of tensorial description are two essential features of the present modelling, leading to three predictions not captured by scalar and/or VP models. First, we show that normal stresses that depend on the preparation of the material [START_REF] Labiausse | Shear induced normal stress differences in aqueous foams[END_REF] can persist as residues even in steady flow. Second, we predict that velocity flow profiles are either smooth or non-smooth depending on the measured value of the stress tensor. Third, as a consequence of these two results, two-and three-dimensional cylindrical Couette flows of EVP materials are non-unique, even in steady state. To summarise, the effect we describe in this paper is not specific to a given material microstructure, but is more generally a consequence of the material's visco-elasto-plasticity and of the specifically tensorial nature of the Couette flow.

The outline of the paper is as follows. Section 2 reviews and discusses the main open questions found in the literature, which we address here. Section 3 discusses some constitutive equations for EVP materials, in particular the EVP model of [START_REF] Saramito | A new elastoviscoplastic model based on the Herschel-Bulkley viscoplasticity[END_REF]. Section 4 presents the solutions of this model and some of the main features which are absent from VP models: the effect of initial conditions, memory effects, non-uniqueness, and non-smooth solutions; we also explain how such EVP model can be compared to actual experimental data. Section 5 examines how variations of the EVP model parameters affect these flow features, and provides an approximate analytical expression for the localisation length (eq. 15). Section 6 discusses the consequences of our results. Section 7 summarises our findings.

Open questions

During the last ten years there has been extensive debate on the Couette flow of various complex fluids, raising some theoretical questions (see eg. Schall and van Hecke (2010)):

• What is the physical origin of localisation?

• Where does the material localise?

• Why do some experiments report smooth profiles at the localisation position and some others experiments report non-smooth ones?

Before we examine the status of these questions, we first need to clarify the vocabulary and hypotheses used in the literature:

• We assume here that the materials of interest can be described using continuum mechanics: this implies that there is a representative volume element (RVE). The RVE should be smaller than the scale of the global flow, but large enough so that one can define: (i) variables such as stresses, strains, velocity gradient; and (ii) parameters, such as elastic modulus, yield strain and viscosity. Note that in an EVP material stress, strain and deformation rate are variables which can be independently defined, and measured in situ if the material can be imaged in full details [START_REF] Graner | Discrete rearranging disordered patterns, part I: robust statistical tools in two or three dimensions[END_REF].

• Since the phrase "shear banding" [START_REF] Vermant | Large-scale structures in sheared colloidal dispersions[END_REF] is used to describe different things, we prefer not to use it. We instead use "smooth" (resp: "non-smooth") to indicate that the velocity gradient is continuous (resp: discontinuous).

• The word "localisation" is a historical term used consistently to mean "coexistence of rigid and flowing regions". We thus use it here.

• In cylindrical Couette flows the localisation radius, denoted as r c , admits several possible definitions (see [START_REF] Gilbreth | Flow transition in two-dimensional foams[END_REF] or section 7.3 of [START_REF] Weaire | The continuum theory of shear localization in two-dimensional foam[END_REF]). Here, we choose to define r c as the position separating the flowing region from the rigid one, or equivalently as the limiting value of the radius at which the deformation rate is zero.

• There are at least two different ways, one in mathematics and one in rheology, used to define the norm of a tensor, and thus the von Mises criterion regarding the yield strain and stress: we choose the mathematical one, presented below (see eq. ( 4) and discussion thereafter).

What is the physical origin of the localisation?

There are several explanations for the origin of the localisation; many of them are reviewed and discussed by [START_REF] Schall | Shear bands in matter with granularity[END_REF]. [START_REF] Coussot | Coexistence of liquid and solid phases in flowing soft-glassy materials[END_REF] used MRI methods to measure the local velocity in three-dimensional cylindrical Couette flows of other EVP materials such as carbopol gel, and more generally yield stress fluids such as bentonite suspensions and cement paste. Despite the apparent simplicity of shear flows, a common description of these experiments is still lacking (see [START_REF] Ovarlez | Phenomenology and physical origin of shear localization and shear banding in complex fluids[END_REF] for a review): on the one hand, thixotropic materials exhibit an intrinsic critical shear rate, i.e. these materials cannot flow homogeneously at a shear rate smaller than a critical value, which is characteristic of the material; on the other hand, non-thixotropic materials may or may not exhibit a critical shear rate (that does not seem to be intrinsic to the material), as discussed below in section 2.3. We focus here on non-thixotropic materials and try to explain this peculiar behaviour.

Intrinsic explanations invoking a non-monotonic constitutive equation

In complex materials which display a shear-induced structural transition, a possible source of localisation is the coexistence of two different shear rates at the same stress, that is, the shear stress versus shear rate curve is multi-valued [START_REF] Huseby | Hypothesis on a certain flow instability in polymer melts[END_REF][START_REF] Berret | Isotropic-to-nematic transition in wormlike micelles under shear[END_REF][START_REF] Porte | Inhomogeneous flows of complex fluids: mechanical instability versus non-equilibrium phase transition[END_REF][START_REF] Decruppe | Insight in shear banding under transient flow[END_REF][START_REF] Bénito | Shear bands in a continuum model of foams: how a 3d homogeneous material may seem inhomogeneous in 2d[END_REF]. In fact, foam experiments [START_REF] Khan | Foam rheology: III. Measurement of shear flow properties[END_REF] and simulations [START_REF] Kabla | Quasistatic rheology of foams I. Oscillating strain[END_REF][START_REF] Okuzono | Intermittent flow behavior of random foams: a computer experiment on foam rheology[END_REF][START_REF] Raufaste | Discrete rearranging disordered patterns: prediction of elastic and plastic behavior, and application to two-dimensional foams[END_REF] suggest that in the quasi-static regime the shear stress versus shear strain curve passes through an overshoot before reaching a plateau, thus being multi-valued [START_REF] Clancy | The response of 2D foams to continuous applied shear in a Couette rheometer[END_REF][START_REF] Weaire | A simple analytical theory of localisation in 2d foam rheology[END_REF][START_REF] Weaire | The continuum theory of shear localization in two-dimensional foam[END_REF]).

This explanation is probably not generally applicable to foams in Couette experiments: foams usually contain at least a few percent liquid, so that their yield strain is lower than for ideally dry foams [START_REF] Marmottant | Discrete rearranging disordered patterns, part II: 2D plasticity, elasticity and flow of a foam[END_REF] and the overshoot becomes undectable [START_REF] Raufaste | Discrete rearranging disordered patterns: prediction of elastic and plastic behavior, and application to two-dimensional foams[END_REF]. Moreover, even in theory and simulations of dry foams, the overshoot disappears if the foam is too disordered [START_REF] Raufaste | Discrete rearranging disordered patterns: prediction of elastic and plastic behavior, and application to two-dimensional foams[END_REF]. This discussion is beyond the scope of the present paper, where we explore the complex behaviour emerging from a simple EVP constitutive equation.

Extrinsic explanations invoking a non-homogeneous stress

In the case of a cylindrical Couette geometry, the curvature induces a spatial heterogeneity of the stress, the inner part being above the yield stress and flowing, and the outer part being below the yield stress and non-flowing [START_REF] Adams | Nonmonotonic models are not necessary to obtain shear banding phenomena in entangled polymer solutions[END_REF]. This explanation of localisation requires a yield stress, and thus arises in VP and EVP models.

An alternative interpretation has been presented by [START_REF] Wang | Impact of boundaries on velocity profiles in bubble rafts[END_REF] for two-dimensional foams. They performed experiments with two configurations in a plane Couette geometry: either bubbles floating on water (bubble raft), or bubbles confined between water and a glass plate. They observed localisation only in the second case, when a glass plate is present. Therefore, the competition between the internal viscosity of the foam and the external friction from the glass plates has been suggested by [START_REF] Wang | Impact of boundaries on velocity profiles in bubble rafts[END_REF] and [START_REF] Janiaud | Two-dimensional foam rheology with viscous drag[END_REF] as a possible cause for localisation, even in the case of circular geometries [START_REF] Clancy | The response of 2D foams to continuous applied shear in a Couette rheometer[END_REF]. [START_REF] Cheddadi | Numerical modelling of foam Couette flows[END_REF] have reconciled these views by comparing in detail the Saramito (2007) EVP model with both steady and transient regimes of the [START_REF] Debrégeas | Deformation and flow of a two-dimensional foam under continuous shear[END_REF] experiment of a two-dimensional foam confined between two glass plates. They have shown that the simple explanation of a non-homogeneous stress is still valid for such systems. The curvature of flow lines and the friction on the glass plates are only two different possible causes for stress heterogeneity. In practice, the friction of glass plates is so small [START_REF] Cheddadi | Numerical modelling of foam Couette flows[END_REF] that it plays a visible role only in the absence of flow line curvature, as in the plane geometry of [START_REF] Wang | Impact of boundaries on velocity profiles in bubble rafts[END_REF]. In two-dimensional cylindrical Couette flows, the available data are not compatible with this friction-dominated interpretation, and are compatible with the stress heterogeneity that arises from geometry.

In what follows, we neglect the friction on plates, and use two-dimensional foams as model systems for two-or three-dimensional non-thixotropic EVP fluids.

What is the localisation position r c ?

In cylindrical Couette flows, [START_REF] Weaire | The continuum theory of shear localization in two-dimensional foam[END_REF] have reviewed the dependency of the localisation position on the cylinder velocity. They find that in the low velocity regime, it is basically undetermined: there is a range of possible values for the localisation position and, moreover, the size of this range diverges when the velocity decreases. There are two questions: is there a hidden variable which fixes the localisation position? for a given experiment, is the position predictable?

To answer these questions, in section 5 we explicitly determine the localisation position, and study its dependence on model parameters and initial conditions.

2.3 Why smooth and/or non-smooth profiles? The first observation of foam Couette flows by [START_REF] Debrégeas | Deformation and flow of a two-dimensional foam under continuous shear[END_REF] reported a smooth velocity profile (Fig. 2a).

However, in 2002, [START_REF] Coussot | Coexistence of liquid and solid phases in flowing soft-glassy materials[END_REF] observed discontinuous shear rate profiles on Couette flows of emulsions and suspensions (Fig. 2b). A discontinuity (denoted γc ) of the shear rate at r = r c was measured: so, the velocity profile is non-smooth at r = r c , since its derivative is related to the shear rate. Such non-smooth profiles were found by others between 2003 and 2008: for worm-like micelles by Salmon et al. (2003a), for lyotropic lamellar systems by Salmon et al. (2003b), for Couette foam flows by [START_REF] Lauridsen | Velocity profiles in slowly sheared bubble rafts[END_REF]; [START_REF] Gilbreth | Flow transition in two-dimensional foams[END_REF]; [START_REF] Dennin | Discontinuous jamming transitions in soft materials: coexistence of flowing and jammed states[END_REF]; [START_REF] Krishan | Viscous shear banding in foam[END_REF], and interpreted theoretically by [START_REF] Denkov | Jamming in sheared foams and emulsions, explained by critical instability of the films between neighboring bubbles and drops[END_REF]; [START_REF] Clancy | The response of 2D foams to continuous applied shear in a Couette rheometer[END_REF]; [START_REF] Weaire | The continuum theory of shear localization in two-dimensional foam[END_REF].

Surprisingly, in 2008-2010, experiments published by [START_REF] Katgert | Rate dependence and role of disorder in linearly sheared two-dimensional foams[END_REF][START_REF] Katgert | Flow in linearly sheared twodimensional foams: From bubble to bulk scale[END_REF][START_REF] Katgert | Couette flow of twodimensional foams[END_REF], Coussot andOvarlez (2010), and[START_REF] Ovarlez | Investigation of shear banding in threedimensional foams[END_REF] showed smooth velocity profiles and continuous shear rates at r = r c : γc equals zero. Note that some papers with contradictory results shared either an author [START_REF] Coussot | Coexistence of liquid and solid phases in flowing soft-glassy materials[END_REF][START_REF] Da Cruz | Viscosity bifurcation in granular materials, foams, and emulsion[END_REF][START_REF] Huang | Flow of wet granular materials[END_REF][START_REF] Coussot | Physical origin of shear-banding in jammed systems[END_REF] or a set-up, the bubble raft [START_REF] Lauridsen | Velocity profiles in slowly sheared bubble rafts[END_REF][START_REF] Gilbreth | Flow transition in two-dimensional foams[END_REF][START_REF] Dennin | Discontinuous jamming transitions in soft materials: coexistence of flowing and jammed states[END_REF][START_REF] Krishan | Viscous shear banding in foam[END_REF][START_REF] Katgert | Rate dependence and role of disorder in linearly sheared two-dimensional foams[END_REF][START_REF] Katgert | Flow in linearly sheared twodimensional foams: From bubble to bulk scale[END_REF][START_REF] Katgert | Couette flow of twodimensional foams[END_REF]. [START_REF] Coussot | Physical origin of shear-banding in jammed systems[END_REF] explained this discrepancy by questioning the quality of the experiments: "previous data on a specific foam [START_REF] Rodts | From discrete to continuum flow in foams[END_REF] were probably affected by experimental artefacts". Similarly, [START_REF] Ovarlez | Investigation of shear banding in threedimensional foams[END_REF] explain that "Our measurements demonstrate that three-dimensional foams do not exhibit observable signatures of [discontinuous] shear banding. This contrasts with the results of [START_REF] Rodts | From discrete to continuum flow in foams[END_REF] and da Cruz et al. ( 2002) which we have shown to pose several experimental problems" and mention that "the case of bubble rafts is still unclear".

Our results emphasise the intrinsic sensitivity of the equations to the sample history [START_REF] Bénito | Shear bands in a continuum model of foams: how a 3d homogeneous material may seem inhomogeneous in 2d[END_REF]; this could explain why experimental artifacts such as impurities [START_REF] Rodts | From discrete to continuum flow in foams[END_REF] or bubble rupture [START_REF] Da Cruz | Viscosity bifurcation in granular materials, foams, and emulsion[END_REF] could yield drastic changes in observations. Even when experimental artifacts are eliminated, our results regarding the effect of trapped stresses due to initial conditions provide a deep reason to explain why, according to the set-up or foam preparation, either a smooth or a non-smooth profile could appear.

3 Constitutive equation

Brief review of EVP models

A number of closely related models have appeared in the literature (see [START_REF] Saramito | A new constitutive equation for elastoviscoplastic fluid flows[END_REF] for a review). Table 1 distinguishes models with respect to their behaviour before and after yielding, their applicability to three-dimensional general flows, and to the existence of a proof that the dissipation is positive.

We now discuss the mathematical formulation and properties of some of these models. The equations are written in two dimensions with polar coordinates to fit with the cylindrical Couette geometry. [START_REF] Schwedoff | La rigidité des liquides[END_REF] rigid solid VE fluid [START_REF] Bingham | Fluidity and plasticity[END_REF] rigid solid Newtonian fluid X X [START_REF] Herschel | Measurement of consistency as applied to rubber-benzene solutions[END_REF] rigid solid power-law fluid X X [START_REF] Oldroyd | A rational formulation of the equations of plastic flow for a Bingham fluid[END_REF] elastic solid Newtonian fluid X X [START_REF] Isayev | Viscoelastic plastic constitutive equation for flow of particle filled polymers[END_REF] elastic solid VE fluid X X [START_REF] Doraiswamy | The Cox-Merz rule extended: a rheological model for concentrated suspensions and other materials with a yield stress[END_REF] elastic solid power-law fluid X Puzrin et al. ( 2003

) elastic solid VE solid X X Saramito (2007) VE solid VE fluid X X Bénito et al. (2008) VE solid VE fluid X X Saramito (2009) VE solid power-law VE fluid X X
Table 1: Summary of the referenced EVP models. The 3D column is marked when the model has been written in a general tensorial sense, e.g. with objective derivatives. The TH column is also marked when the model has a positive dissipation according to the second law of thermodynamics. The italic mark X in the 3D or the TH columns means that the corresponding property was derived after publication of this model. For instance, the Bingham and the Herschel-Bulkley models were first proposed in a one-dimensional simple shear flow context, then extended to three dimensions, and finally found to satisfy the second law of thermodynamics.

the constitutive equation is closed with equations for momentum and mass balance, and appropriate Couette flow boundary and initial conditions. Since external forces and inertia are negligible here (see section 2.1 and [START_REF] Cheddadi | Numerical modelling of foam Couette flows[END_REF]), the momentum balance reduces to

∇ • τ = 0, ( 1
)
where τ is the stress tensor. The material is assumed to be incompressible, so the mass balance is

∇ • v = 0, ( 2
)
where v is the velocity. model (see also [START_REF] Oldroyd | A rational formulation of the equations of plastic flow for a Bingham fluid[END_REF]):

VP model: Herschel-Bulkley

   τ = 2K|D| n-1 D + τ Y D |D| when |D| = 0, |τ d | ≤ τ Y otherwise,
or equivalently: max 0,

|τ d | -τ Y 2K|τ d | n 1 n τ = D, (3)
where D = (∇v + ∇v T )/2 is the deformation rate tensor, |D| = D ij D ij its Euclidian norm, and (∇v) ij = (∂ j v i ) is the velocity gradient tensor. The so-called deviatoric stress τ d is defined according to the spatial dimension in which the flow is investigated:

τ d = τ -1 N I
, where I is the identity matrix, and N = 1, 2 or 3. We focus here on two-dimensional flows and take N = 2. Here, τ Y > 0 is the yield stress, K > 0 is the consistency parameter, and n > 0 is the power-law index.

The von Mises criterion, |τ d | ≤ τ Y , involves the Euclidian norm of the deviatoric part of the stress. In cylindrical coordinates it is written:

|τ d | = 2τ 2 rθ + (τ rr -τ θθ ) 2 2 1 2 . ( 4
)
Note that, favouring the shear stress, some authors use a slightly different definition of the deviatoric stress norm [START_REF] Raufaste | Discrete rearranging disordered patterns: prediction of elastic and plastic behavior, and application to two-dimensional foams[END_REF]:

(τ 2 rθ + (τ rr -τ θθ ) 2 /4) 1 2
. This rheological definition would lead to an alternative, equivalent model: keeping the same equation( 3), it would only multiply the deviatoric stress norm |τ d | by 2 -1/2 , the yield stress τ Y by 2 -1/2 , and the consistency K by 2 (n-1)/2 .

In the circular Couette geometry, V denotes the velocity of the inner cylinder and ∆r = r e -r 0 the width of the gap. Note that η = K(V /∆r) n-1 has the dimension of a viscosity and that η = K when n = 1. The dimensionless Bingham number is

Bi = τ Y ∆r ηV . (5)
It compares the yield stress τ Y with a characteristic viscous stress ηV /∆r. Let us also introduce the dimensionless number Co:

Co = 1 - r 0 r e , (6)
that quantifies the curvature of the circular Couette geometry. When the two radii are close, this number is close to zero. Conversely, when the curvature is extreme, e.g. r 0 becomes small, this number tends to one.

The Herschel-Bulkley model predicts localisation in the circular Couette geometry, as a result of the stress heterogeneity. Its position r c can be numerically computed (see appendix A, equation ( 16)).

The Herschel-Bulkley model reduces to the Bingham one when n = 1, to a power-law fluid when τ Y = 0, and to a Newtonian fluid when n = 1 and τ Y = 0. The shear thinning behaviour is associated with 0 < n < 1 and the (less usual) shear thickening behaviour with n > 1.

3.3 VE model: Oldroyd [START_REF] Oldroyd | On the formulation of rheological equations of states[END_REF] proposed the following viscoelastic model:

η G ∇ τ +τ = 2ηD, (7)
where G > 0 is the elastic modulus and η/G is a relaxation time. The total stress is σ = 2η 2 D + τ where η 2 is a second viscosity, often called the solvent viscosity in the context of polymer solutions. When η 2 = 0, the Oldroyd model reduces to the so-called Maxwell model. The upper-convected tensorial derivative ∇ τ is defined by

∇ τ = ∂τ ∂t + v.∇τ -τ.∇v T -∇v.τ. (8)
The dimensionless Weissenberg number is

W e = ηV G ∆r . (9)
It compares the characteristic viscous stress ηV /∆r with the elastic modulus G. When 1/G = 0, the Oldroyd model reduces to a Newtonian one.

3.4 EVP model: Saramito [START_REF] Saramito | A new constitutive equation for elastoviscoplastic fluid flows[END_REF][START_REF] Saramito | A new elastoviscoplastic model based on the Herschel-Bulkley viscoplasticity[END_REF] and [START_REF] Bénito | An elasto-visco-plastic model for immortal foams or emulsions[END_REF] derived (independently) tensorial elastoviscoplastic models that combine viscoelastic and viscoplastic properties (see Fig. 3 and table 1).

0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 τ Y η G
They satisfy the second law of thermodynamics and match the behaviour of non-thixotropic EVP materials like foams and emulsions: elastic solid before yielding and viscoelastic flow after yielding.

The EVP model presented by [START_REF] Saramito | A new elastoviscoplastic model based on the Herschel-Bulkley viscoplasticity[END_REF] is simple enough to allow for the numerical resolution (with good convergence, see appendix B) of the associated partial differential equations even in intricate two- (Cheddadi et al., 2011a) and three-dimensional geometries, and is thus suitable for practical and industrial purposes. It is written:

1 2G ∇ τ + max 0, |τ d | -τ Y 2K|τ d | n 1 n τ = D. ( 10
)
When 1/G = 0 we obtain the Herschel-Bulkley model, eq. ( 3). Conversely, when n = 1 and τ Y = 0 we obtain the Oldroyd model, eq. ( 7). Finally, when both 1/G = 0, τ Y = 0 and n = 1, the model is Newtonian.

In addition to the independent dimensionless numbers (W e, Bi, Co, n) already introduced for the Herschel-Bulkley and Oldroyd models, we define the elastic yield strain ε Y as

ε Y = τ Y 2G = Bi W e 2 . ( 11
)
This dimensionless parameter is a measure of the softness and deformability of the material. It has been shown to be the main parameter for the characterisation of EVP materials (Cheddadi et al., 2011a) and is often easier to measure from experiments than the yield stress [START_REF] Marmottant | Discrete rearranging disordered patterns, part II: 2D plasticity, elasticity and flow of a foam[END_REF]. The four independent dimensionless numbers (ε Y , Bi, Co, n) completely characterise the problem (table 2), like (W e, Bi, Co, n).

From now on, the Saramito (2009) model will be refered to as "the EVP model", while the Herschel-Bulkley model will be refered to as "the VP model".

4 Overview of the solutions of the EVP model 4.1 Homogeneous flows: transient response to a simple shear [START_REF] Saramito | A new constitutive equation for elastoviscoplastic fluid flows[END_REF][START_REF] Saramito | A new elastoviscoplastic model based on the Herschel-Bulkley viscoplasticity[END_REF] studied simple flows in a geometry without stress heterogeneity, such as uniaxial extensional flow, oscillating shear flow or simple shear flow. To enable comparison with the circular Couette geometry (section 4.2), we first study the prediction of the EVP model in simple shear flow. The fluid is initially at rest: at t = 0, τ = 0. Then for t > 0, a constant shear rate γ is applied.

The solution τ (t) is then computed from eq. ( 10). As long as there is no stress heterogeneity, the EVP model does not predict any localisation. Figs. 4a,b plot the normalised shear stress growth coefficient η + S = τ / γ and the first normal stress difference N 1 = τ 11τ 22 with respect to the applied shear γt. At first, when the stress in the material is still below the yield stress τ Y , the shear stress increases linearly with time while the first normal difference increases quadratically: the material behaves as an elastic solid obeying the Poynting law [START_REF] Höhler | A constitutive equation describing the nonlinear elastic response of aqueous foams and concentrated emulsions[END_REF]. Such a non linear phenomenon has been seen experimentally in foams [START_REF] Labiausse | Shear induced normal stress differences in aqueous foams[END_REF]. After this initial elastic transient, saturation occurs at large applied shear: the stress components tend to a constant value as the applied shear tends to infinity. At the transition to the steady state, one can observe an overshoot of the shear stress that is more pronounced for small values of n. S versus W e. For 0 < n < 1, the shear viscosity decreases monotonically while it tends to a plateau when n = 1. Thus, when n < 1, the material is shear thinning. Observe that the value of Bi controls the viscosity plateau at small values of W e while it has less influence on the viscosity for large values of W e.

0 1 2 3 0 2 4 6 8 t γ (a) η + S (t, γ) K γn-1 n = 1 n = 0.5 n = 0.2 0 1 2 3 4 0 2 4 6 8 t γ (b) N + 1 (t, γ) 2K γn n = 1 n = 0.5 n = 0.2 10 0 10 1 10 -2 10 -1 10 0 10 1 10 2 W e = K γn-1 /G n = 1 (c) η S ( γ) K γn-1 Bi = 4 Bi = 2 Bi = 1 Bi = 0 10 0 10 1 10 -2 10 -1 10 0 10 1 10 2 W e n = 0.5 (d) η S ( γ) K γn-1 Bi = 4 Bi = 2 Bi = 1 Bi = 0
In summary, before yielding, the material behaves as a linear elastic solid while after yielding it is described by a a nonlinear viscoelastic model through the power-law index n.

Cylindrical Couette geometry: startup flow

We can now study how the stress heterogeneity due to the cylindrical Couette geometry modifies the simple shear flow presented in section 4.1.

Let us consider an EVP material described by eqs. (1),( 2),( 10), initially at rest, i.e. such that v = 0 and |τ d | < τ Y . For t > 0, the inner cylinder moves with a velocity V > 0 and the flow develops throughout the gap. The initial spatial distribution of stress reflects the preparation of the material [START_REF] Labiausse | Shear induced normal stress differences in aqueous foams[END_REF]. At first, while |τ d | < τ Y , no plasticity occurs and the equations can be solved explicitly for n = 1 (appendix C). The component τ rr is constant and equal to its initial condition, while the shear stress τ rθ (r, t) is linear in the shear strain:

τ rθ (r, t) = τ rθ (r, 0) - G 2 1 -Co Co 2 (2 -Co) ∆r r 2 V t ∆r , (12)
and the material develops normal stresses which are quadratic in the shear strain:

τ θθ (r, t) = τ θθ (r, 0) + G 2 1 -Co Co 2 (2 -Co) ∆r r 2 V t ∆r 2 . ( 13
)
Below the yield stress, the material exhibits the same elastic behaviour as in simple shear, except for the 1/r 2 spatial heterogeneity of the shear stress. Because of this heterogeneity, the norm of the deviatoric stress |τ d | first reaches the yield stress τ Y at the inner boundary. The plasticity dissipates the applied deformation and a steady state is rapidly reached, as shown by [START_REF] Cheddadi | Numerical modelling of foam Couette flows[END_REF]. When no external force such as friction is present, the saturation propagates instantaneously throughout the gap (see [START_REF] Cheddadi | Numerical modelling of foam Couette flows[END_REF]); therefore, the duration of the transient is of the same order of magnitude as the duration of the elastic regime.

Cylindrical Couette geometry: steady-state solutions

We now review some of the new insights gained from the tensorial EVP model, with an emphasis on the differences with scalar and/or VP models. Fig. 5 shows all the components of some solutions of the EVP model, for a set of parameters used as reference (eq. ( 14) and Fig. 2b), and compares them to the VP model. The exploration of these parameters will be presented in section 5. We define the critical shear rate γc as the jump of γ = 2D rθ at r = r c : it is visible as the slope of the velocity profile (see for instance Fig. 5a).

Range of possible initial conditions

The EVP constitutive equation ( 10) includes derivatives of the elastic stress tensor with respect to time and therefore allows the study of transient flows, as has been done by [START_REF] Cheddadi | Numerical modelling of foam Couette flows[END_REF]. This in turn requires to specify an initial condition that reflects the preparation of the material before the beginning of the experiment. In particular, the tensorial framework allows us to study various initial normal stresses.

Three different constant values for τ θθ are chosen. We denote by EVP 0 the case where the components of the initial stress are all set to zero. We then explore two limit cases, in which the component τ θθ is such that the norm of the initial deviatoric stress tensor is |τ d (r, θ, t = 0)| = τ Y , while the other stress components are zero τ rθ (r, θ, t = 0) = τ rr (r, θ, t = 0) = 0. We thus introduce the two cases EVP ± corresponding to initial stresses τ θθ (r, θ, t = 0) = ± √ 2τ Y . Note that the component τ rr (r, θ, t) remains constant and equal to its initial value zero in the cylindrical Couette geometry (see eq.( 10)). We denote by r + c (resp. r - c ) the critical radius corresponding to EVP + (resp. EVP -).

Results

The results from the EVP model are dramatically different to the results from the VP model, for which the (steady-state) velocity is unique, the normal stress components are zero and the localisation length r c is also unique for a given set of its parameters. Effect of the initial conditions. Surprisingly, Fig. 5 shows that the steady-state solution depends on the initial conditions. We find that they lead to three different steady-state solutions, both for the velocity profile (Fig. 5a) and for the stress (Fig. 5b-d). For a given set of parameters (ε Y , Bi, Co, n), when the initial condition varies the EVP model predicts a continuous set of steady-state solutions.

Non-uniqueness of the localisation length. Depending on the initial condition, the critical radius r c can reach any value in the range [r - c , r + c ]; the highest critical shear rate is reached for the solution with EVP -. We find that r + c = r 0 + 0.33∆r is close to the value predicted by the VP model, but r - c = r 0 + 0.24∆r is significantly smaller, while the critical radius r 0 c = r 0 + 0.28∆r for EVP 0 is in between.

Memory effects. Unlike the VP model, the EVP model exhibits non-zero normal stresses (Fig. 5b): in the flowing region (r < r c ) where plastic rearrangements occur continuously (|τ d | > τ Y ), the normal stress is independent of the initial conditions, and the material progressively loses memory of the initial condition. Conversely, in the non-flowing region (r > r c and |τ d | < τ Y ), the normal stress depends strongly on the initial conditions. The material thus keeps a record of the initial conditions through residual normal stresses.

Smooth and non-smooth solutions. Unlike the VP model, the EVP model can exhibit either smooth (EVP + ), or non-smooth profiles (EVP 0 , EVP -): the velocity gradient and the normal stress can be discontinuous at r = r c (Fig. 5a). Compared to EVP 0 , EVP -exhibits a stronger localisation (r - c < r 0 c ), a higher γc , and a higher jump of the normal stress (Fig. 5b).

Comment on the non-uniqueness of the steady-state solutions

First, we point out that, for each given initial condition, the corresponding time-dependent problem is well-posed, its time-dependent solution is unique, and the associated steady-state solution is unique too. Let us then analyse the sources of the non-uniqueness of the solution.

First, from a qualitative point of view, the non-uniqueness of the steady-state solution can be related to the transient elastic regime: in the extreme EVP -case, the component τ θθ starts with a negative value that quadratically evolves with time (eq. ( 13)); meanwhile, the shear stress develops throughout the gap (eq. ( 12)), and as its contribution to the von Mises criterion dominates that of the normal stresses (see eq. ( 4) and Figs. 5b,c), the yield stress value τ Y is reached while the normal stresses are still negative. In the flowing region near the inner cylinder, the plastic term is not zero anymore and strongly couples the component of the normal stresses to the flow. Their value is prescribed by the flow and no longer depends on the initial condition in the flowing region (Fig. 5b); it differs from the nonflowing region where the normal stresses profiles strongly depend on the initial condition.

When the steady-state regime is reached, these two regions do not join up, which results in discontinuous normal stresses and shear rate. The jump decreases as the initial stress is increased (solutions EVP 0 and EVP + , Fig. 5).

Second, from a formal point of view, the VP model (eq. ( 3)) contains one non-linearity, related to the von Mises criterion. The EVP model (eq. ( 10)) adds a second nonlinearity contained in the Oldroyd derivative (eq. ( 8)). Due to the expression of this Oldroyd derivative, the steady-state EVP equations do not reduce to the VP model. The additional term couples the normal stresses components with the shear stress and the velocity gradient (appendix D). Therefore, even though the velocity gradient reduces here to the shear component, the EVP steady-state solution develops non-zero normal stresses. This additional nonlinearity, interacting with the von Mises criterion, causes non-uniqueness of the stresses in the vicinity of |τ d | = τ Y , i.e. in the vicinity of r = r c .

Comparison with experiments

We explain here how the parameters of the EVP model were chosen in order to fit the experimental data shown in Fig. 2.

Smooth velocity profiles: the Debrégeas et al. experiment

In the foam experiment shown in Fig. 2a, image analysis was used to find the steady-state velocity profiles [START_REF] Debrégeas | Deformation and flow of a two-dimensional foam under continuous shear[END_REF]; it was reanalysed by [START_REF] Janiaud | Foam in a two-dimensional Couette shear: a local measurement of bubble deformation[END_REF] who measured the shear and normal components of the local elastic strain, in both the transient and steady-state regime.

In the experiment, a steady state is reached after a transient regime. The rotation direction is then inverted, and after a second transient another steady state is reached: this is when measurements are recorded [START_REF] Debrégeas | Deformation and flow of a two-dimensional foam under continuous shear[END_REF]. To perform a comparison, we use the experimental steady state as an initial condition, invert the rotation, compute the numerical solution until a first steady state is reached, then again invert the rotation, and compare the corresponding steady state with experiments.

The value of the curvature number is fixed by the geometry:

Co = 1 - r 0 r e = 1 - 71 122 ≈ 0.41803.
We have to adjust the value of the three remaining parameters: ε Y , Bi, n. We focus on the steady-state measurements of elastic strain and velocity.

We start with the yield strain ε Y ; this parameter is independent of the velocity and exerts a large effect on the elastic strain ε (e) = τ rθ /(2G) (see Cheddadi et al. (2011b)). We find that ε Y = 0.175 allows a good fit of the measured value of the shear elastic strain (Fig. 6a). The 0 0.1 [START_REF] Debrégeas | Deformation and flow of a two-dimensional foam under continuous shear[END_REF] analysed by [START_REF] Janiaud | Foam in a two-dimensional Couette shear: a local measurement of bubble deformation[END_REF]; errors bars have been provided by E. Janiaud (private communication). initial normal components of the strain are not affected (Fig. 6b) by the successive rotations in the region r > r c , and remain below the yield strain ε Y : this region undergoes reversible elastic deformation.

ε Y 0.2 0.3 0 r c 0.5 1 r-r 0 ∆r τ 2G (a) τ rθ /( √ 2G) exp. EVP d τ θθ /(2 √ 2G) exp. EVP d τ rr /(2 √ 2G) exp. EVP d 0 0.1 ε Y 0.2 0.3 0 r c 0.5 1 r-r 0 ∆r |τ d | 2G (b) exp. EVP d
Then we have to adjust the values of Bi and n. These parameters have much less effect than ε Y on the profiles of ε (e) , but they do affect the localisation length and the smoothness of the solution (see Cheddadi et al. (2011b)). As this experiment exhibits a very smooth transition from the flowing to the non-flowing region (exponential-like decrease of the velocity), our sensitivity analysis (see below, section 5.4) yields a good agreement with the data when using a small value of the power-law index, n = 1/3.

The localisation length can be measured from the experiment: r cr 0 ≈ 0.3 ∆r; we can use this information in the EVP model, taking advantage of the fact that the localisation is a VP effect and is mostly determined by the underlying VP model. This last model yields a relation between the values of r c , Bi, and n (see appendix A, equation ( 16)): with the chosen value of Co and n = 1/3, it yields Bi ≈ 7.6. This value is slightly adjusted in the EVP model in order to improve the fit: we obtain a better fit of the velocity profile with Bi = 10.

Non-smooth profiles: the Coussot et al. experiment

The experiment shown in Fig. 2b was made with a bentonite-water suspension, using MRI [START_REF] Coussot | Coexistence of liquid and solid phases in flowing soft-glassy materials[END_REF]. MRI measurements provide a precise and sharp velocity profile that could not be obtained with foams (section 6.1). We use it here to explore the solutions of the EVP model, without claiming to explain the particular properties of bentonite (which is thixotropic) or any other specific material.

For this experiment, only the steady-state velocity has been measured, which makes it more difficult to evaluate precisely the parameters of the EVP model. Since the normal stresses are not measured, we take an initial condition set to zero for the sake of simplicity. As the velocity profile is quite abrupt in the vicinity of r = r c , we choose n = 1 for the power-law index. Then, as for the previous experiment, we evaluate the Bingham number Bi under the condition that the critical radius predicted by the underlying VP model matches the value in the experiment (r cr 0 = 0.28 ∆r). From this point, the slope of velocity at r = r c (which is directly related to the critical shear rate) can be adjusted by tuning the yield strain ε Y . The best fit is obtained with ε Y = 0.35 and Bi = 27.

Sensitivity to the parameters

Now, we study how the range [r - c , r + c ] and the critical shear rate γc depend on the parameters of the EVP model. We explore the parameter space with the dimensionless numbers (ε Y , Bi, Co, n). It allows us to probe the effect of the imposed velocity V , through the Bingham number, as it is the only one that depends on V ; the geometry, through the curvature number Co; and the two material parameters: the yield strain ε Y , and the power-law index n. These parameters are varied around a reference set of parameters used for the comparison with [START_REF] Coussot | Coexistence of liquid and solid phases in flowing soft-glassy materials[END_REF] (section 4.4.2, Fig. 2b):

ε Y = 0.35, Bi = 27, n = 1, Co = 1/3. (14)

Effect of the inner cylinder velocity V

The effect of the imposed velocity V is probed through Bi, which is varied in the range [10,60] around the reference value Bi = 27. If V 0 is the velocity associated with this reference value, the range of Bi corresponds roughly to the range [V 0 /3, 3V 0 ] for the velocity. Observe on Fig. 7a that when the Bingham number Bi increases then r c decreases and the size of the zero-velocity zone increases. The critical radius r + c almost matches the critical radius predicted by the VP model. The flow becomes more localised and the gap between r - c and r + c also decreases slowly. The marker represents our reference solution (on Fig. 7a). We

Effect of the curvature Co

The curvature number Co explores the effect of the curvature of the geometry. Recall that Co → 0 corresponds to the plane Couette, with a homogeneous stress throughout the gap, while Co close to one corresponds to cylindrical Couette with a tiny central cylinder, and a highly heterogeneous stress. Observe on Fig. 8 that the size of the zero-velocity zone decreases when Co tends to one. The critical radius r + c of the smooth solution is well predicted by the VP model. Note also that the difference r + cr - c decreases with Co: thus, when r c is close to r 0 , the effect of the initial condition on the localisation length is less visible. Also the maximal discontinuity γc of the critical shear rate increases with Co.

Effect of the yield strain ε Y

The yield strain ε Y has only a small impact on the EVP + solution and a large effect on the velocity profile of the EVP -solution: r - c decreases and γc increases, both quasi-linearly, when ε Y increases. A linear regression through our numerical data leads to:

r - c -r 0 ∆r ≈ -0.23 ε Y + 0.33 - ∆r V γc ≈ 5.5 ε Y . (15)
Note that both the range r + cr - c and the discontinuity γc increase with ε Y . Accordingly, when ε Y vanishes, r - c = r - c and γc = 0: the solution is unique and smooth: there is no elasticity and the model reduces to the VP one.

Effect of the power-law index n

From Fig. 10, the range r + cr - c appears to remain roughly constant while the size of the yielded zone increases with n. The critical radius r + c decreases with n, as for the VP model. The γc discontinuity of the shear rate for the extreme EVP -solution decreases with n. These effects are not intuitive.

Another way to explore the effect of n is to vary it while r c is held fixed. In fact, the critical radius r c is approximately given by the VP model, and can then be predicted as a function of (Bi, n, Co) (appendix A). Inverting this function yields Bi(r c , n, Co): this value of Bi (table 3) is then injected in the EVP model. Using this procedure, the values of n and Bi are varied in such a way that r c remains constant. solution, associated with r + c ≈ r 0 + 0.3∆r while Fig. 11b plots the EVP -non-smooth one, with r - c ≈ r 0 + 0.23∆r. Note that both r + c and r - c are now roughly constant, and that the velocity profile is more curved when n is small. Fig. 11c shows γc : observe that the discontinuity decreases rapidly when n decreases. For n = 1/3, the non-smoothness becomes imperceptible, while the EVP + smooth and EVP - non-smooth profiles are very close and almost indiscernable, as shown on Fig. 11d. Note that even when n = 1/3, there are still multiple solutions, and the range for r c remains as large as [r - cr 0 , r + cr 0 ] ≈ [0.23, 0.3]∆r, despite the fact that these solutions share similar curved profiles.

Discussion

Smooth and non-smooth profiles

The debate about smooth and non-smooth velocity profiles, mentioned in the introduction, is difficult to address experimentally, for two reasons. First, since the steady-state solution is not unique, both smooth and non-smooth profiles can be observed in the same experiment. This depends on the residual stresses due to the initial preparation, which are usually not reproducible and are certainly difficult to suppress [START_REF] Labiausse | Shear induced normal stress differences in aqueous foams[END_REF][START_REF] Raufaste | Discrete rearranging disordered patterns: prediction of elastic and plastic behavior, and application to two-dimensional foams[END_REF], or on the presence of some impurities. This high sensitivity might explain some discrepancies in the literature, and the difficulty to settle the debates regarding localisation.

Second, the experiments in foams are not always precise enough to discriminate between smooth and non-smooth transitions. Fig. 12 compares experimental measurements and the two solutions: EVP + , the smooth one, with dashed lines, and EVP -, the non-smooth one, with solid lines. Experiments were performed with bubble rafts, and in order to reflect the absence of top and bottom plates, the fluid viscosity is introduced in computations as a second Newtonian viscosity η 2 with a viscosity ratio α = η/(η + η 2 ) [START_REF] Saramito | A new constitutive equation for elastoviscoplastic fluid flows[END_REF][START_REF] Saramito | A new elastoviscoplastic model based on the Herschel-Bulkley viscoplasticity[END_REF]. In that case, the smooth solution predicts r + c = r 0 + 0.65∆r and the non-smooth one r - c = r 0 + 0.43∆r together with a critical shear rate γc = 0.33 V /∆r. Observe that neither the smooth or the non-smooth solutions can be distinguished from the experimental measurements.

6.2 Is Couette flow suitable for caracterizing EVP materials?

Flow geometry

We recommend to study flows where residual stresses do not affect the understanding and measurements. Other requirements should include: well-defined boundary conditions; a good separation of scale between the discrete units, the representative volume element (RVE) and the global flow scale; the possibility to have a large variation of the control parameters such as velocity (and in foams, liquid fraction and bubble size dispersity); a variety of measurements providing stringent tests on the EVP model and its parameters. We have shown using experiments and models that a foam flow around an obstacle (Stokes flow) meets these requirements (Cheddadi et al., 2011a).

Couette flows, due to their simple geometry, have a long history of being used to probe the rheology of Newtonian fluids. They are also suitable for complex liquids such as VE or VP materials. However, the present study questions their use in EVP materials, and especially in foams, which are usually excellent model materials to perform in-lab experiments. In fact, in EVP materials, the initial conditions and the preparation method create residual stresses which are difficult to remove and affect the flow, which thus becomes non-unique and poorly controlled. Care is necessary to interpret the results. Future experimental work could try to deepen our understanding of localisation and to test our predictions by working in the following two directions.

Sample preparation

The control of the sample preparation, and of the initial normal stresses, is important.

Obtaining different (although uncontrolled) initial conditions should be reasonably easy by trying different methods to fill the experimental cell.

Choosing these initial stresses is less easy, but possible for instance by performing a high velocity pre-shear in the reverse direction. If the shear rate is high enough so that no localisation is observed in the steady-state regime, the normal stresses reach a value larger than √ 2τ Y . Whatever the initial loading of the material, it has been irreversibly erased by the plastic rearrangements. After such a pre-shear, one should observe smooth velocity profiles, with no critical shear rate. We have numerically checked that after such a high velocity pre-shear the EVP model predicts a smooth solution (similar to EVP + ), even with the EVP -initial loading, that would lead without pre-shear to a non-smooth solution.

Suppressing these initial stresses is approximately (but not completely) possible by a careful sample deposition in the experimental cell [START_REF] Labiausse | Shear induced normal stress differences in aqueous foams[END_REF]. Complete suppression would for instance require applying cycles of oscillating strain of amplitude which begins at around twice the yield strain and then gradually decreases to zero [START_REF] Raufaste | Discrete rearranging disordered patterns: prediction of elastic and plastic behavior, and application to two-dimensional foams[END_REF]; since these cycles have to be performed in each of the three dimensions, such suppression of initial stress would require a dedicated experiment and is in fact easier in simulations.

Finally, some coarsening materials such as foams become gradually isotropic with time, so that their trapped stresses slowly decrease [START_REF] Höhler | Rheology of liquid foam[END_REF].

We emphasize the fact that these two last procedures tend to suppress normal stresses and therefore correspond to the EVP 0 initial loading and non-smooth solutions.

Measurements of parameters and variables

We recommend to perform local measurements of strain or stress. In situ measurement of strain is purely geometric and independent of any knowledge of the sample physics: it can be performed on most systems where the positions of each discrete constituent object is measurable [START_REF] Graner | Discrete rearranging disordered patterns, part I: robust statistical tools in two or three dimensions[END_REF], which includes two-dimensional foams [START_REF] Janiaud | Foam in a two-dimensional Couette shear: a local measurement of bubble deformation[END_REF][START_REF] Marmottant | Discrete rearranging disordered patterns, part II: 2D plasticity, elasticity and flow of a foam[END_REF] or colloids. In situ measurements of stresses are possible with similar methods by measuring the positions of the discrete constituents and having some knowledge of their physical interactions, for instance in two-dimensional foams [START_REF] Marmottant | Discrete rearranging disordered patterns, part II: 2D plasticity, elasticity and flow of a foam[END_REF], or in birefringent materials by using photoelasticity [START_REF] Miller | Stress fluctuations for continuously sheared granular materials[END_REF].

Model parameters are all measurable in principle. Note that the yield strain or stress requires a tensorial measurement, and thus normal differences components [START_REF] Labiausse | Shear induced normal stress differences in aqueous foams[END_REF]. Again, local measurements of strain can help in measuring the yield strain directly [START_REF] Marmottant | Discrete rearranging disordered patterns, part II: 2D plasticity, elasticity and flow of a foam[END_REF].

Summary

We solve here a tensorial model for the cylindrical Couette flow of elastic, viscous, plastic materials. We provide an approximate expression for the rheology versus different material parameters. In turn, our predictions are compared with experiments. We show that there is a complex interplay between elasticity, viscosity and plasticity, which together (but not separately) account for experimental observations. Even in such a simple geometry, the orientational effects are important, so that a tensorial EVP description is necessary to capture many aspects of the physics: reversible elastic deformation (both shear and normal components) below τ Y , following the Poynting law; memory of the preparation of the material through the initial stress condition, and consequently non-uniqueness of the steady flow and persistent residual normal stresses; possible appearance of non-smooth solutions. The normal stresses can be interpreted as an intrinsic structure parameter at the macroscopic level, and therefore independent of the underlying microstructure. These features can be predicted neither by VP models nor by VE models.

We have computed numerically the value of the localisation length versus different parameters so that we can guide experimentalists to design experiments that may or may not exhibit such effects: the effects of the initial conditions are more visible for instance when ε Y is large; or when the dissipation exponent n is large; or when the velocity is large, but small enough to allow for a localisation within the gap; or when the heterogeneity from the geometry is small, but large enough to allow for a localisation within the gap. Altogether, our results provide a validation of the continuous material description, a determination of EVP material parameters, and an in-depth understanding of their complex rheology. Finally, it appears that the steady Couette flow, which has stimulated so many debates, is neither robust nor unique. Despite its apparent simplicity, it involves numerous parameters, such as the initial conditions for the stress tensor, and is difficult to use in practice. Couette flow experiments could be complemented with flows in other geometries, with a stronger dependence in time and/or space.

A Herschel-Bulkley solution in cylindrical geometry

When ε Y = 0 the EVP model reduces to the VP one, and the velocity profile is given by

v θ (r) = √ 2 r 2 c τ Y K 1 n r rc r 1 s 1 s 2 - 1 r 2 c 1 n ds
In the case where the flow is driven by the inner boundary, the critical radius r c is given by the following expression:

τ Y = K √ 2 r 2 c       v(r 0 )
For fixed values of n and r c , the corresponding value of Bi for the VP model, denoted as Bi vp (n, r c ), can be easily computed by using numerical integration. Conversely, for given Bi and n, the value of r c associated to the VP model, denoted as r vp c (n, Bi), can be obtained from a small numerical computation (Fig. 13a). In this geometry, the localisation is observed as long as r c /∆r < r e /∆r = 1/Co. The VP model predicts that r c decreases with the Bingham number: the zero-velocity zone develops and the localisation effect is more pronounced. Observe also that the localisation effect is more pronounced when n decreases at fixed values of Bi and Co. The r c value associated to the VP model is denoted as r c,vp (n, Bi) and its inverse function as Bi vp (n, r c ). At fixed n, r c depends on Bi roughly as a power-law:

B Numerical method

The velocity is approximated by continuous affine finite elements while the stress components are piecewise constant over the mesh. The code is implemented by using the C + + Rheolef finite element library [START_REF] Saramito | Efficient C++ finite element computing with Rheolef[END_REF]. The stopping criteria for a steady solution is satisfied when the residual term is less than 10 -8 . Fig. 14 shows the convergence versus the mesh size h at the vicinity of r = r c for the EVP -non-smooth steady solution presented on Fig. 5.

Observe that the numerical method presents excellent convergence properties, despite the non-smoothness of the solution: the velocity is non-differentiable (Fig. 14a) while normal stress is discontinuous (Fig. 14b).

C Calculations for startup flow

The initial condition is such that v = 0 and |τ d | < τ Y throughout the gap. For t > 0, the inner cylinder moves with a velocity V > 0. The calculations that follow are valid as long as no plasticity occurs.

Eq. ( 23) yields τ rθ = -C(t)/r 2 , where C(t) depends only on time. With n = 1, eq. ( 19) yields τ rr = 0 if we choose τ rr (r, 0) = 0 in the initial condition. As a results, eq. ( 20) implies

D rθ = - C (t) 2Gr 2 .
Using the fact that D rθ = 1/2×r ∂(v/r) ∂r , and the boundary conditions v(r 0 , t) = V ; v(r e , t) = 0, we find The expressions ( 12) and ( 13) of τ rθ (r, t) and τ θθ (r, t) follow from eqs. ( 20) and ( 21).

C (t) = G 2 r 0

D Equations in cylindrical geometry

The velocity field is v = (0, v θ , 0) in cylindrical (r, θ, z) coordinates.

τ θθ = 0, , (21)
with |τ d | = 2τ 2 rθ + 1 2 (τ rrτ θθ ) 2 1/2 . Here, D = (∇v + ∇v T )/2 is the rate of strain tensor, K denotes a generalised viscosity [START_REF] Saramito | A new elastoviscoplastic model based on the Herschel-Bulkley viscoplasticity[END_REF]. The conservation of momentum writes: This system of equations is closed by boundary conditions for the velocity at the inner and external cylinders, respectively r = r 0 and r = r e (Fig. 1a), and by initial conditions for both the velocity v θ and the elastic stress τ .

Figure 1 :

 1 Figure 1: Experimental set-up for a two-dimensional circular shear flow of a foam confined between two horizontal plates. (a) definition of the geometric and kinematic parameters ; (b) picture of the confined two-dimensional liquid foam (from Debrégeas et al. (2001)): the internal radius is r 0 = 71 mm.

Figure 2 :

 2 Figure 2: Smooth and non-smooth velocity profiles. Experimental data are compared with our computations as discussed in section 4.4. (a) Smooth profile: comparison of experimental data on a foam by Debrégeas et al. (2001) with the present computations (ε Y = 0.175, Bi = 10, Co = 0.41803, n = 1/3); (b) Non-smooth profile: comparison of experimental data on a bentonite-water suspension from Fig. 1.a of Coussot et al. (2002), with the present computations (ε Y = 0.35, Bi = 27, Co = 1/3, n = 1). Inset: zoom around r = r c .

Figure 3 :

 3 Figure 3: Schematic representation of the Saramito (2009) EVP model.

Figure 4 :

 4 Figure 4: Start-up shear flow for W e = η γ/G = K γn-1 /G = 1 and Bi = τ Y /(K γn ) = 1: (a) normalised shear stress growth coefficient η + S (t, γ) vs normalised time; (b) normalised first normal stress difference N + 1 (t, γ) vs normalised time. Steady shear flow: shear stress coefficient η S ( γ) vs W e for various values of Bi: (c) n = 1; (d) n = 0.5.

Fig

  Fig. 4c,d plot the steady shear viscosity η S = lim t→∞ η +S versus W e. For 0 < n < 1, the shear viscosity decreases monotonically while it tends to a plateau when n = 1. Thus, when n < 1, the material is shear thinning. Observe that the value of Bi controls the viscosity plateau at small values of W e while it has less influence on the viscosity for large values of W e.

Figure 5 :

 5 Figure 5: Comparison between the steady-state solutions of the VP Herschel-Bulkley model and the EVP Saramito model in cylindrical Couette geometry: (a) velocity, with same experimental data as in Fig. 2b (Coussot et al., 2002) (b) normal stress; (c) shear stress; (d) norm of the deviatoric stress. Lines: solutions of the EVP model for the reference set of parameters (eq. (14)): ε Y = 0.35, Bi = 27, Co = 1/3, n = 1, as on Fig. 2b, with different initial conditions. Thick solid red line (denoted EVP 0 ): τ 0 θθ = 0; thin dashed black line (EVP + ): τ 0 θθ = √ 2τ Y ; thin dotted blue line (EVP -): τ 0 θθ = -√ 2τ Y ; thick solid green line (VP): Co = 1/3, Bi = 27, n = 1. We denote r + c (resp. r - c ) the critical radius of the EVP + (resp. EVP -) solution.

Figure 6 :

 6 Figure 6: The stress tensor: comparison of computations and experiments: (a) components; (b) norm of the tensor. Lines: computations with ε Y = 0.175, Bi = 10, Co = 0.41803, n = 1/3. Symbols: experiments (Debrégeas et al., 2001) analysed by Janiaud and Graner (2005); errors bars have been provided by E. Janiaud (private communication).

Figure 7 :

 7 Figure 7: Effect of the Bingham number Bi. (a) r - c and r + c . Red bars: five examples of the range in which (r c -r 0 )/δr varies. Black cross: parameters as (14) together with τ (t = 0) = 0. Solid black line: value predicted by the VP model with the same parameters except that ε Y = 0. From Coussot et al. (2002) we get V 0 /∆r = 3.25 s -1 . (b) Corresponding values of γc .

Figure 8 :

 8 Figure 8: Effect of the curvature Co: same caption as Fig. 7.

Figure 9 :

 9 Figure 9: Effect of the yield strain ε Y : same caption as Fig. 7. Dashed lines: fits of the EVP -solutions to the data. n 1/3 1/2 2/3 1 Bi 9.54785 12.7851 16.6489 27

Figure 10 :

 10 Figure 10: Effect of the power-law index n: same caption as Fig. 7.

Figure 11 :Figure 12 :

 1112 Figure 11: Effect of the power-law index n at fixed r : (a) EVP + smooth and (b) EVP -nonsmooth velocity profiles, (c) γc ; (d) smooth (dashed and black) and non-smooth (continuous and red) profiles for n = 1/3.

Fig

  Fig.13aplots r c as a function of the Bingham number Bi, for ε Y = 0 (VP model) and Co = 1/3, and different values of the power-law index n. Note that the value Co = 1/3 corresponds to the geometry of the experiment by[START_REF] Coussot | Coexistence of liquid and solid phases in flowing soft-glassy materials[END_REF] as shown on Fig.2b. In this geometry, the localisation is observed as long as r c /∆r < r e /∆r = 1/Co. The VP model predicts that r c decreases with the Bingham number: the zero-velocity zone develops and the localisation effect is more pronounced. Observe also that the localisation effect is more pronounced when n decreases at fixed values of Bi and Co. The r c value associated to the VP model is denoted as r c,vp (n, Bi) and its inverse function as Bi vp (n, r c ). At fixed n, r c depends on Bi roughly as a power-law:

Figure 13 :Figure 14 :

 1314 Fig. 13b represents β n vs n. A nonlinear regression leads to the dotted line of equation:β n ≈ -0.38 n 2 + 0.88 n -1.02,(18)

Table 2 :

 2 Four dimensionless numbers (ε

	Herschel and Bulkley (1926) proposed a power-law variant of the viscoplastic Bingham (1922)

Y , Bi, Co, n), or equivalently (W e, Bi, Co, n), which completely characterise the

[START_REF] Saramito | A new elastoviscoplastic model based on the Herschel-Bulkley viscoplasticity[END_REF]

model in the circular Couette geometry.

  The elastoviscoplastic constitutive equation in cylindrical coordinates writes:

	1 2G 1 2G 1 2G	∂τ rr ∂t ∂τ rθ ∂t ∂τ θθ ∂t	+ max 0, -2D rθ τ rr + max 0, -4D rθ τ rθ + max 0,	|τ d | -τ Y 2K|τ d | n |τ d | -τ Y 2K|τ d | n 2K|τ d | n |τ d | -τ Y	1 n 1 n 1 n	τ rr = 0, τ rθ = D rθ ,	(19) (20)

Couette flows for various EVP materials such as emulsions or liquid foams, including the non-reproducibility that has been reported in flows of foams. We then discuss the suitability of Couette flows as a way to measure rheological properties of EVP materials.

Acknowledgements

We warmly thank S. Cox and G. Ovarlez for detailed comments on the manuscript, M. Dennin and E. Janiaud for providing data, S. Cohen-Addad, B. Dollet and C. Gay for fruitful discussions, the French Groupe de Recherches "Mousses et Emulsions", C. Gay and C. Philippe-Barache for funding the travels for collaborations. FG thanks D. Weaire for intense discussions and questions which stimulated this work.