Shear jamming and fragility of suspensions in a continuum model with elastic constraints
SShear jamming and fragility of suspensions in a continuum model with elasticconstraints
Giulio G. Giusteri ∗ Dipartimento di Matematica, Università degli Studi di Padova, Via Trieste 63, 35121, Padova, Italy
Ryohei Seto † Wenzhou Institute, University of Chinese Academy of Sciences andGraduate School of Simulation Studies, University of Hyogo (Dated: February 4, 2021)Under an applied traction, highly concentrated suspensions of solid particles in fluids can turnfrom a state in which they flow to a state in which they counteract the traction as an elasticsolid: a shear-jammed state. Remarkably, the suspension can turn back to the flowing state simplyby inverting the traction. A tensorial model is presented and tested in paradigmatic cases. Itseffectiveness in reproducing the phenomenology of shear jamming shows that it is advantageous tolink this effect to the elastic response supported by the microstructure rather than to a divergenceof the viscosity.
Introduction. —The rheology of highly concentratedsuspensions of solid particles dispersed in a viscous fluidfeatures a number of surprising phenomena [1–3] amongwhich shear jamming raises important questions for itsinterpretation and challenges for its mathematical model-ing. If the concentration of particles is not very high, thesuspension presents a fluid-like behavior. At high con-centration, one can instead observe a sudden solidifica-tion that occurs after some strain in a shear deformation,whence the name of shear jamming.Under those conditions, the constant stress applied tothe suspension is balanced by the elastic response of thesolidified medium, arguably sustained by the networkof contacts developed among the solid particles duringthe initial flow [4–7]. The shear-jammed material is ina fragile state: if the applied stress is removed no mo-tion arises, but if we reverse the stress, pushing in asufficiently different direction, the suspension flows againand stops only once a certain strain is accumulated. Thishistory-dependent or protocol-dependent rheological re-sponse marks a key difference with states of isotropicjamming, in which the suspension is so concentrated tobe unable to flow, irrespective of the direction of the ap-plied forces.We present a constitutive model that is able to repro-duce the phenomenology of shear jamming and fragility.We base its construction on the understanding that shearjamming corresponds to the onset of an elastic response,related to geometric constraints at the micro-scale, ratherthan to a boost in dissipative phenomena, as implied bymodels containing a divergence of the viscosity. Thanksto this perspective, we devised an effective and yet math-ematically simple model, that features a small number ofparameters clearly associated with physical processes.We define a tensorial model at the outset, without go-ing through the process of designing a scalar model—typically tailored to a restricted set of motions—and thenextending it. In this way, our model is readily applicable to flows in general two- and three-dimensional geome-tries, both boundary- and pressure-driven ones.In the last decade, a few models to describe the physicsof dense suspensions have been proposed, with a focuson capturing discontinuous shear thickening and ratio-nalizing the role of frictional contacts between the solidparticles [8–14]. A common feature of such models is theattempt at building a direct link between the emergentbehavior and the microstructure of the fluid, as charac-terized by the analysis of data coming form detailed sim-ulations at the micro-scale level [15, 16]. In relation toshear jamming, these models are not readily applicableto the macro-scale simulation of this phenomenon, sincethey identify the jammed state with a divergence of theviscosity [17] rather than the onset of elastic responses,so that mathematical singularities appear in the equa-tions. On the other hand, we can find important resultsin the construction of continuum models that focus onthe macro-scale dynamics [18] and utilize a two-materialapproach, by following the coupled evolution of homog-enized fluid and solid phases. Such models are quite ef-fective in reproducing certain observations, but featureseveral parameters that require calibration and involve aset of equations with considerable complexity.Our model is intended to capture the suspension be-havior at a rather large scale, when a single-fluid model isappropriate. We do not take explicitly into account parti-cle migration phenomena [19–21] and rate-dependence ofthe viscosity or stress-induced solidification [17, 22, 23].These can be incorporated for specific applications as ex-tensions of the model but are not essential to reproduceshear jamming. In fact, we stress that it is possible togive a very good qualitative description of shear jammingassuming a constant viscosity, by associating jammingwith the appearance of elasticity.The phenomenon of shear jamming can be viewed asthe emergence of solidity due to the evolution of the sus-pension microstructure. The activation of frictional con- a r X i v : . [ c ond - m a t . s o f t ] F e b tacts between the particles leads to the presence of perco-lating stable formations that span macroscopic portionsof the system [5, 6, 24–26]. This microscopic non-localityof the internal interactions marks the transition from aregime in which momentum is transferred slowly and dif-fusively (viscous fluid) to a regime in which momentumtravels fast and elastically across the system (jammedsolid). When the elastic response is very stiff, one caneven approach the macroscopic non-locality representedby rigid-body motions.Another important aspect brought at the forefront byshear jamming is the material memory. We obviously ob-serve a long-lasting memory of what would be the relaxedconfiguration in the jammed solid regime, but there isalso a memory in the microstructure evolution that gov-erns the type and amount of deformation possible withinthe fluid regime, in which no persistent elasticity is de-tected [16, 27]. Both these aspects need to be capturedin an effective continuum theory and we propose to usetensorial models for all of them. As we shall see, thekinematic descriptors of the system that are useful forour purposes are the velocity field and its gradient, tocapture the viscous dissipation, and a tensorial measureof the strain induced on the material by the motion. Thelatter quantity keeps track of the microstructural defor-mation and, by limiting its evolution with unilateral con-straints in the appropriate space of tensors, we can cap-ture the transition between the fluid regime and the solidone, meanwhile preserving the characteristic reversibilityrepresented by the fragility of the shear-jammed state. The tensorial model. —A crucial role in shear jammingis played by the history of the deformation, since it in-duces some organization of the suspension microstruc-ture, eventually responsible for the solid-like behavior.Alongside the evolution equation for the velocity field u of a fluid with mass density ρ , ρ (cid:18) ∂ u ∂t + ( u · ∇ ) u (cid:19) = div T , (1)driven by the Cauchy stress tensor T , we consider theevolution equation for the deformation gradient tensor F [see 28, Chap. 3.2] in spatial coordinates : ∂ F ∂t + ( u · ∇ ) F = ( ∇ u ) F . (2)Equation (2) is an exact kinematic relation between thevelocity and the displacement of fluid elements and doesnot contain any constitutive assumption.From F we define B ≡ FF T and L ≡
12 log B (3)where log denotes the matrix logarithm. This is well-defined because the left Cauchy–Green tensor B is sym-metric and positive definite for any physical motion. These kinematic quantities track the local strain by fac-toring out rigid rotations, that should not affect the elas-tic response. The tensor L is the spatial counterpart ofthe Hencky strain and a generalization of the scalar strainmeasured in simple shear flows.Another important feature of L is that it is trace-less whenever det B = 1 . This is always the case forus, because we assume incompressibility of the material,namely div u = 0 at all times. We write the stress ten-sor as a pressure term plus the traceless extra stress S , sothat T = − p I + S . The extra stress is the sum of a viscousdissipation plus an elastic response. The dissipative termtakes the familiar form η D , wherein the effective viscos-ity η of the suspension multiplies the symmetric part ofthe velocity gradient D ≡ ( ∇ u + ∇ u T ) / .Regarding the elastic contribution to the stress, we as-sume that there exists a predetermined region N in thespace S of local strains (symmetric and traceless tensors)within which the material is elastically neutral . It meansthat, at each point x and instant t , if L ( x , t ) is in N thereis no elastic response. Moreover, the elastic response willbe proportional to a suitable measure of how far L ( x , t ) isfrom N . The overall isotropy of the suspension suggeststo take N to be a ball centered at the null tensor, namely N ≡ (cid:8) M ∈ S : (cid:107) M (cid:107) ≤ r (cid:9) where (cid:107) M (cid:107) ≡ (cid:113) tr( M T M ) / and r > , the radius ofthe region N , is a dimensionless material parameter thatindicates how much the suspension needs to be sheared toachieve a jammed microstructure. The value of r wouldtypically depend on the volume fraction of solid particles.The case r = 0 describes isotropic jamming, under whichthe suspension behaves always like an elastic solid.Since N is a closed convex subset of S , a projectionoperator Π :
S → N is well-defined and, for any M ∈ S ,the tensor Π( M ) is the element of N closest to M . Sucha projection can be easily expressed as Π( M ) ≡ (cid:40) M if (cid:107) M (cid:107) ≤ r,r M / (cid:107) M (cid:107) if (cid:107) M (cid:107) > r. (4)To reflect the fact that an elastic response is activatedwhenever the logarithmic measure of strain L leaves theneutral region N , we assume an extra stress of the form S = 2 η D + 2 κ (cid:0) L − Π( L ) (cid:1) , (5)where the material parameter κ > represents an elasticstiffness. This is a “soft” way of constraining the strain(as opposed to keeping it always within N ) that is ableto better reproduce some details of the elastic effects ob-served in the proximity of jamming [7], while being easierto handle in simulations. With the present model, thestrain of the jammed material tends to remain close tothe boundary of N if the applied stress is driving it out-wards. Conversely, the suspension can flow again as soonas the stress drives L towards the interior of N . In thisway we can capture both shear jamming and the fragilityof the jammed state.The inclusion of additional dissipative phenomena,that may appear at the onset of jamming, can be achievedby letting η depend on a parameter like λ ≡ (cid:107) L − Π( L ) (cid:107) .For instance, one could assume η = η + η j Θ( λ ) , where η is the viscosity of the suspension in the fluid regimeand η j ≥ is an additional dissipation observed in thejammed regime. Θ( x ) is the Heaviside function, that is for x > and 0 otherwise.The logarithmic measure of strain L is convenient,especially in the presence of the incompressibility con-straint. Indeed, being L traceless, a stress proportionalto it acts tangentially to that constraint. Moreover, itdrives the system along geodesics in the group of trans-formations with unit determinant, that is a Lie groupwith nonzero curvature. Planar extensional flows. —We highlight the basic fea-tures of the model in an idealized case, for which analyti-cal computations can be carried out. Under the deforma-tion associated with planar extensional flows, the currentposition of a particle that occupies the place ( x , y ) attime is given by ϕ ( x , y , t ) = ( x e ε ( t ) , y e − ε ( t ) ) and itsspatial inverse is ϕ − ( x, y, t ) = ( xe − ε ( t ) , ye ε ( t ) ) , where ε ( t ) is an arbitrary function of time and measures thestrain of the material. We immediately obtain F = (cid:18) e ε ( t ) e − ε ( t ) (cid:19) and B = (cid:18) e ε ( t ) e − ε ( t ) (cid:19) . In this case, the computation of the matrix logarithm isstraightforward and yields L = (cid:18) ε ( t ) 00 − ε ( t ) (cid:19) . The velocity is u ( t ) = ( ˙ ε ( t ) x, − ˙ ε ( t ) y ) and the symmet-ric part of the velocity gradient, the usual measure of therate of deformation, is D = (cid:18) ˙ ε ( t ) 00 − ˙ ε ( t ) (cid:19) = ∂ L ∂t . (6)We stress that the second identity in (6) is not valid for ageneric flow (it does not hold, e.g., in simple shear); whenrotation and vorticity are present, they affect the defor-mation history in a nontrivial way, and D and L cannotremain aligned. This fact corresponds to the well knownpresence of normal stress differences in simple shear flowsof viscoelastic fluids.We consider the extensional flow in a cross channelwith hyperbolic boundaries that allow for a perfect slipof the fluid (Fig. 1, left). A pressure difference applied toinlets and outlets of the channel at distance (cid:96) from thecenter generates normal tractions ± τ . In a slow-velocityregime, the linearized flow equations give the pressure Figure 1. We can imagine planar extension in a cross channel(left), with hyperbolic boundaries (red solid lines) that allowfor perfect slip, to which we apply a pressure difference τ between inlets and outlets (blue dashed lines). The linearizedequations for the proposed model reduce in this case to ascalar ordinary differential equation for the strain function ε .It corresponds to that of a damped oscillator with a potentialenergy V featuring a flat region for ε ∈ [ − r, r ] (right). field p ( x, y, t ) = ρ ¨ ε ( t )( y − x ) / . The balance of stressat the inlets or outlets yields the following equation: ρ(cid:96) ε + 2 η ˙ ε − τ = − κ [ ε + r ] if ε < − r, +0 if − r ≤ ε ≤ r, − κ [ ε − r ] if ε > r. (7)This is a scalar ordinary differential equation for thestrain ε ( t ) , equivalent to that of a damped oscillator withpotential energy that features a flat region for ε ∈ [ − r, r ] and parabolic branches outside that interval (Fig. 1,right). This entails transitions between a viscous fluidbehavior, for ε ∈ [ − r, r ] , and that of a viscoelastic solid. Clogging and unclogging. —Let us now consider how themodel performs in simulating a paradigmatic pressure-driven flow through a contraction. In this planar flow,the width of the contraction is one fourth of the chan-nel width (that equals its depth) and of the contractionlength. Periodic boundary conditions for u and F areassumed at the ends of the channel, while a pressure dif-ference ∆ p between the two openings is driving the flow.We introduce a dimensionless form of the evolutionequations by defining a reference pressure difference P and a reference length L given by the width at the in-let. A reference time scale is then t ≡ L (cid:112) ρ/P . Fromthese, we set the Reynolds number Re and the Weis-senberg number Wi (ratio of elastic to viscous forces)according to Re ≡ L √ ρPη and Wi ≡ κLη (cid:114) ρP , and so that the dimensionless stress reads ˜ T = − ˜ p I + 2Re ˜ D + 2WiRe (cid:0) L − Π( L ) (cid:1) . (8)In what follows we consider all quantities as dimension-less but drop the tildes for simplicity. Figure 2. The present model can reproduce many qualita-tive features of shear jamming in a complex flow througha contraction. We varied the (dimensionless) pressure dif-ference ∆ p ( t ) imposed between the left and the right open-ing and measured the flow rate q ( t ) through a cross sectionof the channel and its time integral Q ( t ) = (cid:82) t q ( s ) ds . Pa-rameters in the simulation: channel width and contractionlength is , contraction width is . , Re = 20 √ . ≈ , Wi = 100 / √ . ≈ , r = 1 . / √ , and η j = 0 as we keptalways a fixed viscosity. At startup, ∆ p is set positive, the flow accelerates andthe flow rate reaches a maximum at about t = 1 (Fig. 2).The deformation induced by the flow gives rise to shear-jammed domains that grow from the boundaries towardsthe center of the contraction. This activates an elasticresponse within the material (Fig. 3b) that hinders theflow. When we remove the pressure difference, from t =8 to t = 10 (and again from t = 14 to t = 16 ), thestored elastic energy is completely released with a smallrecoil and then the flow stops (Fig. 2). Nevertheless, themicrostructure remembers to be close to shear jammingand when the pressure difference is turned on again ( t =10 ) only a small fluid displacement is produced, sincewe assist at a rapid reactivation of the elastic responsedistributed along the contraction (Fig. 3c).On the other hand, when ∆ p is reversed to a nega-tive value, the flow lasts for a longer time and the fluiddisplaced through the contraction before shear jammingsets in again is about twice as much as that displacedin the first part of the experiment (Fig. 2). The shear-jammed domains, where the elastic response is active, aredestroyed and rebuilt with a different spatial distributionby the reverse flow (Fig. 3d). The elastic stress at t = 25 is concentrated in a smaller region and more intense thanthat at t = 14 , since it balances an equivalent pressuredrop. By mimicking randomness in the suspension mi-crostructure [22] with spatial fluctuations of the initialdeformation gradient, we obtained a more realistic nu-cleation of the shear-jammed regions. Movies with the Figure 3. (a) The entire system is shown with the dis-cretized mesh. (b–d) The clogging of the channel is due tothe presence, within the contraction to which images are re-stricted, of shear-jammed domains. These are characterizedby a non-vanishing elastic response measured by the param-eter λ ≡ (cid:107) L − Π( L ) (cid:107) . There is a clear difference between theshear-jammed state achieved at t = 8 or t = 14 with a positive ∆ p , that pushes rightward, and the one obtained at t = 25 after reorganization with a negative ∆ p , that pushes leftward. full time-evolution of pressure, flow, and elastic responseare available [29]. Conclusions. —We have shown how the knowledgeabout a complex collective phenomenon, acquired bymeans of experiments and simulations, can be transferredinto a rather simple model of the macroscale physics thatwe observe. By relating shear jamming to the activationof an elastic response and not to a divergence of the vis-cosity, we developed a tensorial model able to reproducethe qualitative features of shear jamming.Such a model can be applied to generic flows and ge-ometries in both two and three dimensions, because itrests on physical considerations that are not peculiar toa specific experimental setup. It becomes particularlyuseful to simulate the flow of suspensions in applications,where the focus is on the emergent collective physics andnot on its microscopic origins.Notably, we simulated a material able to switch, in areversible way, between a fluid-like and a solid-like be-havior. This feature, essential to capture shear jamming,can suggest effective ways to deal also with yielding phe-nomena. While we kept the model as simple as possi-ble, many extensions can be implemented to reproduce arate-dependent behavior.G.G.G. acknowledges the support of the Italian Na-tional Group of Mathematical Physics (GNFM-INdAM)through the funding scheme