Sylvain B Énito

Charles-Henri Bruneau

Thierry Colin

Cyprien Gay

François Molino

An elasto-visco-plastic model for immortal foams or emulsions

Keywords: numbers: 83, 10, GrConstitutive relations 83, 80, IzEmulsions and foams in Rheology 83, 50, AxSteady shear flows, viscometric flow 83, 85, LqNormal stress difference measurements

A variety of complex fluids consist in soft, round objects (foams, emulsions, assemblies of copolymer micelles or of multilamellar vesicles -also known as onions). Their dense packing induces a slight deviation from their prefered circular or spherical shape. As a frustrated assembly of interacting bodies, such a material evolves from one conformation to another through a succession of discrete, topological events driven by finite external forces. As a result, the material exhibits a finite yield threshold. The individual objects usually evolve spontaneously (colloidal diffusion, object coalescence, molecular diffusion), and the material properties under low or vanishing stress may alter with time, a phenomenon known as aging. We neglect such effects to address the simpler behaviour of (uncommon) immortal fluids: we construct a minimal, fully tensorial, rheological model, equivalent to the (scalar) Bingham model. Importantly, the model consistently describes the ability of such soft materials to deform substantially in the elastic regime (be it compressible or not) before they undergo (incompressible) plastic creep -or viscous flow under even higher stresses.

I. INTRODUCTION: ELASTICITY AND PLASTICITY IN FOAMS AND EMUL-SIONS A. From crystals to foams and emulsions

Historically, descriptions of deformations in crystalline, solid materials are based on a decomposition in terms of elastic and plastic components. Conceptual and technical problems arise in this process. On the one hand, general elastic formulations use continous deformations [START_REF] Landau | Theory of Elasticity[END_REF] with respect to a reference state associated with the ordered structure of minimal energy. On the other hand, plasticity is related to the existence of defects in the crystalline structure, called dislocations [START_REF] Landau | Theory of Elasticity[END_REF]2,[START_REF] Kittel | Physique de l'état solide[END_REF], which are set into motion above some local threshold stress. Elementary motion steps constitute discrete relaxation events, which result in a drift of the reference state.

In foams and dense emulsions, the local arrangement of elementary objects (droplets, foam cells) can be highly disordered. In the framework of crystals, this corresponds to the limit of a very high concentration of dislocations. Hence, a small increment of stress gives generally access to a large number of relaxation events. This limits the accessible range of purely elastic deformations between successive elementary relaxations.

Like in crystals, each relaxation event is associated with a topological flip: the stucture of the network locally changes. In foams and emulsions, such events are known as "T 1 processes": nearest neighbour links are exchanged between two pairs of cells. This class of materials thus exhibits an original interplay of elasticity (geometry and continuity), and plasticity (topology and discreteness).

B. Immortal vs. aging fluids

Foams and emulsions usually undergo spontaneous evolution such as coarsening (due to coalescence or ripening) or drainage [START_REF] Weaire | The Physics of Foams[END_REF][START_REF] Coussot | Rhéophysique des pâtes et des suspensions[END_REF][START_REF] Coussot | Rheometry of pastes, suspensions and granular materials[END_REF]. Such changes may induce a few topological rearrangements and cause the material rheological properties to slowly evolve [START_REF] Cipeletti | Slow dynamics in glasses, gels and foams[END_REF][START_REF] Michel Cloitre | Rheological aging and rejuvenation in microgel pastes[END_REF] -a phenomenon known as aging. If the material is subjected to a weak external stress (far below the plastic threshold), such rearrangements may also statistically induce some creep which would otherwise not occur [START_REF] Cohen-Addad | Origin of the slow linear viscoelastic response of aqueous foams[END_REF].

In other materials made of soft, round objects, the relevant molecular processes are slow and no subtantial aging is observed. Among such materials -which may be called 'immortal' -are copolymer micelles [START_REF] Eiser | Flow in micellar cubic crystals[END_REF][START_REF] Eiser | Nonhomogeneous textures and banded flow in a soft cubic phase under shear[END_REF] and some foams and emulsions formulated in such a way as to make the diffusion of the dispersed phase and the rate of film ruptures imperceptible within the experimental time-scale.

C. A brief history of flow localization

The structural characteristics summarized above lead to interesting non-linear mechanical behaviours in which a rich interplay exists between structural and mechanical responses.

One of the most extensively studied problems concerns flow localization, studied in various materials, from micellar solutions to granular flows.

In the thoroughly studied system of surfactant solutions self-organized as giant micelles, the flow curve was originally observed to exhibit a plateau-like behaviour under controlled shear rate [START_REF] Berret | Isotropic-to-nematic transition in wormlike micelles under shear[END_REF].

Structural observation followed, demonstrating shear-banding. Thus, in the stress plateau region, two different organizations of the material coexist: an isotropic region, similar in structure to the original solution, and a strongly birefringent region, in which the micelles are aligned to a high degree with the flow direction [START_REF] Decruppe | Insight in shear banding under transient flow[END_REF].

This situation was initially interpretated in terms of out-of-equilibrium phase transitions in the material, leading to a steady-state coexistence between two structurally homogeneous domains [START_REF] Porte | Inhomogeneous flows of complex fluids: Mechanical instability versus non-equilibrium phase transition[END_REF].

Theoretical descriptions attempted to capture the onset of this transition and the stability of the coexistence in terms of interfacial dynamics and mechanical instability [START_REF] Lu | Effects of nonlocal stress on the determination of shear banding flow[END_REF].

A more detailed investigation of the birefringent phase has more recently revealed strong spatial [START_REF] Lerouge | Interface instability in shear-banding flow[END_REF] and temporal [START_REF] Bécu | Spatiotemporal dynamics of wormlike micelles under shear[END_REF] variations. That disagrees with the initial simple picture.

The birefringent phase consists in numerous transient, narrow zones of very large shear.

The term "fluid fracture" [START_REF] Berret | Evidence of shear-induced fluid fracture in telechelic polymer networks[END_REF] has been proposed to describe these individual events, which have been observed in different systems with similar rheological properties but different internal structures (connected microemulsions [START_REF] Molino | A transient network of telechelic polymers and microspheres: structure and rheology[END_REF], copolymer cubic-phase [START_REF] Eiser | Flow in micellar cubic crystals[END_REF][START_REF] Eiser | Nonhomogeneous textures and banded flow in a soft cubic phase under shear[END_REF]).

The understanding of the shear-banding phenomenology has thus begun to shift from a phase transition picture to a fracture picture, and the interest is now focusing on the highly localized events that initiate the transition, and on their relation to the structural properties of the material.

Similarly, in the case of foams, a similar phenomenology of localization has been observed both in experiments [START_REF] Debrégeas | Deformation and flow of a two-dimensional foam under continuous shear[END_REF] and simulation [START_REF] Kabla | Local stress relaxation and shear banding in a dry foam under shear[END_REF], and interpreted [START_REF]by reflections on the dynamics of earthquakes[END_REF] in terms of the interaction between individual events via (mostly elastic) deformations of the surrounding material [START_REF] Falk | Dynamics of viscoplastic deformation in amorphous solids[END_REF][START_REF] Picard | Elastic consequences of a single plastic event : a step towards the microscopic modelling of the flow of yield stress fluids[END_REF][START_REF] Picard | Slow flows of yield stress fluids: complex spatiotemporal behavior within a simple elastoplastic model[END_REF]. The challenge now consists in understanding the self-organization of dispersed relaxation events into a localized fracture-like behaviour.

Thus, the emergence of a fluid fracture from these discrete, elementary relaxation events appears as a well-defined problem. In that respect, two main problems remain open: (i) the role of structural disorder, and (ii) non-local effects between individual events, mediated through elastic stresses in the material Both problems relate to the unknown typical length scale on which the discrete system should be averaged for a descriptions in terms of a continuous model. More precisely, problem (i) addresses this length scale from the limit of smaller length scales where disorder is relevant. Conversely, problem (ii) addresses it from the larger length scale limit: a suitable constitutive equation, incroporated into the general framework of continuum mechanics, provides the tools for generating such non-local effects.

D. Ingredients of our model

In the present work, we focus on the second problem discussed above, and construct an example of a rheological model inspired by the behaviour of such "immortal" fluids. It is characterized by four main features.

1. The flow properties are motivated and discussed in terms of microscopic considerations (T 1 processes).

2. In order to incorporate the non-local elastic effects mentioned above, our model is written in a fully tensorial form, whether in two and in three dimensions.

3. It is a commonly observed feature that such soft materials deform substantially before they display plasticity. In other words, their yield stress is comparable to their elastic modulus (unlike for classical, hard crystals). Correspondingly, the present model implements a consistent description of the elastic properties of the material up to finite deformations (i.e., beyond the usual approximations valid at small deformations).

4. As we shall see now, an amorphous, elastic material undergoing plasticity loses the memory of past events. The initial reference state thus has no physical relevance. Correspondingly, our model is developed in the Eulerian formalism (attached to the current reference state).

II. LOSS OF MEMORY AND CONSEQUENCES

FIG. 1: When a block of foam undergoes very large deformations (symbolized by the evolution of the coordinate system between the three left-hand side drawings), one might think that its local structure correspondingly evolves towards a highly stretched configuration. In fact, it always remains similar to its initial state, as illustrated on the right-hand side drawing.

In foams or emulsions and in crystals alike, large deformations of the sample do not imply large deformations of individual objects, since topological rearrangements relax local stresses. In both types of systems, two objects that are initially in contact can end up at a large mutual distance when many topological rearrangements have occurred. In such a situation, there does not exist any kind of elastic restoring force between both objects.

Hence, the distance between them is irrelevant to the current mechanical behaviour. As a consequence, the deformation from the initial state, which keeps track of such large distances, is mechanically irrelevant, even though it has a clear experimental meaning. In other words, the material has lost the memory of such large deformations.

As we shall see, this is the reason why:

• we use Eulerian coordinates;

• we define a quantity called "stored deformation";

• we discuss the impact of local disorder.

A. Eulerian description FIG. 2: Coordinate systems for a deformed foam. Let us attach a system of coordinates to the initial state of a foam (left-hand side drawing). After some deformation, we may describe the quantities attached to the foam in terms of either of two coordinate systems (right-hand side drawings): either the initial (now deformed) coordinate system, or a new (undistorted) coordinate system defined on the current state of the foam. In the case of elastic deformations, when material keeps trace of its initial configuration, such a choice does not have noticeable consequences. But in the case of deformations that imply plastic events which progressively erase the memory of the initial state (see Figure 1), the choice of the initial, much deformed coordinate system would not be physically (or computationally!) particularly relevant.

In order to describe the deformations of a material sample and the evolution of the physical quantities attached to it, two types of coordinate systems are commonly used: either "Lagrangian" coordinate systems attached to the initial state of the sample, or "Eulerian" systems attached to its current configuration, which coincide to first order in a small deformation context.

In the case of an elastic material, that keeps the memory of its initial state, the choice of one or the other does not have consequences other than computational (physicists most usually use Lagrangian coordinates [START_REF] Landau | Theory of Elasticity[END_REF]).

