Understanding and predicting viscous, elastic, plastic flows.

Foams, gels, emulsions, polymer solutions, pastes and even cell assemblies display both liquid and solid mechanical properties. On a local scale, such "soft glassy" systems are disordered assemblies of deformable rearranging units, the complexity of which gives rise to their striking flow behaviour. On a global scale, experiments show that their mechanical behaviour depends on the orientation of their elastic deformation with respect to the flow direction, thus requiring a description by tensorial equations for continuous materials. However, due to their strong non-linearities, the numerous candidate models have not yet been solved in a general multi-dimensional geometry to provide stringent tests of their validity. We compute the first solutions of a continuous model for a discriminant benchmark, namely the flow around an obstacle. We compare it with experiments of a foam flow and find an excellent agreement with the spatial distribution of all important features: we accurately predict the experimental fields of velocity, elastic deformation, and plastic deformation rate in terms of magnitude, direction, and anisotropy. We analyse the role of each parameter, and demonstrate that the yield strain is the main dimensionless parameter required to characterize the materials. We evidence the dominant effect of elasticity, which explains why the stress does not depend simply on the shear rate. Our results demonstrate that the behaviour of soft glassy materials cannot be reduced to an intermediate between that of a solid and that of a liquid: the viscous, the elastic and the plastic contributions to the flow, as well as their couplings, must be treated simultaneously. Our approach opens the way to the realistic multi-dimensional prediction of complex flows encountered in geophysical, industrial and biological applications, and to the understanding of the link between structure and rheology of soft glassy systems.

Materials such as pastes or polymer solutions display both solid-like and liquid-like behaviors; they are successfully described by visco-elastic (VE) or visco-plastic (VP) models. However, we still lack testable predictions of the timeand space-dependent flows of soft glassy materials [6,7], that are made of disordered assemblies of deformable, rearranging units [6,19]. It had been suggested that fluctuations remain relevant even at large scale, in which case detailed statistical theories of long-range correlations and avalanches would be required [20,21]. This view is challenged by recent experiments which suggest that even these materials can be treated as continuous materials described by tensorial equations [6,14]; thus in principle partial differential equations could lead to the long-awaited predictions.
Based on our experience with foams, we believe that the reason of the difficulty comes from the fact that these materials are simultaneously viscous, elastic, and plastic (VEP). Under small deformation, a foam reversibly comes Send offprint requests to: pierre.saramito@imag.fr, francois.graner@curie.fr a Present address: Laboratoire de Physique de la Matière Condensée, UMR 6622 CNRS and Univ. Nice-Sophia Antipolis, Parc Valrose, F-06108 Nice Cedex 2, France back to its shape; at large deformation, it can be irreversibly sculpted and gets a new shape; under an increasing deformation rate, it irreversibly flows, with an increasing viscous stress [22][23][24][25][26]. Existing continuous models of foam flows include for instance either a phenomenological scalar description [17,27], or a complete tensorial description of the elasticity [9] or plasticity [19] based on the micro-structure. Overall, continuous VEP models tend to successfully reproduce some experimental measurements, such as elastic and loss moduli, or compliance (see for instance refs. [3,9,11,25]).
We want here to understand and predict VEP flows.
Our approach is to test whether a continuous VEP model To really test a model, and also lead to practical applications, we address the full spatial and orientational heterogeneity of a flow. We need to investigate a controlled,  (8)(9)). This latter point, overlooked by most models, is crucial to test whether the elastic stress could be entirely determined by the shear rate, leading to a VP rheology (such as Bingham [28] or Herschel-Bulkley [29]). Recent experiments [15,16] suggest that it is not the case: in fact, the stress, the shear rate and the elastic deformation should be treated as independent variables, so that a full VEP treatment is required.
One such flow is the well documented flow around an obstacle [12] (Fig. 2). It displays a strong spatial heterogeneity, simultaneous VEP behaviors, a large range of elastic deformations, several elongation and rotation rates, and various relative orientations of the relevant tensors [19]. It enables to follow a bubble at different stages, while it stretches, then while it relaxes: thus, even in a steady flow, transient effects and relaxation times are apparent.
It is classically used as a stringent test to discriminate between different models [30] (Fig. 3).
We use foams as model systems of VEP materials. Experiments with foams or emulsions, especially in 2 dimensions, enable an easy, simultaneous visualisation of the micro-structure (bubbles or droplets, which act as tracers of velocity and deformation) and the large scale (global flow heterogeneities).
To compare mesurements from discrete experiments with continuous predictions from partial differential equations, we use the experimental tools we have developed [19]. V, E and P contributions are expressed in the same units, favoring a unified description of solid and liquid behaviors; each of them is valid in all regimes (so that e.g. elastic contributions can be measured even out of the elastic regime). These are local (in situ) measurements which link the foam structure and rheology.
Our plan is as follows. In order to make this paper It has already been used for space-and time-dependent predictions [11]. It has first been used to solve simple cases, such as steady uniaxial elongation. It has also applied to oscillatory regimes, and calculations of G ′ , G", rigidity and loss moduli.
It has then been implemented to calculate both timedependent and steady shear Couette flows, which depend on one space variable (circular [32] or planar [33] geometries). It involves two strong non-linearities, intrinsic to VEP flows, and thus independent of the model: one because the plasticity appears above a yield point (eq. 3), and the other because of the advection of elastic stress (eqs. (8)(9)). Despite these unavoidable difficulties, the model has been solved [34,35]. The resulting velocity, elasticity and plasticity fields agreed with experimental measurements.
This improved our understanding of Couette flows [35].
The localization of the velocity field results from the stress heterogeneity, so that the circular geometry by itself can induce localisation. In planar geometry, where the stress is a priori homogeneous, localisation necessarily relies on another cause of heterogeneity, such as an external friction.
Initial normal stresses can be preserved even in a steady flow, so that there are residual normal stresses which depend on initial conditions linked with the foam preparation method: the steady flow is not unique. Despite its simplicity, a Couette flow displays specific VEP features [34]. For instance, at the boundary of the localised region, the discontinuity of velocity gradient depends on the residual normal stresses, and thus on initial conditions. For all these reasons, and because the range of experimental data is limited, Couette flows have only a limited ability to discriminate between models, or between parameter values.

Constitutive equations: solid mechanics
We start with constitutive equations specific of a semifluid semi-solid material (see Fig. 1). In order to emphasize the dominant role of elasticity, we express them here in terms of deformations and their rates, as is usual in the context of solid mechanics: Eq. (1) is a constitutive equation for the Cauchy stress σ tot , with pressure, viscous and elastic terms. Here p is the hal-00565805, version 1 -14 Feb 2011 pressure, I the identity tensor;ε is the total deformation rate tensor, η 1 is the viscosity of the material apparent at small deformation (in the foam, it includes the dissipation inside the soap films); ε e is the elastic deformation tensor, µ the shear modulus.
Eq. (3) is a plasticity equation, which specifiesε p by stating that plasticity increases above the yield strain ε Y : under a strong shear rate, the elastic deformation |ε e | can become higher than ε Y . Here |ε e | is the norm of ε e , which we define as |ε e | = [(ε e yx ) 2 + (ε e xx − ε e yy ) 2 /4] 1/2 to facilitate the comparison with a scalar shear (fixed shear direction xy with a given amplitude ε), i.e. ε e xx − ε e yy = 0 and ε e yx = ε, so that |ε e | = ε [6]; the externally measured scalar shear is then γ = 2|ε e |. Another acceptable definition would be the euclidian norm of the deviatoric elastic strain tensor [35], In soft disordered materials, plasticity is related with local rearrangements. In foams, these happen when bubbles swap neighbors and are called "T1" processes [22,23,36]. They create a transient local deformation. Here λ is the relaxation time of the material after such a local deformation [37]. We can construct the dissipation due to plasticity, which has the dimension of an effective viscosity η 2 = λµ: it determines the loss modulus at large amplitude. The softness and deformability of the mate-rial appears in the value of ε Y , of order of unity, so that both the elastic and plastic behaviors are experimentally observable; its glassy (i.e. disordered) nature implies that λ, µ, ε Y are isotropic [6].

