Continuum modeling of soft glassy materials under shear

Soft Glassy Materials (SGM) consist in dense amorphous assemblies of colloidal particles of multiple shapes, elasticity, and interactions, which confer upon them solid-like properties at rest. They are ubiquitously encountered in modern engineering, including additive manufacturing, semi-solid flow cells, dip coating, adhesive locomotion, where they are subjected to complex mechanical histories. Such processes often include a solid-to-liquid transition induced by large enough shear, which results in complex transient phenomena such as non-monotonic stress responses, i.e., stress overshoot, and spatially heterogeneous flows, e.g., shear banding or brittle failure. In the present article, we propose a pedagogical introduction to a continuum model based on a spatially resolved fluidity approach that we recently introduced to rationalize shear-induced yielding in SGMs. Our model, which relies upon non-local effects, quantitatively captures salient features associated with such complex flows, including the rate dependence of the stress overshoot, as well as transient shear-banded flows together with non-trivial scaling laws for fluidization times. This approach offers a versatile framework to account for subtle effects, such as avalanche-like phenomena, or the impact of boundary conditions, which we illustrate by including in our model the elasto-hydrodynamic slippage of soft particles compressed against solid surfaces.


I. INTRODUCTION
Soft Glassy Materials (SGMs) encompass a broad variety of colloidal particles densely packed into an amorphous microstructure showing solid-like properties at rest [1,2].These particles, which can be either soft and deformable or hard, form a jammed assembly with glassy-like mechanical properties characterized by (i) a linear viscoelastic response where the elastic contribution is dominant [3], and (ii) time-dependent properties referred to as "aging" in the literature [4,5].Moreover, under a sufficiently large external stress or strain, particles can rearrange.For vanishingly low shear rates, such rearrangements take the form of local plastic events such as T1 events in foams and emulsions [6], or sheartransformation zones in colloids [7].Eventually, for sufficiently large accumulated deformation, these plastic events, which act as a mechanical noise, lead to the fluidization of the sample.Remarkably, such a shearinduced solid-to-liquid transition displays generic features that are quite insensitive to the sample microstructure [8][9][10][11].For instance, under a constant applied shear rate γ, the stress σ builds up and reaches a maximum before relaxing towards a steady-state value.Such a nonmonotonic response, known as a stress overshoot [12], coincides with the yielding of the sample, which may either flow homogeneously, or rather display a spatially heterogeneous yielding process [13].In the latter case, flow heterogeneity occurs due to localized, brittle-like failure [14], or results from a more ductile process in which an ar-rested region coexists with a fluidized one, referred to as a "shear band," whose lifespan depends on the volume fraction and on the particle interactions [15][16][17].
Various modelling efforts have been undertaken over a broad range of spatial scales, from that of the building block, thanks to, e.g., Molecular Dynamic simulations [18][19][20], to mesoscale or macroscopic continuum approaches in which the SGM microstructure is accounted for only by a few parameters [12,21], up to typically 10, in order to capture more subtle effects such as nonisotropic resistance of the sample inherited from shear history [22,23].Here, we shall focus on a continuum approach traditionally referred to as "fluidity models" [24,25] in which the microscopic properties of the sample are expressed by the fluidity f , a local quantity, which stands for a rate of plastic events.
Recently, fluidity models were derived theoretically by Bocquet et al. [26] as the continuum limit of a microscopic equation for the probability distribution originally proposed by Hébraud and Lequeux [27].Such an approach showed that fluidity models naturally encompass non-local effects in steady state via a so-called "cooperativity" length scale that quantifies the extension of the region that is impacted by a neighboring plastic rearrangement [28][29][30][31].The approach of Ref. [26] was extended to transient flows by some of us in Ref. [32].This extended non-local version of the fluidity model has been used to quantitatively capture the key features of the yielding transition of a soft glass [33][34][35].
In the present Perspective, we first summarize these findings to illustrate the power of the non-local fluidity model.Second, we extend our approach and include the elasto-hydrodynamic slippage of soft particles compressed against a solid surface in order to describe recent results from the literature [36].