However, in the case of a material that progressively looses the memory of its initial configuration, such as a foam or an emulsion via rearrangements (T 1 processes described later in this paper), it would be physically irrelevant (and computationally tedious) to refer to the initial sample state. One therefore generally uses Eulerian coordinates in such situations.

For instance, the Navier-Stokes equations are usually expressed in a fixed coordinate system (Eulerian approach).

B. From deformation to 'stored deformation'

Experimentally, the accessible variables are (1) the deformation or deformation rate, as measured or imposed at the sample boundary, and (2) the stress (at least in some systems such as foams, where it can be extracted from the shape of the individual objects, or in photoelastic systems). In order to set up a spatio-temporal numerical scheme, one needs not only continuity and force balance equations, but also a specific evolution relation of the local stress in terms of the deformation rate. Our purpose in this paper is to propose a model example for this missing ingredient, and to explore its properties.

The very notion of a fixed reference state, and of a global deformation from this state, being of no use, the stress, as an index of the local elastic deformation, represents only the recoverable (or 'stored') part of the deformation. It can be defined through the following thought experiment, described on Fig. 3: a fragment of the material is cut, in the deformed state, and allowed to relax; the "stored deformation" is defined as the inverse of the deformation observed during this relaxation.

In the case of foams, it was shown a few years ago [START_REF] Aubouy | A texture tensor to quantify deformations[END_REF] that it is possible to construct a deformation tensor from the experimentally observed inter-bubble (centre-to-centre) vectors, which indeed faithfully represents the stress [START_REF] Asipauskas | A texture tensor to quantify deformations: the example of two-dimensional flowing foams[END_REF].

C. Disorder, stored deformation, and stress

The knowledge of the local stored deformation is exactly equivalent to the knowledge of the stress. They are related through a specific, material-dependent, constitutive relation, namely the elasticity, be it linear or not.

Generically, disordered systems are locally frustrated and contain internal stresses even in the absence of external applied stress. In other words, stored deformations are nonzero.

One might think it possible to relax stored deformations by cutting the material into pieces and sewing them together again, as described above. In fact, the relaxed pieces do not fit together nicely, even after adjusting local orientations: the field of stored deformations cannot be reconstructed from a displacement field. If one were to sew all pieces together again, one would need to stretch each of them appropriately, thus reconstructing a frustrated stress field corresponding to a state with zero external applied stress.

The notion of reference state is always clear in a local context. But it cannot be extended to any macroscopic part of a disordered material. Indeed, it would not be extensive: one half of a relaxed sample generally does not match the relaxed state of the same sample half.

D. Evolution of the stored deformation

Let us return to the construction of the local evolution of the stress in terms of the applied deformation rate.

The main point is that topological events participate in the applied deformation but relax part of the corresponding stored deformation. The evolution of the stored deformation must thus involve both an entrainment part and a plastic part.

The entrainment part is purely kinematic, driven by the velocity gradient. The plastic part reflects the T 1 relaxation processes triggered at large stored deformations. It always [45] tends to lower the stored deformation. It reflects the rate at which the material looses memory of the local reference state which is implicit in the stored deformation.

III. CHOICE OF A RHEOLOGICAL MODEL

Let us now choose a rheological model. The considerations of Section II D the rheological model of a foam must incorporate a spring (which represents elasticity) in series with a creeping, plastic part. In order to choose this plastic part, let us now review a few common rheological models. They correspond to the generic form given by Figure 4, where we added an optional viscous element in parallel with the other two elements. Note that such a viscous element impacts the stress response of the material, but not the local dynamics of stored deformation (at least not directly [46]).

Such models are listed on Table III together with an indication of their creep and relaxation properties. We now review some of these models, which have been used in the context of foams or similar materials.

A. Burger model for weak applied stresses

The rheology of dry liquid foams under weak stresses is well-described [START_REF] Cohen-Addad | Origin of the slow linear viscoelastic response of aqueous foams[END_REF] viscous element) corresponds to the spontaneous T 1 processes (which are responsible for aging phenomena) being slighlty biaised by the ambiant stress.

B. Models with non-zero threshold

Although well suited to describe the response of foams at weak stresses, the Burger model does not incorporate the finite plasticity threshold. Let us now review some models that do incorporate the threshold, even though they oversimplify the short-time response under weak stresses.

More precisely, these models are designed to provide most, if not all, of the following features: (i) a simple elastic response to small stresses, (ii) yielding above a stress threshold, and (iii) a viscous response at large, constant deformation rates.

The first family of such models includes the simple elasto-plastic model (a spring in series with a solid friction element), with an optional viscous element in parallel (Table III).

The simple elasto-plastic model does not provide feature (iii). The viscous element was incorporated by Marmottant and Graner [START_REF] Marmottant | An elastic, plastic, viscous model for slow shear of a liquid foam[END_REF] to account for this feature, observed in foams.

The second family, also used in the context of foams or rheologically similar materials (see Table III), includes a viscous element coupled to the solid friction element, i.e., the Bingham model (Takeshi and Sekimoto [START_REF] Takeshi | Internal stress in a model elastoplastic fluid[END_REF], and present work), also with an optional viscous element in parallel (Saramito [START_REF] Saramito | A new constitutive equation for elastoviscoplastic fluid flows[END_REF]).

In all cases except the pure elasto-plastic model, the viscous element provides feature (iii). assumed to respond linearly. The ability of the material to relax from above the threshold stress when applied deformation is stopped is indicated. The value of the stress threshold is compared to that of the elastic modulus in order to estimate the deformation at the threshold. Models have been labeled as "linear" when the deformation at the threshold is small, whether explicitely or implicitely. Some references are given when such models have been used in the context of foams or other soft, disordered materials. In some instances (marked with label "+ viscous"), an additional viscous element was introduced in parallel with the other elements altogether. This element is then indicated with dashed lines in the corresponding diagramme in the left-hand side column.

C. Each foam has its own rheological model

The reason why different models have been suggested is that the rheology of a foam varies with several parameters, among which surface tension, bubble size, polydispersity, surfactant properties. Let us concentrate on the effect of volume fraction and continuous phase viscosity, see Figure 5.

Elasticity and volume fraction

At gas volume fractions below the close-packing threshold, the foam flows and displays no elasticity.

Slightly above the threshold, elasticity is weak at small stored deformations as bubbles move rather freely between neighbouring bubbles, and strengthens at larger stored deformations as they come into closer contact with neighbours.

At large gas volume fractions, the elastic modulus is expected to be large even at small deformations, as the bubbles are already in close contact.

Plastic threshold and volume fraction

In order for neighbouring bubbles to undergo a plastic rearrangement (such as a T 1 process), they need to deform more importantly when the gas volume fraction is larger. As a result, the plastic threshold σ y is also expected to be larger.

Large elastic deformations

The plasticity threshold of foams usually corresponds to moderate deformations, that are beyond the small deformation regime. For instance, the threshold deformation for a polydisperse foam under shear in the dry limit (volume fraction approaching unity) is on the order of 30 to 50% [START_REF] Kraynik | Microrheology of random polydisperse foam[END_REF].

The ratio σ y /G is an indication of whether elastic non-linearities appear prior to the onset of plasticity. It is not clear to us whether this ratio increases or decreases with gas volume fraction.

Relaxation and volume fraction

Under stationary conditions, elasticity is inactive and both viscous elements play similar roles. They can be distinguished, however, in transient responses.

At moderate volume fractions, we expect any relative motion of bubbles to generate viscous dissipation: the general viscous element should dominate over the viscous friction element.

At higher gas volume fractions, bubbles interact so intimately that most viscous dissipation can be expected to arise during plastic events. In other words, the viscous friction element should dominate over the general viscous element.

The relative weight of the viscous friction element and the general viscous element is apparent on Figure 5.

Their ratio also impacts the ability of the material to relax when the applied deformation rate is suddenly brought to zero. Indeed, in this respect, only the viscous friction element is relevant. A model without such a viscous friction element [START_REF] Marmottant | An elastic, plastic, viscous model for slow shear of a liquid foam[END_REF] does not display relaxation in common situations (like oscillatory measurements) where the deformation rate may reverse.

Viscous vs. plastic behaviour

Besides, the visosity of the continuous phase impacts the relative importance of the viscous elements and the solid friction element. This impacts the stationary response of a foam, which typically changes from mainly plastic to mainly viscous as the deformation rate is increased. This transition is expected to occur at lower deformation rates if the continuous phase viscosity is increased.

D. Choice of the Bingham model

Let us now choose a specific model in order to develop a fully tensorial version of it.

Except for sollicitations at very low deformation rates, the viscous elements impact the rheological response. Should we keep them both?

The general viscous element provides an additional contribution to the stress. It does not represent any technical difficulty. Besides, it does not impact the behaviour of the system under imposed applied deformation conditions.

By contrast, we believe that the viscous friction element provides essential features such as the ability to display relaxation.

For these reasons, in the remaining part of this paper, we focus on the Bingham model

Volume fraction

Viscosity of continuous phase (see Fig. 6), which is the simplest one to provide all three desired properties together with relaxation. We hope that it may also apply to a broad range of materials made of densely packed, soft, essentially round objects.

We shall keep in mind that the parameters of the model (G, σ y and η on Fig. 6) will depend on such physical quantities as the volume fraction and the viscosity of the continuous phase, as discussed in Section III C and illustrated by Fig. 5. Let us now check (with its scalar version) that the Bingham model (Figure 6) behaves as expected for this class of materials. In this model, the evolution of stress (for positive values of the stress [START_REF]) values of the stress, Eq. (1)[END_REF]) is given by:

σ = G D - (σ -σ y ) η/G θ[σ -σ y ] (1)
where θ(x) = 1 when x ≥ 0 and θ(x) = 0 otherwise, and where D is the applied deformation rate, i.e., the symmetric part of the velocity gradient tensor:

D = ∇ v + ∇ v T 2 (2)
Equivalently, the evolution of the spring elongation ε = σ/G is given by:

ε = D - (ε -ε y ) η/G θ[ε -ε y ] (3)
where ε y = σ y /G.

Let us now consider successively three simple experiments: (a) quasistatic imposed deformation; (b) constant, imposed deformation rate; (c) constant, imposed stress. The corresponding evolution of the main rheological variables is schematically represented on Figure 7.

Quasistatic imposed deformation

Under low applied deformation, the stress depends linearly on deformation. This low deformation regime is valid as long as the resulting stress is smaller than the yield value σ y .

At larger deformations, the stress remains constant and equal to σ y , even under arbitrary large (but constant) deformations.

The reason why foams and emulsions display such a solid friction behaviour is that the T 1 processes are very similar to the relaxation of surface bumps involved in the friction between rough solids: in both cases, part of the mechanical work done by the imposed stress is dissipated in discrete relaxation events which enable discrete deformation steps. As a result of these events, the work is proportional to the total deformation (rather than to the velocity, as in a viscous fluid).