Conservation equations: fluid mechanics
Generic conservation equations for an isothermal flow are expressed in terms of the velocity v and its derivatives, as is usual in the context of fluid mechanics: div v = 0.
Eq. (4)  where k is a constant. We find experimentally (see Section 3.3) that k is small enough that in the present flows the effect of f ext is not measurable, so that we neglect it.
Eq. (5) describes the incompressibility of the flow: it applies to slow flows when the compressibility modulus is much higher than the shear modulus, as is the case in foams [22][23][24].

Closing equations: linking solid and fluid mechanics
To close the system of eqs. (1-3) and eqs. (4)(5) requires to couple the deformations to the velocity.

hal-00565805, version 1 -14 Feb 2011
First, the total deformation rate equals the symmetrized velocity gradient:ε Second, the time variation of the elastic deformation tensor ε e accounts for its advection by the flow velocity v.
The model should be objective, that is, the expression of the equations should remain the same for an observer who has a movement of translation or rotation with respect to the experiment. The advection of the elastic deformation tensor is thus described with a frame invariant tensorial derivative [38]:ε The objective derivative is [38,11]: where Here W (v) = 0.5(∇v − ∇v T ) is the antisymmetric part of the velocity gradient, and a ∈ [−1, 1] is the so-called "a parameter" [11,38], which effect is discussed in section 3.3.

Experimental methods
Experimental set-ups have been described in refs. [14,19,39]. Bubble distributions are monodisperse in size (area 16 mm 2 ) and disordered in geometry (shape, number of neighbors). While bubbles pass through the field of view, no rupture is observed, and coarsening is negligible.
The wet foam is prepared by direct bubbling into the 1 m long channel. At its entrance in the channel, it displays normal differences in elastic deformation, ε e xx − ε e yy (axis 2 in Fig. 7).
For the dry foam, bubbles pass first through a chamber (in which the foam drains): this chamber enables to vary the liquid fraction over more than three decades, and homogeneizes the foam while relaxing its normal differences in elastic deformation [39].
Experimental measurements treat solid and liquid behaviors with a unique set of mutually compatible tools [19]. We derive the continuous description directly from averages over almost a thousand of images of discrete measurements performed on all bubbles which can be automatically identified using image analysis, that is, which do not touch the obstacle. The entrance velocity V is measured ±2% as the average over the side of the field of view.
The texture (bubble size, elongation and packing) and its variation (bubble stretching and rearrangements) enable to measure in situ the velocity gradient (not shown), the elastic deformation ±2% and the plastic deformation rate ±7% [19]. We plot here deviatoric terms, see section 3.3.

Resolution
Eqs. (1)-(5) can be solved in 2D or in 3D. Their main difficulties are intrinsinc to VEP flows, independently of the details of the model. They reduce to a set of three partial differential equations with three unknowns (ε e , v, p) and hal-00565805, version 1 -14 Feb 2011 the coupled system is highly non-linear : its numerical resolution needs to be performed carefully.
Here we solve these equations in 2D with a finite element algorithm first used for a simple Couette circular geometry [35], extended here to handle more complex flow domains [34]. As in 1D, the stationary solution is obtained by solving the time-dependent problem with a second order time-splitting algorithm, already used for VE [40], generalized here to VEP; it allows to treat separatly the two main non-linearities of the equations, namely the plasticity term in eq. (3) and the stress transport term in the objective derivative (eq. (8)). Unlike in 1D [41], the stress transport term needs to be treated specifically by upwinding techniques; we chose a robust method, the discontinuous Galerkin scheme [42]. In addition, the nonlinearity linked to plasticity needs a much more careful discretization than for the Couette resolution [35] to ensure a proper decreasing of the residue of the stationnary problem. The spatial discretization is performed with a mixed finite-element method as in [30]. In order to get a general method suitable for any geometry, the domain is discretized with triangles.
The calculation domain is a channel, 15 R = 22,5 cm upstream and 30 R = 45 cm downstream of the obstacle.
The mesh, made of 1100 triangles, is locally refined near the obstacle (see Fig. 5a). We start from a foam at rest and enforce the entrance velocity V . Unlike for most liquid flows, but in agreement with foam flow experiments [14], we use slip boundary conditions. Careful tests have been performed [34] in order to ensure that the mesh is suffi-ciently refined and that discretization does not affect the results presented here. Iterations are performed (Fig. 5b) until the residue of the stationary problem is less than 10 −7 (see Fig. 5c). Calculations in 2D run in half a day on a Intel T7300 Core 2 Duo processor (2 GHz, 4 Mb cache, 32 bits). The 2D algorithm has been validated by reproducing the 1D Couette calculation [35], which runs in a few minutes.

Choice of parameters
The parameter with the main effect is ε Y : a change of ible. It means that here the limit k → 0 is regular. This case is similar to cylindrical Couette geometry, but different from planar Couette geometry [35].

hal-00565805, version 1 -14 Feb 2011
We choose η 1 /η 2 = 0.1 according to the Couette case [35], since it lies in the middle of a range where its exact value barely affects the flow, even up to a factor of 10.
We choose the co-rotational derivative [38], with a = 0, so that ε e is deviatoric. In that case, the term β a in eqs. The relaxation time λ is related to the Weissenberg number We = λV /R. We distinguish three velocity regimes.
At high velocity, We of order 1 or higher, the material rheology can display non-linearities in addition to those already present in the model. Since foams rupture at high velocities [39,43], this regime would be easier to investigate with other materials. In the low velocity range, We greater than 10 −2 but smaller than 1, the exact value of We does not affect the flow. This is the case for both experiments considered here, as well as for several foam flow experiments reported in the literature. Further decreasing We over two or three decades would lead to the ultra slow range, where the fore-aft asymmetry strongly increases.
This is done in very few well-controlled experiments [39] The limit We → 0 at constant ε Y is singular: it implies a divergence of the Bingham number Bi = 2ε Y /We [11]. In fact, the ultra slow regime is not quasi-static [25,27,44]. and does not exactly match quasi-static simulations [27].

Representations
Results are plotted either as maps or graphs.

Wet foam flow
First, as a preliminary characterization, we study the flow of a wet foam (7% liquid fraction, Fig. 2a). We calculate the measurable fields: velocity v, elastic deformation ε e , plastic deformation rateε p . We use the parameters obtained in ref. [35] and rescale them to the geometry of the present set-up. We investigate separately the effect of each parameter (see Section 3.3). We check that, with a yield strain ε Y = 0.1 ± 0.02, the calculations agree well simultaneously with all available experimental data ( Fig. 3 and Figs. [6][7][8]. Such value of ε Y is the expected order of magnitude for a foam with this liquid fraction [19]. Other parameters have less effect. We use λ = 0.2 ± 0.1 s [14] and η 1 /η 2 = 0.1. We take k = 0 without signif- breaking the up-downstream (fore-aft) symmetry (Fig. 3).
A velocity overshoot, the so-called "negative wake", is clearly visible behind the obstacle (Fig. 3). This characteristic feature of VEP flows is barely affected even if we vary V across the low velocity regime, confirming experimental observations [14]. For instance, dividing V by 20 barely changes the overshoot (Fig. 3).
This strongly contrasts with VP flows, which are always fore-aft symmetric [30]. VE flows represent a mixed situation, where the negative wake has already been evidenced, both experimentally [45] and numerically [46,47].
In fact, it occurs for low extensional viscosity fluids and models (e.g. FENE-CR [47] but not Oldroyd-B [31]), at elongational rates large enough in comparison with the inverse relaxation time, so that the elastic deformation does not vanish. However, at the low velocity investigated here, the VE flow is completely fore-aft symmetric (Fig. 3), and even indiscernable from viscous flows, whatever the viscosity. To interpret this set of observations, it seems that the overshoot appears when the elastic deformation ceases to follow passively the total deformation rate. This can occur if there is a mechanism which saturates the value of elastic deformation, which is the case in some VE models at high velocity, and in any VEP model because of plasticity.
To summarize this first test, by adjusting only one parameter, which value has the expected order of magnitude, we can adjust both qualitative and quantitative features hal-00565805, version 1 -14 Feb 2011 of all available data. We reproduce the observed negative wake and evidence the specificity of VEP flows.

Dry foam flow
Second, we turn to prediction. Since the parameter which has the most significant effect on the flow is the yield strain [14] (section 3.3), we choose to predict the flow for a twice larger value, ε Y = 0.2.
These predictions are plotted on Fig. 4a, top and Fig. 4be, lines. The overshoot on the symmetry axis (Fig. 4b) is larger and closer to the obstacle than for ε Y = 0.1: this reflects a larger effect of the elastic deformation. The elastic deformation field extends more than the velocity gradient ( Fig. 4a, top), which itselfs extends more than the plastic deformation rate field (Fig. 11, top): this confirms that the three fields are physically independent quantities [19].
To check this prediction, we then perform a dry foam experiment, since decreasing the liquid fraction of a foam is expected to affect especially the yield strain. With 1.2% liquid fraction, Fig. 2b), we observe that ε e max , the maximum value of ε e measured on the experiment, is twice as much that of the wet foam. The effect of elasticity is even stronger and the agreement with our prediction even better, without adjusting any parameter. Measurements confirm the predicted spatial distribution, magnitude, direction, anisotropy of fields.
The bubble velocity (Fig. 4) passes the three most stringent tests. First, the position and magnitude of overshoot on the symmetry axis (Fig. 4b). Second, the graph along the axis 5 (Fig. 4c). And third, the exact position of the arrest points, defined in the referential of the foam, as points where v x − V = 0: close to axis 5, on y = 0 and ±5 cm, see Fig. 4a. Other axes confirm the agreement (Fig. 9).
The bubble elastic deformation too agrees remarkably well in spatial distribution, amplitude and direction (Figs. 4a,d,e and Fig. 10). Interestingly, its orientation does not directly correlate with that of streamlines, or equivalently of the velocity vectors. This confirms that the elastic deformation should be treated as a variable independent of the velocity; we have also checked (data not shown) that it does not directly correlate with the total deformation rate. Note that residual normal stresses are visible beyond the obstacle. Unlike in the Couette case, here they are reproducible and their origin is understood: they are a direct effect of the obstacle, and do not depend on the foam preparation method.
The plastic deformation rate is calculated as the total deformation rate minus the elastic deformation rate (eq. 2). Its predicted spatial distribution and directions agree with that of the experimental measurements, which represent the time-averaged orientation, frequency and anisotropy of the bubble rearrangements (Fig. 11).

Discussion and conclusion
To summarize, a continuous description of viscous, elastic, plastic material with physically meaningful parameters can reproduce and even successfully predict a tensorial, spatially developed flow of a disordered rearranging structure.

hal-00565805, version 1 -14 Feb 2011
We analyse and interpret the effect of each parameter separately. We emphasize the dominant role of elasticity and thus identify the yield strain as the most important parameter. The flow does not reduce to VE or VP separately, so that we emphasize the specific complexity of VEP materials.
Our method opens the way to computing two-or threedimensional flows under any type or amplitude of deformation. It applies to those depending on one space variable: for instance a flow through a channel [15,16], or during simultaneous squeezing and shearing [6]. It also applies to those depending on two space variables: for instance a flow through a hole in 2D [19,48], or in 3D with axisymmetry.
At the expense of longer calculations, it can even apply to flows which depend on three space variables, for instance through a twisted or branching pipe.  [29] instead of a Bingham one, as observed in [17,50]. The friction on walls too scales non-linearly with high velocity [26]. Non-linear elasticity at large deformation, although seldom reached in foams [19], can too be taken into account [9,24]. At low velocity, plasticity seems to appear progressively: some bubbles begin to rearrange below the yield strain [19].
More generally, the model can be adapted specifically to any given VEP material of known properties. The value of the parameters of eqs. (1)-(5) depends on the microstructure, its disorder and its physico-chemical properties. This is where the present approach can be enriched by statistical models based on the micro-structure [7,15,16,20,21,24,44,51].     (bottom) Plots along axes 1-5. Same caption as Fig. 8.