II. PHENOMENOLOGY OF SHEAR START-UP EXPERIMENTS
In a shear start-up experiment, one imposes a constant shear rate γ at time t = 0 upon an SGM initially at rest. Figure 1 illustrates the general phenomenology of shear start-up through a selection of experimental results on Carbopol microgels [37][38][39].As recalled in the introduction, a stress overshoot is classically observed: at short times t, the shear stress σ grows linearly with the strain γ = γt [see red dashed line in Fig. 1(a)], which is typical of an elastic response.At longer times, the stress progressively deviates from a linear response and reaches a maximum σ M at time t M .As the stress maximum is reached, the material has become strongly anisotropic, and subsequent stress relaxation processes lead the material to flow on even longer time scales.This global behaviour is typical of ductile-like yielding.
To get more insight into the local structure of the flow during the solid-to-liquid transition, velocity profiles measured using ultrasonic velocimetry are displayed in Fig. 1(b) and 1(c).In the case of the present microgels, the material is homogeneously strained prior to the stress overshoot, and the stress maximum corresponds to the point when the microgel fails at the shearing surface [see Fig. 1(b)].Such failure is followed by an elastic recoil [see the negative velocities for the velocity profile with • symbols in Fig. 1(b)], then by a fully arrested regime with v = 0 across the whole sample except for a thin, unresolved lubrication layer at the moving wall.Yet, on time scales much longer than the time t M of the stress maximum, a fluidized region, i.e., a shear band, of width b grows from the moving wall and coexists with the arrested, solid-like material until a fully homogeneous flowing state is reached at a well-defined fluidization time τ f [see Fig. 1(c)].Finally, whatever the complexity of the stress relaxation, the stationary velocity profiles of the present microgels are ultimately homogeneous with insignificant wall slip.It is also essential to note that the fully flowing material is well described by the widespread Herschel-Bulkley (HB) constitutive law relating the stress σ and the shear rate γ in the system at steady state [17], σ = σ y + A γn , where σ y is the yield stress, A the consistency and n the shear-thinning index.
From the experimental results displayed in Fig. 1, one may argue that the yielding transition can be considered as a dynamical first-order phase transition, where one phase (the fluid-like phase) nucleates into the other phase (the solid-like phase).As we shall see, this idea underlies most of the following discussion, which focuses on two relatively simple yet fundamental questions: (i) How does the stress overshoot σ M depend on the applied shear rate γ? (ii) How does the fluidization time τ f depend on the applied shear rate γ or stress σ?Our ultimate goal is to obtain a general framework that describes quantitatively the yielding transition and its mesoscopic features.

III. CONTINUUM MODELING
We start by considering a two-dimensional shear geometry where the SGM is confined between two infinite parallel plates separated by a distance L. The flow is assumed to be one-dimensional along the direction x, i.e., it is described by a velocity field v = (v x , v y ) with v x (x, y, t) = v(y, t) and v y (x, y, t) = 0, where y denotes the velocity gradient direction and t the time.The wall at y = 0 moves with a constant velocity v 0 imposed at the initial time t = 0, while the wall at y = L remains fixed with zero velocity.As in experiments, the SGM is assumed to be initially at rest so that v(y, 0) = 0 for y ∈ [0, L].For the sake of simplicity, we introduce the dimensionless stress and shear rate, Σ = σ/σ y and Γ = γ/(σ y /A) 1/n , such that the SGM in steady state follows the dimensionless HB law: Following Refs.[26,28], in order to describe the local behaviour of the SGM, we introduce the fluidity f (y, t) of the SGM as the relevant order parameter in the system, as well as a characteristic length scale, called the cooperativity scale ξ, which controls spatial dynamics of the fluidity.Qualitatively, the fluidity corresponds to the rate of plastic events at a given time and position in the system.When the SGM flows in steady-state under an applied stress Σ, Bocquet et al. [26] linked the fluidity to elasto-plasticity at mesoscale through the following equation: where and Θ is the Heaviside function.As noted in Ref. [26], Eq. ( 2) can be associated to the functional derivative of (4) From the above equation, it is tempting to consider F [f ] as a free energy functional for the fluidity f .Based on this idea, we proposed in Refs.[32,33] to extend the approach and formulate the dynamics of the system using F [f ].In order to model shear start-up, i.e., a constant velocity v 0 imposed at the moving wall at t = 0, we take the shear rate Γ as the imposed control parameter.In this case, the quantity f ≡ f / Γ should be proportional to the number of plastic events occurring at some position y over the time scale Γ−1 .Such a number may increase or decrease locally depending on the dynamics of the system induced by the external driving Γ.Since Γ is constant, the temporal variation ∂ f /∂t is nothing but the fluidity variation due to t ≡ Γt, i.e., ∂ f /∂t = ∂f /∂ t.Our first modeling step is to assume that the fluidity dynamics is given by: where κ[f ] plays the role of a "mobility".Suppose now that the system can be decomposed into two different regions: a fluidized region where f > 0 and a solid-like region where f = 0. Our second important assumption is to require that both regions correspond to stationary solutions of Eq. ( 5).Assuming κ[f ] to be an analytic function of f , the simplest choice is κ[f ] ∼ f .This implies that the formation of a shear band in the system coexisting with a solid-like region with exactly f = 0 can be described by the superposition of two stationary states of the dynamics: one corresponding to f > 0 [see Eq. ( 2)] and the other one to f = 0.Moreover, it can be easily understood that, if the solid-like region is described by a small yet non-vanishing fluidity f > 0, then it cannot remain solid forever and it will eventually flow, i.e., it is unstable.This situation therefore corresponds to transient shear banding as we shall describe in the following.The third and last modeling step is to couple Eq. ( 5) with an equation for the time evolution of the stress Σ(t), which we suppose spatially homogeneous.This can be done by decomposing the total strain Γ = Γ el + Γ pl into an elastic contribution Γ el = τ Σ, where τ is a characteristic time inversely proportional to the elastic modulus, and a plastic contribution Γ pl such that Γpl = Σ f , where . . .denotes spatial average, and by using the well-known Maxwell model [40].The evolution equation for the shear stress then reads: with Γ = Γel + Γpl .This last equation is coupled to Eq. ( 5), which can be rewritten as: with m given by Eq. ( 3).Finally, we must specify boundary and initial conditions.Assuming that the external driving is acting at the boundary y = 0, we choose the following boundary conditions.For m 2 > 0, we impose a "wall fluidity" f w = f (0, t) = m 2 (Σ) at the moving wall, and ∂ y f (L, t) = 0 at the fixed wall.When m = 0, we assume ∂ y f (0, t) = 0 = ∂ y f (L, t) at both walls.Moreover, we take the initial fluidity profile to be homogeneous and very small, i.e., f (y, 0) = f 0 with f 0 1.As soon as Σ > 1, we expect a shear band to develop from y = 0 with a size b (t) that increases with time and whose dynamics is set by the spatio-temporal evolution of the fluidity.
Figure 2 provides examples of numerical resolutions that illustrate the general phenomenology of shear startup in our fluidity model.In particular, the time evolution of Σ(t), obtained from the numerical integration of Eqs. ( 6) and ( 7) for ξ = 0.04, f 0 = 10 −4 , and Γ = 2, shows a stress overshoot very similar to that observed in experiments [Fig.2(a)].Moreover, the fluidity profiles f (y, t) displayed in Fig. 2(b) for two specific times t present a sharp interface between a fluidized shear band for y < b , where f ∼ m 2 , and a solid-like region for y > b , where f (y, t) = f 0 .Note that the shear band grows in size because of the instability of the solid-like region, while retaining a sharp interface at f f 0 .This results from our requirement that κ[f ] = f .More precisely, it is possible to estimate analytically the steepness of the interface at f f 0 , shown in dash-dotted lines in Fig. 2(b), as ∂ y f | y= b /m 2 = (m/5ξ 2 ) 1/2 .Deep into the region f = f 0 1, i.e., far enough from the interface, the fluidity increases algebraically in time as f (y, t) f 0 1 + f 0 t 0 m(s)ds .However, at the boundary of the fluid-like region for y b (t), an instability occurs with an exponential growth of the fluidity.Dimensional considerations suggest that the instability extends over a scale of order ξ/m 1/2 and grows with a characteristic time scale m −3 .This implies that the size of the shear band b ( t) satisfies the equation [35]: We will come back to the dynamics of the transient shear band below when discussing the fluidization.

IV. STRESS OVERSHOOT
We now investigate how the stress maximum Σ M scales with Γ, which is an observable classically extracted from experiments [38,41,42].Concomitantly to the stress overshoot, the size b of the shear band increases with time.As detailed in Ref. [35], at short time t, the band dynamics is dominated by the diffusion term f w ξ 2 ∼ m 2 ξ 2 and b grows as (m 2 ξ 2 t) 1/2 , while for large enough b , it follows Eq. ( 8).The overall process is illustrated in Fig. 2(c) by plotting the distance to the yield stress, Σ − 1, as a function of the effective shear rate L Γ/ b , which corresponds to the average shear rate in the fluidized band.Such a representation of the flow dynamics clearly highlights the separation between two different dynamical regimes: a short-time "unsteady" regime that strongly depends on the applied shear rate Γ and where the data fall well below the equilibrium HB curve, and another "quasi-steady" regime at longer times, where all data nicely collapse on the HB curve, including during the transient shear-banding regime.
Based on the previous observations, the scaling of Σ M with Γ can be computed using Eqs.( 6) and (7).In a nutshell (see Refs. [34,35] for full details), considering that the stress grows linearly up to the stress maximum, i.e., Σ( t) t/τ , that the stress is large enough that m = (Σ − 1) 1/2n /Σ 1/2 Σ 1/2n−1/2 , and that f ∼ b m 2 , the condition dΣ/d t = 0 at the stress maximum leads to: where tM is the strain at the stress maximum Σ M .Further analysis of the two dynamical regimes then yields: where B and C are two numerical prefactors, and The first term on the r.h.s of Eq. ( 10) dominates for small Γ, when the shear band grows due to the diffusion term f w ξ 2 , while the second term dominates for large Γ, when the shear band increases according to Eq. ( 8).Note that Eqs. ( 10) and ( 11) hold both for transient and stable shear bands.As shown in Ref. [35], an extensive survey of the existing numerical and experimental data shows excellent agreement with Eq. (10). Figure 3 illustrates this agreement by comparing the model predictions to experiments on Carbopol microgels.In both cases, two power-law regimes can be identified in Σ M − 1 vs Γ, and the exponents are consistent with the values β = 1/3 in the "diffusive" regime and with α = 4/17 in the "asymptotic" regime at large Γ predicted for a shear-thinning index n = 1/2.Finally, we emphasize that Eq. ( 10) depends on ξ with a singular limit for ξ → 0, and that the above results depend on the choice κ[f ] ∼ f .The good agreement between Eq. ( 10) and experimental data therefore provides strong support for the present formulation of the fluidity model, which constitutes a remarkable, non-trivial result.

V. INCLUDING ELASTO-HYDRODYNAMIC (EHD) INTERACTIONS INTO THE MODEL
Based on experiments and numerical simulations, Cloitre, Bonnecaze and collaborators [43][44][45] have shown that the flow of SGMs constituted of dense assemblies of deformable particles, such as microgels, emulsions or glasses of elastomeric particles, is controlled by elastohydrodynamic (EHD) interactions, which result from the lubrication flows of solvent within the thin films between the particles.In particular, a recent study [36] has shown that EHD interactions impact the scaling of the stress overshoot in a non-trivial way.We herewith discuss an easy way to include such EHD effects in our continuum model through a simple modification of the Maxwell equation ( 6) for the stress evolution.We propose to add a contribution Γ EHD from EHD interactions to the total strain, Γ = Γ el + Γ pl + Γ EHD , which is related to the shear stress through ΓEHD = Γ0 Σ 2 , where Γ0 is a reference shear rate below which EHD effects become significant.This specific choice of scaling for the EHD interactions is justified by [44,46].The resulting modified Maxwell model reads: First, EHD interactions modify the steady-state rheology.Indeed, with f = m 2 Θ(Σ − 1) and dΣ/d t = 0, we get: The inset in Fig. 4 compares the steady-state flow curve predicted from Eq. ( 13) with n = 1/2 and Γ0 = 0.04 to experimental data on microgels obtained under two different boundary conditions [46].While the experimental flow curve for rough shearing surfaces (brown squares) nicely follows the HB law, the flow curve measured for smooth surfaces (yellow circles) presents a kink for Γ 0.05 that is usually interpreted as the hallmark of predominant slippage at the walls [10].Interestingly, including EHD interactions in our model using Γ0 = 0.04 allows us to nicely predict the steady-state flow curve for the smooth surface: in spite of some deviations at extremely low shear rates, EHD contributions produce deviations from the HB behaviour when Γ Γ0 , while leaving the HB flow curve essentially unaltered for Σ 1, very much like experimental observations.This suggests that Γ0 most probably embeds some non-trivial dependence on boundary roughness, which remains to be modeled theoretically.
Therefore, for Γ < Γ0 where EHD effects dominate, the stress maximum no longer depends on the HB exponent, but rather simply on the EHD scaling as Σ M ∼ ( Γ/ Γ0 ) 1/2 .For Γ Γ0 , however, the scaling of Eqs. ( 10) and (11) is recovered.This is confirmed in Fig. 4 by the numerical integration of the full dynamical equations.Note that when EHD interactions dominate, the exponent 1/2 is observed for the stress maximum Σ M rather than for Σ M − 1, which indicates that the yield stress is no longer a "reference" stress for the stress overshoot.

VI. TRANSIENT SHEAR BANDING AND FLUIDIZATION TIME
The fluidity model can be further used to compute the duration of the transient shear banding regime, i.e., the fluidization time T f as a function of Γ, which constitutes an important prediction for experiments and applications of SGMs.Indeed, during the fluidization process, the system satisfies the balance Γ = f Σ = m 2 b Σ.This allows us to compute m as a function of b and Γ.Using Eq. ( 8), one then predicts that, for small enough Γ: in excellent agreement with experimental data [33,37].Note that the scaling exponent for T f vs Γ is independent of the HB exponent n.
In the case of stress-induced fluidization, i.e., when forcing at a constant external stress Σ, Eq. ( 8) can still be used upon identifying t = m 2 t, which results from the fact that m 2 ∼ Γ for small Γ.Thus, Eq. ( 8) generally predicts T f ∼ (ξm 9/2 ) −1 , which leads to T f ∼ 1/[ξ(Σ−1) 9/4n ] for small imposed values of Σ − 1.Therefore, the present fluidity model predicts that the ratio of the scaling exponents under imposed Γ to that under imposed Σ is given by the HB exponent n, as observed in experiments on Carbopol microgels [17,33].
Finally, as examined in details in Ref. [35], one may introduce long-range correlations in the fluidity through noise-like dynamics and investigate how the above predictions depend on boundary conditions.In brief, the transient shear-banding scenario and the scaling of T f given by Eq. ( 15) are very robust to fluidity correlations when the fluidity at the moving wall is fixed through f w = f (0, t) = m 2 .However, when rather fixing the fluidity gradient at the moving wall through ∂ y f (0, 0) = 0, long-range spatial correlations conspire with the boundary condition to promote the emergence of a completely different fluidization scenario, where the growth of the shear band is prevented, leading to a stress increase that is initially smoother, but later characterized by an abrupt drop, resembling brittle-like failure [14] and similar to the one discussed in recent theoretical and numerical works [18,47,48].

VII. SUMMARY AND OPEN QUESTIONS
We started this Perspective paper by asking how two classical observables that characterize shear start-up in SGMs, namely the stress overshoot Σ M = σ M /σ y and the fluidization time, depend upon the applied shear rate Γ.We have shown that a dynamical fluidity model allows one to predict the way rheological variables should be analyzed, i.e., Σ M − 1 vs. Γ and T f vs. Γ or m.The corresponding scaling exponents and their dependence on the HB exponent are in excellent agreement with experiments on Carbopol microgels.The model is versatile enough to include EHD interactions that, when dominant, change the scaling of the stress overshoot to Σ M ∼ Γ0 whatever the underlying HB behaviour.Overall, the model predictions hinge on some basic ingredients: (i) the "mobility" function κ[f ] in Eq. ( 5) that allows the coexistence of a fluidized band and a solidlike region, (ii) the Maxwell equation ( 6) for the stress evolution, and (iii) the boundary conditions, which are crucial, for they discriminate between ductile-like and brittle-like types of fluidization.Focusing on items (ii) and (iii) above, we highlight two important open problems.
First, future work should analyze situations where the forcing has some increased complexity.In particular, it would be very interesting to explore the present fluidity model with time-dependent protocols such as the shearrate ramps that are widely used by rheologists.Whether or not this model may predict rheological hysteresis in SGMs and its dependence with the shear-rate sweep rate [49][50][51], in the two cases of transient and permanent shear banding, is an outstanding task.
Second, accounting precisely for boundary conditions is key for further theoretical advances.As already noted in Ref. [35], since a simple change of boundary conditions may suppress the nucleation of the fluid-like phase at the moving wall, boundary conditions appear to control the shear-induced solid-to-liquid transition in SGMs.However, we still miss physical insight into the microscopic dynamical processes at play at the walls.Here, the pro-posed phenomenological treatment of EHD interactions and the observation that EHD parameters must depend on boundary conditions call for more modelling effort.This can open the way to obtain a realistic fit of experimental results once the various parameters for the continuum modelling are extracted from experiments.

FIG. 1 .
FIG. 1. Phenomenology of shear start-up experiments in Carbopol microgels.(a) Stress σ as a function of strain γ = γt recorded after a shear rate γ = 5, 1, 0.2, and 0.03 s −1 from top (darker color) to bottom (lighter color) is applied at time t = 0.The red dashed line highlights the linear response at short time σ = G0γ with G0 = 300 Pa.Inset: same data plotted as a function of time t using semilogarithmic scales.(b,c) Velocity profiles v normalized by the velocity of the moving plate v0 as a function of the distance y to the moving wall normalized by the gap size L and recorded (b) at short times around the stress overshoot under γ = 0.1 s −1 , and (c) at long times during the transient shear-banding regime under γ = 0.7 s −1 .In both cases, the inset shows the corresponding stress response σ(t) and the colored symbols show the times at which the velocity profiles in the main graph are recorded.The colored lines are guides to the eye in (b) and fits to the velocity profile in the shear band in (c).The gray dashed line in (c) shows the velocity profile expected for a Newtonian fluid in the absence of wall slip.