Under imposed deformation, supposing that the deformation value is reached through a quasi-stationnary process, the stress will remain constant as soon as the elastic stress associated with a deformation increment is exactly compensated by the stress relaxed through the plastic processes. In this situation, even if the imposed deformation can be arbitrarily large, the stored deformation remains equal to the value that corresponds to the yield stress: any extra applied deformation is relaxed through T 1 processes. The noisy aspect of the stress plateau reflects the disorder of the material -and hence, of the distribution of available relaxation processes.

Constant imposed deformation rate

Under constant applied deformation rate (see Figure 7b), the stress rises linearly at short times, as long as it is smaller than σ y . At later times, it eventually stabilizes above the threshold, and its final value σ ∞ increases with the deformation rate (affinely in the Bingham model).

Quasistatic imposed stress

Under constant applied stress, by contrast (see Figure 7c), the system displays two different behaviours. At small stress values (σ < σ y ), it behaves elastically. It switches to a flow behaviour at larger stress values (σ > σ y ).

( Note that the gas diffusion mode and the liquid permeation mode are plastic, even for small magnitudes: they are accompanied by dissipation and they lead to situations that are stretching deviatoric isotropic stretching gas diffusion permeation liquid

( stable even in the absence of any extra applied stress. Indeed, in situation [START_REF] Landau | Theory of Elasticity[END_REF], the extra amount of gas occupies the extra volume. Similarly, in situation (4), the amount of gas in the bubble and the bubble volume have remained constant, so the force balance within the material has not been altered despite the loss of liquid.

In the present work, for the sake of simplicity, we restrict to changes in the applied stress i.e., modes (0) → (2) and (0) → (3). For small magnitudes, these two modes are elastic: if the applied stress is removed, the system returns to state (0). For larger amplitudes mode (0) → (3) is plastic. Indeed, under large applied deviatoric stresses, rearrangements of the liquid films between bubbles (like the T 1 processes described in Section VI A below) lead to stress-free states that differ from the initial state, even though they may locally be very similar, if not identical, to state (0).

Density, velocity and stress

With this choice of modes (0) → (2) and (0) → (3) rather than modes (0) → (1) or (0) → (4), the weight fraction of each phase remains constant, and the material can safely be handled like a one-phase fluid, which substantially simplifies its description.

The overall material density, ρ, can be expressed in terms of the contributions from both phases:

ρ = ϕ ρ liq + (1 -ϕ) ρ gas (4) 1 ρ = ϕ w ρ liq + 1 -ϕ w ρ gas (5)
where ϕ (resp., ϕ w ) is the volume fraction (resp., weight fraction) of the liquid phase.

Because we ignore the gas diffusion and liquid permeation modes depicted on Fig. 8, the velocity of the material at larger length scales suffers no ambiguity, as it is the same in the dispersed phase (bubbles or droplets) and in the continuous (liquid) phase:

v = v dispersed phase = v continuous phase (6)
As for the stress, it varies strongly at the microscopic scale within such a structured medium as a foam: compressive within a bubble and across the gas/liquid interfaces, tensile along these interfaces (surface tension), tensile within the liquid Plateau borders and vertices.

In the present work, the stress variable σ represents the sum of these contributions averaged at some larger length scale, where such structural details are smoothed out:

σ = σ gas + σ gas/liq interf. + σ liq (7)

G. Compressibility

In Section III D, we motivated our choice of the Bingham model to describe the rheological behaviour of an emulsion or a foam. The above considerations, however, were based essentially on scalar arguments.

One of the main properties that reflect the actual, tensorial nature of deformations and stresses in the material is compressibility, which is the ability of the material to adapt its volume when the pressure (i.e., the isotropic part of the stress) is changed (see mode (0) →

(2) in Fig. 8). A material is considered incompressible when the elastic modulus involved for changes in volume is much larger than the elastic modulus involved for deformations at fixed volume. In practice, the deformations of such a material therefore obey a fixed volume condition (mode (0) → (3) in Fig. 8).

Among the materials we address in the present work, emulsions can be considered incompressible for all practical purposes since both phases are liquid.

Foams can also be considered incompressible as long as their constitutive bubbles are not too small. Indeed, their compression modulus is then typically equal to the pressure in the dispersed, gas phase, while their shear modulus (and other moduli corresponding to deformations at constant volume) are on the order of the interphase surface tension divided by the typical bubble radius. Hence, under atmospheric pressure and with usual liquids, gases and surfactant molecules, foams with bubbles not smaller than 0.1 mm in size can safely be considered incompressible. By contrast, the compressibility of foams made of micron-sized (or even smaller) bubbles cannot be neglected.

As for bubble monolayers, when regarded as two-dimensional foams, their apparent compressibility depends on the boundary conditions. When such a monolayer is squeezed between two solid plates, it can be considered incompressible under the same conditions concerning the typical bubble size as a three-dimensional foam. In other situations, a bubble monolayer may have at least one free interface, for instance when it floats on a bath of the liquid, continuous phase, and/or when its upper surface is in contact with the atmosphere.

In such a situation, each bubble is free to slightly deform in the vertical direction in order to better accomodate in-plane stresses. As a result, the bubble monolayer may appear compressible as seen from above, even though the total bubble volume in fact remains essentially constant.

In the present work, in order to be able to describe the rheology of all such systems, we consider a compressible material (modes (0) → (2) and (0) → (3) both available in Fig. 8). Nevertheless, specific properties or mathematical formulations suitable for the incompressible case (mode (0) → (3) only) are provided whenever appropriate.

IV. GENERAL FORMULATION FOR MATERIALS CAPABLE OF CREEP A. Evolution of the stored deformation

After the scalar description given in Section III E we now turn to a tensorial version of the model. In particular, instead of the scalar deformation ε, we will now use a deformation tensor e, to be defined in detail in Section V.

As mentioned at the beginning of Section III, we are interested in materials that display some elasticity and are capable of creep (see Figure 4). The evolution of their deformation can be decomposed into an elastic part and a creeping part [START_REF] Mora | Étude de l'intégration temporelle du tenseur taux de déformation. Application à la modélisation de l'élastoplasticité en grandes transformations[END_REF]. An example of such a decomposition was provided above for the (scalar) Bingham model, see Eq. (3). More generally, the evolution of the stored deformation can be written in the form:

ė = kinematics(e, ∇ v) -creep(e, D), (8)
where D is the applied deformation rate as defined by Eq. ( 2).

In the Bingham model, the creep term depends only on the stored deformation e. In other models [START_REF] Marmottant | An elastic, plastic, viscous model for slow shear of a liquid foam[END_REF], it additionally depends on D. This point will be discussed later, in Section X B.

B. Stress evolution

For a material capable of creep (see Section IV A), once the evolution of the stored deformation is known, the stress evolves as prescribed by elasticity (see Figure 9).

Hence, in the elastic regime, since stored deformation evolves purely kinematically, the stress evolution results from the kinematics and elasticity alone. In the presence of plasticity, although the evolution of the stored deformation includes an additional term (see Eq. 8), the stress evolution is still deduced therefrom in the same way, namely via elasticity.

C. Consistency of some commonly used stress evolutions

It is common habit to express the stress evolution directly in the form

dσ dt = g(σ, ∇ v), (9)
At first sight, this might seem exactly equivalent to the evolution for the stored deformation as given by Eq. ( 8).

In fact, if no special care is taken, such an expression does not usually correspond to the composition of the stored deformation evolution with elasticity as illustrated on Figure 9.

The validity of such an evolution equation for stress is then implicitely restricted to the domain of small stored deformations, where elasticity is linear.

As an example, as shown in another work [START_REF] Benito | [END_REF], generalized Maxwell models that involve the Gordon-Schowalter derivative interpolation [START_REF] Saramito | A new constitutive equation for elastoviscoplastic fluid flows[END_REF][START_REF] Gordon | [END_REF] suffers such restrictions (except for the special cases of upper and lower-convected derivatives).

V. RELAXED, REFERENCE STATE, AND ELASTICITY

x X F t+dt F t ∇ v dt FIG.
10: Variation of the local material stored deformation (see Fig. 3) that corresponds to an infinitesimal displacement field. As described by Eq. ( 11), tensor F describes the transformation from the relaxed, reference state (sphere, generic point X) to the current state (ellipsoid, generic point x). F is chosen to be a pure deformation, while ∇ v dt may include some rotation. In this example, the evolution of tensor F is purely kinematic (see Appendix A), as the material is supposed to behave in a reversible, elastic manner here. See Fig. 29 for details on the relaxed state.

Let us now construct a consistent framework for the material elasticity, with the help of Figure 10, which pictures the evolution of the stored deformation due to an applied flow.

A. Stored deformation

In such a plastic material, as pointed out above, there is no point to define the deformation with respect to some remote, initial reference state. At any time, however, every region of the material is stretched (or not). In other words, its conformation differs (or not) from the conformation it would have in the absence of stress from the neighbouring regions: the stress locally induces (via its elasticity) a stored deformation, which can be visualized in the form of an ellipse (or, more generally, of an ellipsoid). To obtain the ellipse (see Figure 3), one needs to cut out a piece of material, draw a circle on it while it is thus relaxed, and put it back in place.

The circle (or, more generally, the sphere) in the relaxed state is described by:

X T • X = R 2 , ( 10
)
where X is a vector whose origin is the center of the sphere, and whose end is a generic point on the sphere (see Figure 10). The ellipsoid in the stretched material can be described by some vector x, which depends linearly on X if the sphere radius R is infinitesimal, and can be expressed in terms of the transformation from the relaxed state to the current state:

x = F • X (11)
The transformation from X to x (see Figure 10) can be chosen as a pure deformation (see Fig. 3). Tensor F is thus symmetric, and in its principal axes, we have:

F =      λ 1 0 0 0 λ 2 0 0 0 λ 3      (12)
where, in the case of an incompressible material, det

F = λ 1 λ 2 λ 3 = 1.
The equation for the ellipsoid is obtained from equations ( 10) and ( 11):

x T • F -2 • x = R 2 (13)

B. Finger tensor and associated deformation

As mentioned in the caption of Figure 3, the relaxed state local orientation can be chosen arbitrarily. This choice thus must not have any incidence on the material elasticity. We therefore need a variable that reflects local deformation without being sensitive to local orientation. A commonly used such variable is the Eulerian Finger tensor B = F • F T . In the present case where the local orientation is chosen such that F be a pure deformation (i.e., represented by a symmetric tensor),

B = F 2 =      β 1 0 0 0 β 2 0 0 0 β 3      (14)
where

β i = λ 2 i .
From the Finger tensor, we can also construct a deformation[48] e = 1 2 (B-I):

e =      1 2 (β 1 -1) 0 0 0 1 2 (β 2 -1) 0 0 0 1 2 (β 3 -1)      (15)
C. Evolution of the stored deformation

Let us now describe how a piece of material under stress is further deformed when an infinitesimal displacement field vdt is applied to the material. The corresponding deformation,

∇ v dt, where (∇ v) ij = ∂v i ∂v j , (16)
weakly deforms the ellipsoid (see Figure 10). As shown in Appendix A, one can express the equation for the deformed ellipsoid and derive the evolution equation for tensor B:

dB dt -∇ v • B -B • ∇ v T = 0 ( 17
)
where d/dt = ∂/∂t + v • ∇ is the particulate derivative and where the entire left-hand side is the upper-convected objective derivative of B.

The evolution of the associated deformation, defined by Eq. ( 15), also involves its upperconvected objective derivative:

de dt -∇ v • e -e • ∇ v T = D, ( 18
)
where D is the applied deformation rate (see Eq. 2).

Equations ( 17) and ( 18) provide the variation of the stored deformation in the case of a purely elastic behaviour, i.e., in the absence of any plasticity. We shall soon discuss the origin of plasticity and indicate how it may enter such evolution equations. But let us first give a clear formulation of how the stress relates to the stored deformation.

D. Elasticity

Since the material is supposed isotropic, the principal axes of the stress coincide with those of the Finger tensor, and the stress can be expressed as a linear combination[49] of three powers of B, for instance:

σ = a 0 I + a 1 B + a 2 B 2 (19)
where a 0 , a 1 and a 2 are scalar, isotropic functions of B (i.e., of its scalar invariants). Equivalently, the deviatoric (traceless) part of the stress

σ = dev(σ) = σ - I d tr(σ) = σ + pI (20)
(where p is the pressure and d is the dimension of space) can be expressed as

σ = a 1 dev(B) + a 2 dev(B 2 ) =      σ1 0 0 0 σ2 0 0 0 σ3      (21)
and the pressure is given by:

p = -a 0 - a 1 d tr(B) - a 2 d tr(B 2 ) ( 22
)
More specifically, we will assume that the material is hyperelastic, i.e., that the stress results from the differentiation of an elastic potential E(B), here defined as being an elastic energy per unit mass of the material. It can then be expressed [START_REF] Bertram | Elasticity and Plasticity of large deformations[END_REF] as:

σ = 2ρ dE dB • B (23) σ = 2ρ dev dE dB • B (24)
with:

ρ dE dB = a 0 2 B -1 + a 1 2 I + a 2 2 B ( 25
)
In the case of an incompressible material (see Section III G for examples), the stored deformation (expressed in terms of B or e) does not depend on pressure: it only depends on the deviatoric part σ of the stress. Conversely, the deviatoric stress σ can be expressed in terms of B or e:

σ = F (B) = G(e) (26)
In this incompressible limit, the pressure p (more precisely, its term a 0 ) varies very strongly when det B departs from unity. In practice, one usually considers that p is not known explicitely: in a practical situation, it must be obtained from the boundary conditions, while strictly enforcing the constraint det B = 1.

VI. PLASTICITY

Now that we have described the evolution of the stored deformation in the material (Eq. ( 17) or ( 18)) and the material elasticity (Section V D), let us turn to the description of its plastic properties.

As mentioned earlier, a plastic term must be added to Eq. ( 17) or [START_REF] Berret | Evidence of shear-induced fluid fracture in telechelic polymer networks[END_REF] in order to reflect the deformation due to topological events:

dB dt -∇ v • B -B • ∇ v T = -2 D B p ( 27
) de dt -∇ v • e -e • ∇ v T = D -D B p ( 28
)
This equation is the tensorial version of the decomposition of the applied deformation into a kinematic part and a creep part expressed by Eq. ( 8).

In the present section, we recall the microscopic origin and discuss the mathematical properties of the plastic term D B p in the evolution equation.

A. Threshold for a single T 1 process in two dimensions

Figure 11 provides a two-dimensional illustration of a T 1 process in an emulsion or a foam.

Four objects (hereafter called quadruplet) are initially organized in a diamond conformation.

As compared to the initial conformation at rest (point R 1 ), a weak applied stress induces an elastic deformation, symbolized by a continuous branch (thick line). Once the applied stress exceeds a certain value (symbolized by σ y on Fig. 11), the system undergoes a discrete flip while the four objects swap neighbours. This is the T 1 process, symbolized by a jump onto the other branch.

More precisely, let us suppose that the quadruplet is symmetrical like that on Fig. 11, with axes ξ and ζ. The stress orientation that is most favourable for the T 1 process to occur is that in which the stress axes are aligned with those of the quadruplet. Thus, Fig. 11 is drawn in terms of Σ ξξ and Σ ζζ , with the non-diagonal stress component Σ ξζ being equal to zero.

ζ ξ d ξ -d ζ C 2 R 1 R 2 C 1 σ y -σ y Σ ξξ -Σ ζζ
The condition for the quadruplet to remain on branch 1 of Figure 11 can be expressed in the form:

f eigen (Σ) ≤ 0 ( 29
)
where

f eigen (Σ) = Σ ζζ -Σ ξξ -σ y T 1 (30)
Here, σ y is a real, positive number. When the pressure is varied or when the stress is not aligned with the quadruplet (non-zero shear stress component, i.e., Σ ξζ = 0 in the present, 2d case), it more generally takes the form:

σ y T 1 = σ y T 1 (p, Σ ξζ ) (31)
σ y T 1 is expected to increase with the applied pressure p = -Σ ξξ -Σ ζζ . Indeed, increasing the pressure reduces the typical object size, thus enhancing surface tension effects. It should therefore increase the stress threshold in such surface tension sensitive systems as foams or emulsions.

In the case of an incompressible material, σ y T 1 depends only on Σ ξζ , not on p. More generally (in two and three dimensions), the stress threshold function f eigen depends only on the deviatoric part of the stress and can be written in the form:

f eigen (Σ) = g eigen [dev(Σ)] (32)
As schematized on Fig. 12, let us now discuss how the T 1-threshold discussed above impacts the plasticity threshold and the plastic flow.

B. Plasticity threshold

Instead of a single T 1-capable quadruplet (Fig. 12, left), let us now consider a mesoscopic element of material, containing quadruplets with many different orientations (Fig. 12, centre).

If we neglect the effects of disorder (see next paragraph for a short discussion), the plastic threshold for this chunk of material is reached when the most favourably oriented quadruplet reaches its own threshold:

f (σ) = max α f eigen (σ α )
where α runs over all orientations in space, and where σ α is the representation of tensor σ in the axes of a quadruplet oriented according to α. Note that the most sensitive quadruplets are those oriented along the stress: the most favourable orientation α coincides with that of the stress tensor (possibly up to some permutation of the axes).

In the case of the particular, two-dimensional threshold given by Eqns. ( 31) and ( 30), we obtain:

f (σ) = |σ (1) -σ (2) | -σ y (p) (33)
where σ (1) and σ (2) are the eigenvalues of the stress tensor and where p = -σ (1) -σ (2) is the pressure.

In the case of an incompressible material, the plastic threshold depends on the sole deviatoric part of the stress, like in Eq. ( 32):

f (σ) = max α g eigen [dev(σ α )] = g[dev(σ)] (34)
C. Disorder and plasticity threshold

The local disorder of soft object positions and interactions implies that the stress acting on a particular quadruplet slightly differs from the ambient mesoscopic stress, and the threshold of some of the quadruplets is lower than expected for a particular orientation of the applied stress.

As a result, the threshold value σ y in Eq. ( 33) above is in fact slightly lower than the microscopic value σ y T 1 in Eq. ( 30). Similarly, expression [START_REF] Gordon | [END_REF] for f (σ) is slightly overestimated: σ α should in fact be understood as the α-oriented representation of a locally disorder-enhanced version of the mesoscopic stress σ.

D. Plastic deformation rate

In a continuum model, the rate at which topological events occur is expressed as the plastic flow D B p and depends on the stress σ, or equivalently on the material deformation (tensor B) via its elasticity (see Eq. 24). The plastic flow therefore has the following form:

D B p (σ) = h(B), (35)
where h is a tensorial function. Since the plastic flow results from T 1 events, whose principal axes almost[50] coincide with those of the stress (see paragraphs VI A and VI B), function h is an isotropic function of B: it is such that the principal axes of D B p also coincide with those of B. Thus, in the same axes as those of B:

D B p =      δ 1 0 0 0 δ 2 0 0 0 δ 3      (36)
Equivalently, since powers of tensor B constitute a basis for symmetric tensors that have the same principal axes, D B p can be decomposed for instance in the following way:

D B p = b0 I + b1 B + b2 B 2 , (37)
where scalar functions b0 (B), b1 (B) and b2 (B) depend solely on the scalar invariants of tensor B. Let us emphasize the fact that in general, all three functions should be non-zero above the plasticity threshold. Restricting D B p to be proportional to B, for instance, would mean that all three directions would flow indepedently of each other. This would be a very specific choice, certainly not relevant for most materials. It would be the plastic equivalent to taking a Poisson ratio ν = 0 in linear elasticity, a property valid for only a restricted class of materials.

Function h (or equivalently eigenvalues δ i in Eq. 36 or coefficients b0 , b1 and b2 in Eq. 37), must also obey some constraints based on physical grounds, which are discussed in the next few paragraphs. Due to plasticity, the deformation between two real configurations (see Figure 13) does not reflect solely the increment dB dt dt in stored elastic deformation, but also the irreversible drift DB p dt of the relaxed local configuration, where the drift rate DB p is related to the plastic term D B p defined by Eq. ( 35) through covariant transport via tensor F :

E. Incompressible plastic flow rate

DB p = F -1 • D B p • F -1 (38)
In the present work, we do not consider all possible plastic deformation modes. Among the situations depicted on Fig. 8, conformations (0), ( 1) and ( 4) are equilibrium situations (at least on time scales where the dispersed phase does not diffuse and where the continuous phase does not permeate). Hence, modes (0) → (1) and (0) → (4) correspond to plastic evolutions of the material. In the present work, as stated in Section III F 1, such evolutions that involve mutual diffusion between both phases are not addressed.

The plastic evolutions of interest for us correspond to stress-induced evolutions, where the relaxed conformation is always characterized by the same amount of material per bubble or droplet (inside and around it). Among the equilibrium conformations (0), ( 1) and ( 4), only conformation (0) is thus eligible.

We are thus interested in elastic deformations (combinations of modes (0) → ( 2) and (0) → (3)) followed by plastic rearrangements that bring the local conformation back to (or at least towards) situation (0). In other words, the local relaxed conformation is always locally similar to situation (0), even though some bubbles or droplets may have swapped positions.

In particular, the evolution of the local relaxed state must not be accompanied by any change in volume. This can be expressed through a condition on the drift rate DB p : tr( DB p ) = 0 [START_REF] Goyon | Microfluidic velocimetry reveals spatial cooperativity in the flow of soft glassy materials[END_REF] Using Eq. ( 38), this can be expressed in terms of D B p :

0 = tr(B -1 • D B p ) (40) = δ 1 β 1 + δ 2 β 2 + δ 3 β 3
Expression (40) thus expresses incompressibility for the plastic flow rate D B p . Note that in order to satisfy the incompressibility condition given by Eq. ( 40), the plastic flow given by Eq. ( 37) can be written in the form:

D B p = b 1 B • dev(B) + b 2 B • dev(B 2 ) ( 41
)
Scalar functions b 1 and b 2 , like b0 , b1 and b2 in Eq. ( 37), depend solely on the invariants of tensor B.

In this case, the coefficients of Eq. ( 37) can be obtained through the relations:

b0 = b 2 (42) b1 = - b 2 2 [tr(B)] 2 + b 2 6 tr(B 2 ) - b 1 3 tr(B) (43) b2 = b 1 + b 2 tr(B) (44)
and are related to one another through: At this point, it may be useful to precisely delineate two different types of incompressibility, which can be discussed using Fig. 8.

As stated in paragraph VI E above, in the present work we are only interested in locally relaxed states of type (0), and the plastic flows we consider correspond to relaxed conformations that evolve among such states of type (0), and the density of the locally relaxed state is conserved: the material is plastically incompressible.

By contrast, as mentioned in paragraph III F 1, the actual state of the system (with local stresses) may have a local density that differs from that of the relaxed state, due to evolutions of type (0) → (2). In other words, the material is assumed to be compressible, that is: elastically compressible.

[51]

It follows that with our assumptions, the material may locally change volume (and the density ρ then departs from its initial value ρ 0 ) when the flow is not divergence-free (∇ • v = trD = 0). Nevertheless, because the plastic evolution is assumed incompressible, the material returns to its initial density ρ 0 as soon as the local stress vanishes.

For some systems, elastic compressibility is negligible and the material can be considered elastically incompressible, as explained in paragraph III G. The material then locally never changes volume, and the flow is then divergence-free (∇ • v = trD = 0).

G. Thermodynamic constraints

The total work developed by the stress in the material is given by tr(σ • D). In the present system, the internal energy U is purely elastic: U = E. The first principle of thermodynamics can thus be written as

ρ dE = tr(σ • D) dt + δQ, ( 46
)
where δQ is the heat uptake per unit volume.

Besides, the dissipated power P dissip per unit volume results from viscosity, especially during the relaxation of individual T 1 processes. It constitutes the only source of entropy in the system, which is non-negative according to the second principle:

P dissip = T dS created ≥ 0 ( 47
)
It also constitutes the only source of heat in the system. If we assume that the system remains at a constant temperature, this condition can be written as:

P dissip dt + δQ = 0 (48)
The Clausius-Duhem inequality is readily derived from the above equations:

tr(σ • D) -ρ dE dt ≥ 0 (49)
Let us express this inequality in a different way.

From Eq. ( 23), we get:

tr(σ • D) = 2ρ tr B • dE dB • D (50)
We also have:

dE dt = tr dE dB • dB dt (51)
where dB/dt is given by Eq. ( 27). Hence, from Eqs. ( 49), ( 50) and (51), we obtain the constraint on the plastic flow D B p that corresponds to the Clausius-Duhem inequality:

P dissip = 2ρ tr D B p • dE dB ≥ 0 ( 52
)
H. Plastic deformation rate: a summary

In the last few paragraphs, we showed that the plastic deformation rate must obey constraints that express the fact that:

dB dt -∇ v • B -B • ∇ v T = -2 D B p ( 53
)
One needs the elastic law (Eq. 19 or 23)

σ = 2 ρ dE dB • B (54)
and the usual force balance equation

∇ • σ + ρ f = ρ d v dt , ( 55
)
where f represents external forces per unit mass. The evolution of the density ρ obeys the usual mass conservation equation

∂ρ ∂t + ∇ • (ρ v) = dρ dt + ρ trD = 0 ( 56
)
We provide later (see Eq. 94) a version of this equation that includes liquid permeation.

Let us now mention the special case of an initially non-dilated material, and discuss dissipation at weak applied stresses.

A. Initially non-dilated material

Let us assume that the material initially has a uniform density i.e., with stored deformation

ρ(t 0 , r) = ρ 0 ( 57
)
and that it verifies everywhere det B(t 0 , r) = 1, (58)

In such a case, our assumptions stated in paragraph III F 1 imply that the material density is related to the determinant of tensor B in a very simply way. Indeed, the derivation shown in Appendix B implies that at any later time, throughout the sample:

ρ = ρ 0 √ det B ( 59
)
With such assumptions, Eq. ( 59) can therefore be used to replace ρ within Eq. ( 55), and Eq. ( 56) is then not useful any more.

B. Viscous losses under weak stresses

The Bingham model addressed in the present work reduces to a simple elastic system when subjected to weak stresses (below the yield stress). As a result, vibrations may be present in the material, at arbitrary high frequencies. Such vibrations may be undesirable. Not only do they make numerical simulations of the above system of equations more complicated, if not impossible, but they do not faithfully reflect the damping observed in real materials. This problem does not arise under large stresses, as the plastic flow rate introduces dissipation.

One of the simplest ways to introduce some dissipation at weak stresses is to add a viscous term to the stress. Eq. ( 54) thus becomes:

σ = 2 ρ dE dB • B + V (D, B) (60)
where the viscous term V is a symmetric tensor which depends linearly[52] on the applied deformation rate D. As the material is locally anisotropic due to the stored deformation, V additionally depends on B.

The general form of function V (D, B) is a sum of terms whose principal axes are those of tensor B:

l 1 (B) tr(l 2 (B) • D) (61)
and of terms that depend tensorially on D:

m 1 (B) • D • m 2 (B) + m 2 (B) • D • m 1 (B) (62)
Here, all functions l i and m i are isotropic scalar functions of tensor B.

If we neglect any effect of the material deformation, we may simply add a linear viscous term, and Eq. ( 54) then becomes:

σ = 2 ρ dE dB • B + η D (63)

VIII. HOMOGENEOUS CONSTANT SHEAR FLOW: ONE EXAMPLE

Let us now use the constitutive Eq. ( 27) to obtain the evolution of the system in a very common type of flow: a shear flow at a constant shear rate γ starting at time t = 0.

For non-homogeneous flows, the constitutive equation must be combined with the classical mechanics equations for continuum media, as mentioned in Section VII. Here, for simplicity, we assume that the material is homogeneous and incompressible (uniform density ρ 0 at all times) and that the flow remains homogeneous: no shear-banding, etc.

A typical result from Eq. ( 27), coupled with an elastic law (here incompressible, see Eq. 21) is the shear stress as a function of the shear deformation γ since time t = 0 when shear started. An example of such a mechanical response (with parameters chosen as described in Section VIII B) is provided on Fig. 14. The equations of our model provide:

1. the material response in its elastic state (first part of the curves on Fig. 14);

2. the threshold that marks the onset of plasticity (point where curves split apart depending on the shear rate);

3. the transient response that results from plasticity (second part of the curves);

4. the stationary response, as a function of the applied shear rate γ (see Fig. 16).

A. Method

Let us recall the evolution equation ( 27) for the material deformation (tensor B):

dB dt -∇ v • B -B • ∇ v T = -2 D B p (64)
This equation must be understood as written in a basis attached to the (fixed) laboratory frame. Let x be the axis of velocity, y the axis of the velocity gradient and z the vorticity axis. The symmetry of the shear flow implies that the corresponding material deformation in Section VIII B. The response increases during the initial period of time, when the material deforms purely elastically (as expected, the response is then independent of the shear rate). The increase then levels off once the plasticity threshold has been reached, and saturates at a stationary value. For a plot of the final, stationary value on the applied shear rate, see Fig. 16.

B has two principal axes within the xy plane (in directions X and Y yet to be determined, see Fig. 17) and one principal axis along z. Note that D B p and σ also has the same principal axes X, Y and z as B (see Eqs. 23 and 41).

As shown in Appendix C 1, Eq. (64) provides a system of differential equations for β 1 , β 2 and for angle θ between axes x and y (see Fig. 17):

u = cos(2θ) (65) dβ 1 dt = γβ 1 √ 1 -u 2 -2δ 1 (β 1 , β 2 ) ( 66
)
dβ 2 dt = -γβ 2 √ 1 -u 2 -2δ 2 (β 1 , β 2 ) (67) du dt = γ√ 1 -u 2 1 -u β 1 + β 2 β 1 -β 2 (68)
β 3 = 1 β 1 β 2 (69)
where Eq. ( 69) results from the assumed material incompresibility. Once the evolution of tensor B is known, the stress is obtained through Eq. ( 24).

As for dissipation, given by Eq. ( 52), it involves both the gradient of the elastic energy (see Eq. 25) and the plastic deformation rate given by Eq. ( 41).

The dissipation per unit volume can now be expressed as:

P dissip = 2ρ tr D B p • dE dB (70)

When D B

p is known explicitely, it can be calculated more directly:

P dissip = (a 1 + a 2 β 1 ) δ 1 + (a 1 + a 2 β 2 ) δ 2 + (a 1 + a 2 β 3 ) δ 3 (72)
Figure 15 displays both the injected power, σ12 γ, and the dissipated power, P dissip , in the course of a shear experiment. Notice that as long as the material remains in the elastic regime, no dissipation occurs. The injected power is being stored entirely as elastic energy.

As dissipation starts, an overshoot of dissipated power is observed. Asymptotically, both quantities converge towards the same value while the elastic deformation and energy of the material reach their stationary values.

B. Chosen elasticity, threshold and plastic flow

Let us choose a very common (Mooney-Rivlin) type of incompressible elasticity, which has been shown to suitably approximate the non-linear elasticity of foams [START_REF] Höhler | Rheology of liquid foam[END_REF][START_REF] Höhler | Constitutive equation to describe the nonlinear elastic response of aqueous foams and concentrated emulsions[END_REF]. The corresponding elastic energy [START_REF] Bertram | Elasticity and Plasticity of large deformations[END_REF] can be expressed as:

ρ 0 E = k 1 2 (I B -3) + k 2 2 (II B -3) (73)
where

I B = tr(B) (74)
II B = 1 2 [tr 2 (B) -tr(B 2 )] = tr(B -1 ) (75)
Correspondingly, the coefficients of Eq. ( 21) can be expressed as:

a 1 = k 1 + k 2 I B ( 76
)
a 2 = -k 2 (77)
As in refs. [START_REF] Höhler | Rheology of liquid foam[END_REF][START_REF] Höhler | Constitutive equation to describe the nonlinear elastic response of aqueous foams and concentrated emulsions[END_REF], we choose the values of k 1 and k 2 in terms of the shear modulus G as:

k 1 = 1 7 G (78)
k 2 = 6 7 G (79)
Note that in Section IX, we explore other values for k 1 and k 2 (still in the framework of a Mooney-Rivlin elasticity).

For the plastic deformation rate D B p , as expressed by Eq. ( 41) so as to obey plastic incompressibility, we choose the following coefficient values: This choice implies, in particular, that the plasticity threshold corresponds to the following condition:

b 1 = (trB -4) θ(trB -4) (80)
tr(B -4) = 0 (82)

We now discuss qualitatively the evolution of tensor B. We then show the effects of the material elasticity and discuss the corresponding rheological response.

C. Three-dimensional evolution of the material stored deformation

With the choices made in Paragraph VIII B for the material elasticity and plasticity, the evolution of the system given by Eqs. (65)-( 69) is depicted on Fig. [START_REF] Berret | Evidence of shear-induced fluid fracture in telechelic polymer networks[END_REF]. For three different values of the shear rate γ, it represents the trajectory of the material stored deformation as for β 1 and β 2 , successively in the elastic regime and in the plastic regime. The plasticity threshold, given by Eq. ( 82) is also represented, as well as the locus of the stationary states.

In the elastic regime, the form of Eqs. (C2)-(C3) imply that β 3 remains equal to its initial value, β 3 = 1, i.e., that the product β 1 β 2 remains equal to unity. This can indeed be checked from Eqs. (C11)-(C12), and it reflects the fact that as long as no plastic events have occured in the material, the absence of material deformation in direction z (planar shear in the x-y plane, see Eq. C1) implies that the stored deformation is not modified in direction z.

In other words, in the elastic regime, the fact that no plasticity occurs in direction z (i.e., δ 3 = 0) implies that in the same direction, the stored deformation remains constant, i.e.,

β 3 = const.
Conversely, in the stationary regime, β 3 remains constant (like most other quantities) and this imposes that no plasticity occurs in direction z, i.e., δ 3 = 0.

Meanwhile, in the transient regime, the stored deformation β 3 in direction z has evolved from its initial value towards its new, stationary value, despite the absence of any velocity gradient in this direction. This is made possible by the plastic events (δ 3 = 0) which allow internal relaxation within the material even in the absence of any flow in this direction.

D. Shear thinning

The stationary shear stress is represented on Fig. 16 as a function of the shear rate γ.

Above the yield stress, the stress increases when γ is increased, as expected. Notice that it increases in a subliminear way. We believe that the main reason for this shear-thinning behaviour could be the fact that the plastic flow rate we chose, as expressed in termes B, is cubic rather than linear for large stored deformations: where, again, θ(x) = 1 when x ≥ 0 and θ(x) = 0 otherwise.

D B p = (trB -4) θ(trB -4) B 2 - B 3 tr(B) (83

E. Normal stress differences

The stress tensor σ is obtained from the Finger tensor B through Eq. ( 21). In the case of a homogeneous shear flow (see above), it can be expressed not only in the XY z basis (Eq. 21) but also in the xyz basis associated with the flow:

σ =      σ11 σ12 0 σ12 σ22 0 0 0 σ3      =      c 2 σ1 + s 2 σ2 cs(σ 1 -σ2 ) 0 cs(σ 1 -σ2 ) s 2 σ1 + c 2 σ2 0 0 0 σ3      (84)
where s (resp., c) denotes the sine (resp., the cosine) of angle θ between axes x and X (see Fig. 17).

The first and second normal stress difference can be expressed as:

N 1 ≡ σ11 -σ22 = (σ 1 -σ2 )u (85)
N 2 ≡ σ22 -σ3 = 1 -u 2 (σ 1 -σ3 ) + 1 + u 2 (σ 2 -σ3 ) (86)
with σi given by Eq. ( 21):

σi

-σj = a 1 (β i -β j ) + a 2 (β 2 i -β 2 j ) (87)
Hence, for instance:

N 1 = u a 1 (β 1 -β 2 ) + a 2 (β 2 1 -β 2 2 ) (88) N 2 = a 1 2 β 1 + β 2 - 2 β 1 β 2 + u(β 2 -β 1 ) + a 2 2 β 2 1 + β 2 2 - 2 (β 1 β 2 ) 2 + u(β 2 2 -β 2 1 ) (89)
In the case of Mooney-Rivlin, coefficients a 1 and a 2 are given by Eq. ( 76) and the normal stress differences can be expressed as:

N 1 = u(β 1 -β 2 ) k 1 + k 2 β 1 β 2 (90) N 2 = k 1 2 β 1 + β 2 - 2 β 1 β 2 + u(β 2 -β 1 ) + k 2 2 2β 1 β 2 - β 1 + β 2 β 1 β 2 + u β 2 -β 1 β 1 β 2 (91)
F. Discussion of the stress response

Let us now comment briefly on the results presented on Fig. 21.

Normal stress differences N 1 and N 2 increase gently and monotonically to reach their stationary values.

The material is under traction in the direction of vorticity (σ 3 > 0), as well as in the direction of β 1 , while it is under compression in the direction of β 2 (not shown). This is consistent with the fact that it is stretched in the direction of vorticity (β 3 > 1).

The salient feature of these results is that the shear stress σ12 presents an overshoot during its transient, a behaviour which is intrinsically unstable in a homogeneous material, and could trigger flow localization. Such an overshoot for the shear stress has been observed in sheared foams [START_REF] Khan | Foam rheology: III. measurement of shear flow properties[END_REF].

IX. IMPACT OF ELASTICITY AND PLASTICITY ON THE STRESS RE-

SPONSE

Our dynamical equation was formulated in terms of tensor B (see Eq. 64). From the resulting evolution of B, when can then derive the stress evolution using the elastic law (Eq. [START_REF] Molino | A transient network of telechelic polymers and microspheres: structure and rheology[END_REF]. In this section, we explore how the choice for the elasticity and the plasticity impact the stress response.

A. Impact of elasticity

Note that in this paragraph the deformation history will be the same in all cases, and only the stress in the material will differ. Just like in Section VIII, we restrict ourselves to a Mooney-Rivlin type of elasticity, which is complex enough for illustrative purposes, but which has no specific properties with regards to our problem.

The elasticity of dry foams has been shown [START_REF] Höhler | Rheology of liquid foam[END_REF][START_REF] Höhler | Constitutive equation to describe the nonlinear elastic response of aqueous foams and concentrated emulsions[END_REF] to be well captured by such a Mooney-Rivlin elasticity, using parameter values k 1 = 1/7 and k 2 = 6/7. These values correspond to Figs. [START_REF] Porte | Inhomogeneous flows of complex fluids: Mechanical instability versus non-equilibrium phase transition[END_REF][START_REF] Lu | Effects of nonlocal stress on the determination of shear banding flow[END_REF][START_REF] Lerouge | Interface instability in shear-banding flow[END_REF][START_REF] Bécu | Spatiotemporal dynamics of wormlike micelles under shear[END_REF][START_REF] Berret | Evidence of shear-induced fluid fracture in telechelic polymer networks[END_REF][START_REF] Molino | A transient network of telechelic polymers and microspheres: structure and rheology[END_REF][START_REF] Debrégeas | Deformation and flow of a two-dimensional foam under continuous shear[END_REF][START_REF] Kabla | Local stress relaxation and shear banding in a dry foam under shear[END_REF][START_REF] Falk | Dynamics of viscoplastic deformation in amorphous solids[END_REF] where we took the modulus value G = 1.

Exploring different values for k 1 and k 2 (still with G = 1), we show in Fig. 23 the time evolution of shear stress, stress in the vorticity direction and normal stress differences. The mechanical behaviour remains similar in all these cases. In particular, the shear stress always presents an overshoot. This overshoot is more pronounced for k 1 = 1 and k 2 = 0. In all cases, the normal stress differences do not present any overshoot and reach monotonically their stationary values. As for the stress along the vorticity direction, It presents an undershoot when k 1 > k 2 , and its stationary value becomes negative at some point between the situation in the center of Fig. 23 (k 1 = k 2 = 1/2) and the situation on the right-hand side (k 1 = 1 and k 2 = 0).

B. Impact of plasticity

Let us now briefly explore the impact of the plastic deformation rate D B p (B) through the influence of parameter b 2 in Eq. [START_REF] Emile | Swelling of a single foam film under slipping[END_REF].

In Fig. 24, we take b 1 = (trB -4)θ(trB -4) like in Figs. 14-22. We then observe the (unspectacular) effect of choosing b 2 = b 1 rather than b 2 = 0. Note that in this case, as opposed to the role of elasticity in paragraph IX A, the time evolution of tensor B is affected by the choice of the plastic law.

C. Discussion

In paragraphs IX A and IX B above, we illustrated the fact that the stress response of the material can be altered through a change in the material elasticity and plasticity alone.

Our approach emphasizes the fact that the choice of the convective derivative is somewhat arbitrary. A material is often considered as rather "upper-convected" or rather "lowerconvected" because its stress evolution follows more closely the eponymous convective derivative. In fact, always choosing the upper-convected derivative (and tensor B for a measure of the deformation) would be equivalent, provided the elastic and plastic laws be suitably adjusted. A forthcoming paper [START_REF] Benito | [END_REF] will discuss in detail the consistency of commonly used rheological laws (see Section IV C above) and explore these issues more in depth.

X. CONCLUSION AND PERSPECTIVES

In order to describe the rheological behaviour of foams and emulsions (and potentially other materials), we have developped a continuous framework to describe the evolution of the stored deformation tensor. It includes elasticity, up to the large deformations commonly encountered in such systems, and plasticity.

Within this framework, we showed that it is possible to play with various expressions for the elasticity and the plasticity. We hope that it is thus possible to adequately describe the rheological behaviour of a large range of foams and emulsions with dispersed phase volume fraction close to unity.

In less concentrated foams and emulsions, volume fraction should be added as an extra field to account for osmotic compressibility of the dispersed phase. The issue whether the model may be adapted to describe usual, aging complex fluids, is left for future investigations.

As mentioned in paragraph III A, the Burger model (see Fig. 25) describes the response of dry foams under weak stresses adequately. By contrast, the Bingham model (used as a basis for our formulation, see Fig. 26) captures the large stress behaviour. It would be interesting to combine both models, as represented in Fig. 27. We did not elaborate on this combined model, as it is technically more complex than our Bingham-like model, and we did not want to focus on short timescales or weak stresses, where both models have differing behaviours.

A. Determining the elasticity from experiments

Beyond all these rheological and phenomenological models, it is important to make the connection with the microscopic scale.

From this point of view, the main issue is the definition of the deformation in the material.

In the present work, deformation is built on a thought experiment: we cut out a piece of material at sufficiently large a length scale for disorder to be smoothed out. By contrast, a way to define deformation was introduced by Aubouy et al. [START_REF] Aubouy | A texture tensor to quantify deformations[END_REF] for all systems in which the individual objects and their mutual contacts are experimentally accessible, such as 2D foams and emulsions (and 3D foams and emulsions when tomography will have become a routine technique for imaging such systems on rheologically relevant timescales).

They first construct a symmetric tensor M from the centre-to-centre vectors [START_REF] Graner | Statistical tools to characterize discrete rearranging patterns, in 2 or 3 dimensions: cellular materials, assemblies of particles[END_REF] for pairs of first neighbours. This tensor is a dilation (proportional to the unit tensor when the material is at rest). They then define the deformation U as the logarithm of M.

One might think that their tensor U is just one deformation among many others. But with respect to the particular systems that they consider (2D foams in the dry limit under quasistatic deformation), it plays a very special role: it is the deformation for which the elastic law remains linear up to large deformations [START_REF] Asipauskas | A texture tensor to quantify deformations: the example of two-dimensional flowing foams[END_REF] (in fact, up to the onset of plasticity).

This result confers some weight to the initial, scale-independence arguments [START_REF] Aubouy | A texture tensor to quantify deformations[END_REF] for the choice of the logarithm as the link between their texture (dilation) tensor M and their deformation measure U = log M.

In our perspective, as we will show in a more elaborate manner elsewhere [START_REF] Benito | [END_REF], the choice of a deformation measure corresponds to a choice for the convective derivative. The deformation measure U [START_REF] Aubouy | A texture tensor to quantify deformations[END_REF][START_REF] Graner | Statistical tools to characterize discrete rearranging patterns, in 2 or 3 dimensions: cellular materials, assemblies of particles[END_REF] could therefore lead directly to a continuum formalism well-suited for (at least) two-dimensional foams in the elastic regime.

B. Plasticity and mechanical noise

In the present work, the plastic deformation rate D B p was assumed to depend only on the local stress via the local stored deformation, see Eq. ( 35). This dependence is sufficient to account for the non-local elastic effects observed in foams [START_REF] Debrégeas | Deformation and flow of a two-dimensional foam under continuous shear[END_REF][START_REF] Kabla | Local stress relaxation and shear banding in a dry foam under shear[END_REF], and mentioned in paragraph I C. Indeed, when plastic events occur (non-zero D B p ), this impacts the local stored deformation via Eq. ( 27). This, in turn, affects the stress due to elasticity (Eq. 23). Equation (55) then implies that the stress is modified in the surrounding material. Thus, within this continuum model in which the plastic deformation rate depends solely on the stored deformation (or stress), plasticity at one location alters the stress (and the stored deformation) elsewhere in the material, thus possibly contributing to triggering plasticity there.

Nevertheless, stress may not be the only factor that determines the rate at which T 1 processes occur. For instance, in a recent work by Marmottant and Graner [START_REF] Marmottant | An elastic, plastic, viscous model for slow shear of a liquid foam[END_REF], the plastic deformation rate (D B p in our notation) is proportional to the total deformation rate D when it has the same sign as the stored deformation [START_REF]a formulation specific to 2D materials[END_REF]. Incidentally, such a choice implies that relaxation[55] cannot take place in the system.

The fact that D B p depends on D can be interpreted physically as the fact that mechanical noise may well help triggering plastic events. As an alternative implementation of the impact of mechanical noise, one might take the plastic deformation rate D B p as slightly enhanced in the presence of non-zero total deformation rate (although mainly determined by the stored deformation). For instance, if we assume that the mechanical noise is isotropic, we may include a multiplicative factor of the form (1 + |D 2 | τ 2 ) into D B p , where plasticity would be significantly enhanced for deformation rates around and above τ -1 .

We did not enter such subtelties in the present work, and restricted our study to a purely stress-dependent plastic deformation rate, as expressed by Eq. [START_REF] Höhler | Rheology of liquid foam[END_REF].

C. Plasticity and fluidity or fragility

In the present model, the evolution of the stored deformation is given by Eq. ( 27

In reality, the plastic flow rate D B p may well depend on yet other variables than B and D.

A few years ago, in their study of sheared two-dimensional foams between two solid plates, Kabla and Debrégeas noticed that T 1 events appeared preferentially in regions where the stress tensor was most disordered [START_REF] Kabla | Local stress relaxation and shear banding in a dry foam under shear[END_REF].

Very recently, Goyon and coworkers [START_REF] Goyon | Microfluidic velocimetry reveals spatial cooperativity in the flow of soft glassy materials[END_REF] studied the flow of three-dimensional emulsions in milli-fluidic channels. Because of the large aspect ratio of the channel section, they were able to minimize edge effects and to obtain two-dimensional velocity profiles for the stationary flow of such emulsions between two walls. One of their key results is that even though they are able to measure locally the relation between the shear stress (which they derive from the applied pressure) and the shear rate (which they observe), they do not obtain a single mastercurve σ = σ( γ) when changing the applied pressure.

By contrast, Eq. ( 27), combined with Eq. (93), would predict, in stationary flow, that the stored deformation could be expressed in terms of the velocity gradient, B = B(∇ v), and similarly for the stress, via elasticity. Hence, in the present situation, they would predict:

σ = σ( γ)
The observed behaviour [START_REF] Goyon | Microfluidic velocimetry reveals spatial cooperativity in the flow of soft glassy materials[END_REF] contradicts this prediction. It thus shows that the plastic flow rate must depend on additional variables. The observations have been shown [START_REF] Goyon | Microfluidic velocimetry reveals spatial cooperativity in the flow of soft glassy materials[END_REF] to be compatible with the existence of diffusive, scalar quantity Γ, which they call fluidity.

It will be a challenge, in future studies, to identify the microscopic origin of such a quantity, which should be truly tensorial in more general situations than plane shear. In the quasistatic limit, fluidity should probably be related to the frozen stress disorder identified by Kabla and Debrégeas [START_REF] Kabla | Local stress relaxation and shear banding in a dry foam under shear[END_REF]. must include a term that reflects permeation. This effect can be incorporated into Eq. ( 56) in a rough manner by adding a diffusion term:

∂ρ ∂t + ∇ • (ρ v) = D permeation ∆ρ (94)
where the diffusion coefficient D permeation depends on the geometrical dimensions of the Plateau borders between foam vertices (which depend mainly on the volume fraction ϕ) but also on the hydrodynamic boundary conditions along the Plateau borders (which depend on various preparation conditions such as surfactant nature and concentration, salt, etc, in a very non-trivial manner).

E. Towards dilatancy

Eq. ( 94) above is only the first step towards a model that would include other effects known to exist in foams and emulsions. Indeed, relaxing the assumptions made in paragraph III F 1 would enable us to include an important phenomenon which is well known in the context of granular media and which has been recently demonstrated in liquid foams [START_REF] Marze | Protein and surfactant foams: linear rheology and dilatancy effects[END_REF],

namely dilatancy. When deformed, the local conformation of the foam, as schematized on Fig. 8, tends to go from conformation (0) to conformation (4). The physical origin of such a phenomenon is not clearly understood yet. It might be related to the thickening of films [START_REF] Emile | Swelling of a single foam film under slipping[END_REF] or Plateau borders [START_REF] Emile | Swelling of a single foam film under slipping[END_REF][START_REF] Terriac | Viscous force exerted on a foam at a solid boundary: influence of the liquid fraction and of the bubble size[END_REF] observed when a Plateau border glides at a solid wall.

Together with the identification of the microscopic origin of fluidity (see paragraph X C above), including the effect of dilatancy will thus be yet another challenge in the everbewildering rheology of foams and emulsions. where s (resp., c) denotes the sine (resp., the cosine) of angle θ between axes x and X (see Fig. 17).

Inserting the above equations into the evolution equation ( 27), we obtain:

˙ c 2 β 1 + s 2 β 2 = 2 γcs(β 1 -β 2 ) -2(c 2 δ 1 + s 2 δ 2 ) (C4) ˙ cs(β 1 -β 2 ) = γ(s 2 β 1 + c 2 β 2 ) -2cs(δ 1 -δ 2 ) (C5) ˙ s 2 β 1 + c 2 β 2 = -2(s 2 δ 1 + c 2 δ 2 ) (C6)
Summing Eqs. (C4) and (C6) yields the evolution of β 1 + β 2 . The difference between Eqs. (C4) and (C6), multiplied by (β 1 -β 2 ) cos(2θ), plus Eq. (C5) multiplied by 2(β 1 -β 2 ) sin(2θ), yields the evolution of β 1 -β 2 . From there, the evolution of β 1 and β 2 is readily obtained, see Eqs. (66) and (67). Using these equations, the difference between Eqs. (C4) and (C6) then yields the evolution of cos(2θ), see Eq. (68)

Elastic shear flow

Let the initial configuration be a material at rest, with β 1 = β 2 = 1. The system of Eqs. (66), (67) and (68) has a singularity at t = 0 since Eq. (68) contains a factor 1 β 1 -β 2 . In fact, the system of Eqs. (C4), (C5) and (C6) can be solved explicitely in the domain where

FIG. 5 :

 5 FIG. 5: Tentative variation of foam rheology as a function of its volume fraction and of the continuous phase viscosity. Four instances of the local aspect of the foam are drawn in the center of the figure, for low (light grey) or high (dark grey) continuous phase viscosity and for dispersed (round) or concentrated (faceted) bubbles. In each case, a tentative corresponding rheological model is schematically represented in terms of one spring, one solid friction element, one viscous friction element, and one general viscous element. The strength of each element is coded as weak (dashed line), medium (thin solid line) or strong (thick line).

FIG. 6 :

 6 FIG.6: With the Bingham model, from which we develop a fully tensorial model in the present work, the response to weak stresses is elastic (modulus G); by contrast, the response to large stresses presents both a yield threshold (σ y ) and a viscous component (η).

FIG. 7 :

 7 FIG. 7: Behaviour of the Bingham model (time t, deformation ε, deformation rate D, stress σ) in three series of experiments: (a) quasistatic imposed deformation (with a zoom on the stress fluctuations due to individual T 1 processes and corresponding elastic loading periods); (b) constant, imposed deformation rate (top: value of plateau stress); (c) constant, imposed stress.

FIG. 8 :

 8 FIG.8: Four changes in the local structure of a foam. (0) Initial structure. (1) Structure obtained through diffusion of gas from neighbouring regions into the bubble: the number of gas molecules in the bubble (represented by black dots) has increased, while the amount of liquid that surrounds the bubble has not changed. (2) Structure obtained by isotropic stretching: the quantity of gas and the quantity of liquid have not changed. This is not a stress-free conformation. (3) Structure obtained by deviatoric stretching: the volume, the quantity of gas and the quantity of liquid have not changed. This is not a stress-free conformation. (4) Structure obtained by permeation of liquid from the vicinity of the bubble towards other places in the foam: the amount of liquid that surrounds the bubble has decreased, while the number of gas molecules has not changed. In the present work, only changes where both phases are transported simultaneously, such as (0) → (2) and (0) → (3), are considered.

FIG. 9 :

 9 FIG.9: Evolution of stored deformation e and stress σ in the elastic regime. The evolution de/dt of the stored deformation e is purely kinematic, as it directly results from the convection of the material the velocity gradient ∇ v. As for the stress σ, it is related to the current stored deformation e through the elastic law. As a result, by composition, the evolution dσ/dt of the stress results from the expression of de/dt and from the elastic law.

FIG. 11 : 2 (

 112 FIG.11:Schematized T 1 process in two dimensions: conformations of four neighbouring objects (quadruplet) subjected to external forces. One pair of objects is oriented along axis ξ and the other one along axis ζ (see inset). The stress principal axes are assumed to be aligned with those of the quadruplet, which therefore adopts a diamond configuration. Moreover, the applied pressure is assumed to be constant. The diagramme indicates how both pairs compare as for the (centreto-centre) inter-object distance (horizontal axis) and as for the external stress that they undergo (vertical axis). By convention, positive values of the stress Σ ξξ (resp. Σ ζζ ) denote traction on the pair that is oriented along axis ξ (resp. ζ). A T 1 process starts from the relaxed situation R 1 , denoted by a black circle. Through compression along axis ξ (i.e., Σ ξξ < 0) or traction along ζ (Σ ζζ > 0) or through a combination of both, the system follows branch 1 (thicker line) towards the critical point C 1 . The T 1 itself is the sudden transition from C 1 (on branch 1) to branch 2 (thinner line). In the case of an incompressible material the quadruplet conformation depends only on Σ ξξ -Σ ζζ , while for a compressible material it additionally depends on Σ ξξ + Σ ζζ = -p.

FIG. 12 :

 12 FIG.12:From T 1 process to plastic flow. Left: the stress needed for four bubbles to undergo a T 1 process is illustrated on Fig.11. Centre: when a chunk of material contains a number of quadruplets of bubbles with various orientations, it may be postulated that the plasticity threshold corresponds to the most favourably oriented quadruplet (this assumption discards the effects of non-homogeneous deformations, and is discussed in the main text). Right: the rate at which T 1 processes occur in the chunk of material when it is subjected to elongation B is expressed in the form of a plastic flow D B p , where the material configuration at time t (dashed line) is mapped onto configuration at time t + dt (full line) through infinitesimal deformation D B p dt.

FFIG. 13 :

 13 FIG. 13: Incompressibility and plasticity. The incompressible evolution DB p of the relaxed conformation implies the usual property tr( DB p ) = 0. Despite this, the corresponding flow D B p in real space is not traceless in general (tr(D B p ) = 0): it is related to DB p through a finite transformation B, see Eq. (38) for details.

FIG. 14 :

 14 FIG. 14: Shear stress σ12 in the course of a shear experiment, as a function of the shear deformation γ, for three different shear rate values γ. The elasticity and plasticity terms are chosen as described

FIG. 15 :

 15 FIG.15: Dissipation P dissip per unit volume and injected power σ12 γ in the course of a shear experiment, as a function of the shear deformation γ. The applied shear rate is γ = 8. The elasticity and plasticity terms are chosen as described in Section VIII B. Note that there is no dissipation (P dissip = 0) in the elastic regime, prior to the onset of plasticity. Part of the injected work is not dissipated (non-zero surface area between both curves) and corresponds to the elastic energy stored in the material.

= 3a 1 FIG. 16 :

 116 FIG.16: Shear stress σ12 and dissipation power P dissip per unit volume in the stationary regime, as a function of the shear rate γ. The elasticity and plasticity terms are chosen as described in Section VIII B.

FIG. 17 :

 17 FIG.17:The homogeneous shear flow imposed in basis xyz attached to the laboratory (z being the direction of vorticity) causes the principal axes of tensor B to tilt (basis XY z) with respect to the shear. The lengths of the ellipse axes represent the magnitudes of eigenvalues β 1 and β 2 .

  x) = 1 when x ≥ 0 and θ(x) = 0 otherwise. Note that in Section IX, we explore the case b 2 = b 1 for comparison.

FIG. 18 :

 18 FIG. 18: Form of the stored elastic deformation in the course of an experiment and in the stationary regime, for three different values of the shear rate γ. The axes are the first two eigenvalues, β 1 and β 2 , of tensor B.

elongation β 1 (FIG. 20 :

 120 FIG. 19: Eigenvalues of the stored elastic deformation in the course of an experiment as a function of the shear deformation γ. The applied shear rate is γ = 8.

FIG. 21 :FIG. 22 :

 2122 FIG.21: Shear stress σ12 , deviatoric stress σ3 in the vorticity direction and normal stress differences N 1 and N 2 in the course of an experiment as a function of the shear deformation γ, for γ = 8.

FIG. 23 :FIG. 24 :

 2324 FIG. 23: Impact of elasticity on the stress response. The shear stress σ12 , the deviatoric stress σ3 in the vorticity direction and the normal stress differences N 1 and N 2 are plotted as a function of the shear deformation γ, for a shear rate γ = 8. Mooney-Rivlin elasticity is chosen with different parameter values. Left: k 1 = 0 and k 2 = 1. Center: k 1 = 1/2 and k 2 = 1/2. Right: k 1 = 1 and k 2 = 0.

FIG. 25 :

 25 FIG. 25: Burger model: schematic diagramme (left) and plot (right) of plastic flow as a function of tensor B, in the stationary regime and under small applied streses. In Burger's model, suitable for foams under low applied stresses (see paragraph III A), the plastic flow increases with the applied stress. If we neglect the short time scale response provided by the Voigt-Kelvin element (spring and viscous element in parallel), then this increase is weak and linear, as determined by the viscous element in series.

BFIG. 26 :

 26 FIG. 26: Bingham model: schematic diagramme (left) and plot (right) of plastic flow as a function of tensor B, in the stationary regime. In a model based on Bingham's model, such as that developed in the present work, the plastic flow is zero up to a certain stress (yield stress). Beyond the yield stress, it increases (linearly with the applied stress).

FIG. 27 :

 27 FIG.27:Combined Bingham-Burger model: schematic diagramme (left) and plot (right) of plastic flow as a function of tensor B, in the stationary regime. In a combined Bingham-Burger model, the plastic flow would increase weakly with the applied stress at low stresses, and more strongly above the yield stress. We believe it would mimic the rheological behaviour of a foam quite adequately.

  as mentioned in paragraph X B above, one can include the effect of, say, mechanical noise by including a dependence of D B p on the applied deformation rate D: D B p = D B p (B, D)

c 2 β 1 + s 2 β 2 cs(β 1 -β 2 ) 0 cs(β 1 -β 2 ) s 2 β 1 + c 2 β 2 c 2 δ 1 + s 2 δ 2 cs(δ 1 -δ 2 ) 0 cs(δ 1 -δ 2 ) s 2 δ 1 + c 2 δ 2

 2110121221101212

  1. the plastic flow depends on the stress or equivalently on the stored deformation (see Eqs. 35 or 37); 2. if applicable, the plastic flow is incompressible (see Eqs. 40 or 41); 3. the associated dissipation is positive (see Eq. 52). In this brief section, let us discuss how the constitutive equation derived in the present work can be inserted into a set of equations and provide a complete, continuous model. Keeping in mind the considerations of paragraph III F, let us now close the evolution equation of the system (Eqs. 27 or 28)

	VII. COMPLETE, CONTINUOUS MODEL

Acknowledgements

We gratefully thank Miguel Aubouy, Annie Colin, Franois Graner, Jean-Charles Razafindrakoto, Reinhard Höhler and Philippe Marmottant for stimulating discussions, and Pierre Rognon for a critical reading of the manuscript.

D. Liquid permeation channel dry region wet region FIG. 28: Foam with inhomogeneities in liquid volume fraction. The gas/liquid interface is flat in inter-bubble films, and curved in vertex regions at the junction between 3 (in two dimensions) or 4 (in three dimensions) bubbles. This curvature implies that the liquid pressure in the vertices is smaller than the gas pressure in the neighbouring bubbles. In wet regions, the interface curvature is less pronounced than in dry regions. As a result, there exists a pressure gradient, and the liquid tends to flow from wet regions towards dry regions. The intensity of the liquid flow depends on the hydraulic resistance in the Plateau borders that convey most of the liquid (depicted as the channel in this two-dimensional drawing).

Permeation of the continuous, liquid phase through the network of channels formed by the bubbles or droplets may result from gravity, due to the density mismatch between both phases, and is called drainage. Even in the absence of gravity, permeation may occur when the volume fraction of both phases is not uniform throughout the foam. Indeed, as depicted on Fig. 28, the pressure in the liquid phase is lower in dryer regions than in wetter regions for otherwise identical pressure values in the gas bubbles. As a result, the liquid permeates from wetter regions towards dryer regions.

Since vector field v represents the bubble (or droplet) velocity while the material density ρ includes the mass of the liquid phase (see paragraph III F), the density conservation Eq. [START_REF]To derive Eq. (17[END_REF] x 29: Augmented version of Figure 10. As mentioned in the caption of Figure 3, the relaxed state local orientation is chosen in such a way that tensor F is a pure deformation. As a consequence, points X and X ′ do not coincide except for particular values of ∇ v.

APPENDIX A: EVOLUTION OF FINGER TENSOR

In order to derive the evolution of tensor B (see Eq. 14) in some known velocity field, we compute how an ellipse with generic point x (see Eq. 13) is convected into a new ellipse (see Fig. 10). The generic point x ′ of the new ellipse is obtained as:

The equation for the new ellipse can be written in two ways:

The evolution of tensor B -1 = F -2 is thus given by:

From Eq. (A2) above, we then compute [START_REF]To derive Eq. (17[END_REF] the evolution of tensor B, see Eq. [START_REF] Bécu | Spatiotemporal dynamics of wormlike micelles under shear[END_REF].

APPENDIX B: DENSITY AND STORED DEFORMATION

Our assumptions stated in paragraph III F 1 imply that the material density is related to the determinant of tensor B in a very simply way. The evolution of det B can always be written in the form:

After multiplying Eq. ( 27) by B -1 , taking the trace and using the plastic incompressibility expressed by Eq. ( 40), we insert the result into Eq. (B1) and obtain:

and finally:

Comparing Eqs. ( 56) and (B3) shows that for any given element of material, the current density and deformation are linked via their initial values:

where r 0 is the position of material point r at time t 0 . The shear velocity gradient, the material deformation B and the flow rate D B p have the the material is purely elastic, with δ 1 = δ 2 = 0. Using the same notation as in Appendix C 1:

Whence:

Plasticity onset

When in the elastic regime, the system is described by Eqs.(...) above. Once the threshold is reached, given by Eq. ( 82), the plastic term comes into play.

Eq. (82) by itself describes the limit of the elastic regime in terms of tensor B for any choice of the system history in the elasic regime (not just the simple shear implemented here). It is represented on Fig. 19.

Stationary shear flow

When the flow is stationary, the system of Eqs (66-68) can be simplified. In particular, Eq. ( 68) then implies:

and finally the tilt θ of basis XY relative to basis xy can be expressed as:

From Eqs (66) and (67), still in the stationary regime, we obtain:

Combining this with Eq. ( 40) which expresses the material incompressibility, we obtain

Again from Eqs (66) and (67) in the stationary regime, we derive the shear rate:

Using Eq. (C15), we then obtain:

Thus, the stationary branch of the flow plots are obtained by choosing pairs of values for β 1 and β 2 that satisfy Eq. (C17), where δ 1 and δ 2 are functions of β 1 , β 2 and

To obtain the relevant values for β 1 and β 2 , we first find the pair that verifies simulataneously the threshold condition given by Eq. ( 82) and the stationary condition given by Eq. (C17).

We then use a differential equation derived from Eq. (C17) to follow the corresponding curve in the β 1 -β 2 plane.

Once the values for β 1 and β 2 are known, the shear rate is derived via Eq. (C19).

We can now slightly simplify Eqs. ( 84)-(91) in the stationary regime.

The first and second normal stress difference can be expressed as:

In the case of Mooney-Rivlin, coefficients a 1 and a 2 are given by Eq. ( 76) and the normal stress differences can be expressed as: