Beyond linear elasticity: Jammed solids at finite shear strain and rate

1.2 MB·

by  Julia Boschan

Beyond linear elasticity: Jammed solids at finite shear strain and rate
Julia Boschan,1 Daniel V˚agberg,1 Ell´ak Somfai,2 and Brian P. Tighe1
1Delft University of Technology, Process & Energy Laboratory,
Leeghwaterstraat 39, 2628 CB Delft, The Netherlands;
2Institute for Solid State Physics and Optics, Wigner Research Center for Physics,
Hungarian Academy of Sciences, P.O. Box 49, H-1525 Budapest, Hungary
(Dated: January 5, 2016)
Theshearresponseofsoftsolidscanbemodeledwithlinearelasticity,providedtheforcingisslow
andweak. Bothoftheseapproximationsmustbreakdownwhenthemateriallosesrigidity,suchas
in foams and emulsions at their (un)jamming point – suggesting that the window of linear elastic
response near jamming is exceedingly narrow. Yet precisely when and how this breakdown occurs
remains unclear. To answer these questions, we perform computer simulations of stress relaxation
6
and shear startup experiments in athermal soft sphere packings, the canonical model for jamming.
1
By systematically varying the strain amplitude, strain rate, distance to jamming, and system size,
0
we identify characteristic strain and time scales that quantify how and when the window of linear
2
elasticitycloses,andrelatethesescalestochangesinthemicroscopiccontactnetwork. Ourfindings
n indicatethatthemechanicalresponseofjammedsolidsaregenericallynonlinearandrate-dependent
a
on experimentally accessible strain and time scales.
J

1
Linear elasticity predicts that when an isotropic solid to nonlinearity and irreversibility in the particle trajec-

]
is sheared, the resulting stress σ is directly proportional tories, andeventuallytosteadyplasticflow[19–24]. Ina
t
f to the strain γ and independent of the strain rate γ˙, series of influential studies, Schreck and co-workers [25–
o
s 29] recently asked how many contact changes a jammed
. σ =G0γ, (1) packing undergoes before its mechanical response ceases
t
a to be linear. To answer this question, they studied the
m with a constant shear modulus G0 [1]. The constitu-
onsetofmixingbetweenexcitedvibrationalmodesinmi-
tive relation (1) – a special case of Hooke’s law – is a
- crocanonical ensembles of N particles, and found that
d simple, powerful, and widely used model of mechanical
trajectories cease to be linear as soon as there is a single
n response in solids. Yet formally it applies only in the
rearrangement (made or broken contact) in the contact
o limitofvanishinglyslowandweakdeformations. Inprac-
c ticematerialspossesscharacteristicstrainandtimescales network. Contact changes occur for perturbation am-
[ plitudes that vanish as 1/N, i.e. essentially immediately
that define a linear elastic “window”, i.e. a parameter
in large systems. Their findings caused the authors to
1 range wherein Hooke’s law is accurate. Determining the
question, if not the formal validity, then at least the use-
v size of this window is especially important in soft solids,
8
where viscous damping and nonlinearity play important
6
roles [2]. The goal of the present work is to determine
0
0 when Hooke’s law holds, and what eventually replaces 10−6
0 it, in packings of soft frictionless spheres close to the N=1024
. (un)jamming transition. p=10−4.5
1
0 Jammed sphere packings are a widely studied model 10−7
of emulsions and liquid foams [3–6] and have close con-
6
1 nections to granular media and dense suspensions [7–9].
v: Linear elastic properties of jammed solids, such as mod- sσ 10−8 γ†
i uli and the vibrational density of states, are by now well res
X understood [10, 11]. Much less is known about their vis- st
linearelasticity
ar coelastic [7, 12] and especially their nonlinear response 10−9 γ˙0=0+ (QS)
[13, 14]. Yet the jamming transition must determine the
γ γ˙0<γ˙†
sizeofthelinearelasticwindow,becausetheshearmodu- ∗
lusG0vanishescontinuouslyatthejammingpoint,where 10−10 γγ˙˙00≈>γγ˙˙††
theconfiningpressurepgoestozero. Indeed,recentstud-
ies of oscillatory rheology [15] and shocks [16–18] have 10−7 10−6 10−5 10−4
strainγ
shown that, precisely at the jamming point, any defor-
mation is effectively fast and strong, and neither viscous
effects nor nonlinearities can be neglected. FIG. 1. Ensemble-averaged stress-strain curves of packings
Becauseelasticityinfoams,emulsions,andotheramor- sheared at varying strain rate γ˙0. Hooke’s law predicts a
phousmaterialsresultsfromrepulsivecontactforces,mi- linearstress-straincurve(dashedline). Thecrossoverstrains
crostructuralrearrangementsofthecontactnetworkhave γ∗ andγ† areindicatedforthedatashearedatslowbutfinite
signaturesinthemechanicalresponse. Namely,theylead rate 0<γ˙0 <γ˙† (open circles).
2
fulness of linear elasticity in jammed solids – not just at Browniandisksthatinteractviaelasticandviscousforces
the jamming point, but anywhere in the jammed phase. when they overlap. Elastic forces are expressed in terms
Subsequently, Van Deen et al. [30] and Goodrich et of the overlap δ =1 r /(R +R ), where R and R
ij ij i j i j
−
al. [31, 32] argued that the situation is not so dire. They denote radii and (cid:126)r points from the center of particle i
ij
demonstrated that coarse grained properties of jammed to the center of j. The force is repulsive and acts along
solids are far less sensitive to contact changes than are the unit vector rˆ =(cid:126)r /r :
ij ij ij
the individual trajectories. Namely, under ensemble av-
(cid:40)
eragingtheslopeofthestress-straincurveisthesamebe- k(δ )δ rˆ , δ >0
f(cid:126)el = − ij ij ij ij (2)
fore and after the first contact change [30], and changes ij (cid:126)0, δ <0.
ij
in the density of states are negligible [32]. These re-
sultsshowthatlinearelasticconstitutiverelationsremain The prefactor k is the contact stiffness, which generally
useful near jamming, but they say nothing about when depends on the overlap
Hooke’s law eventually does break down.
Recentexperiments[13,21]andsimulations[14,24,33] k =k0δα−2. (3)
provide evidence for a two stage yielding process, where
Herek isaconstantandαisanexponentparameterizing
packingsfirstsoftenandonlylaterestablishsteadyshear 0
the interaction. In the following we consider harmonic
flow. Yet it remains unclear precisely how rate depen-
interactions (α = 2), which provide a reasonable model
dence, nonlinearity, and contact changes contribute to
for bubbles and droplets that resist deformation due to
the breakdown of linear elasticity. In order to unravel
surface tension; we also treat Hertzian interactions (α=
these effects, it is necessary to vary strain, strain rate,
5/2), which correspond to elastic spheres.
pressure, and system size simultaneously and systemati-
We perform simulations using two separate numerical
cally–aswedohereforthefirsttime. Usingsimulations
methods. The first is a molecular dynamics (MD) algo-
of viscous soft spheres, we find that Hooke’s law is valid
rithm that integrates Newton’s laws using the velocity-
withinasurprisinglynarrowwindowboundedbyviscous
Verlet scheme. Each disk is assigned a uniform mass
effects at small strain and nonlinear softening at large
m = πR2 proportional to its volume. Energy is dissi-
strain. The size of the linear elastic window displays i i
pated by viscous forces that are proportional to the rel-
power law scaling with pressure and correlates with the
ative velocity ∆(cid:126)vc of neighboring particles evaluated at
accumulationofnotone,butanextensivenumberofcon- ij
the contact,
tact changes.
The basic scenario we identify is illustrated in Fig. 1,
f(cid:126)visc = τ k(δ )∆(cid:126)vc , (4)
which presents ensemble-averaged stresses versus strain. ij − 0 ij ij
Shear is applied via a constant strain rate γ˙ at fixed
0 whereτ isamicroscopicrelaxationtime. Viscousforces
0
volume. We identify three characteristic scales, each of
can apply torques, hence particles are allowed to rotate
which depend on the initial pressure p: (i) For strains
as well as translate.
below γ∗ γ˙ τ∗, where τ∗ is a time scale, viscous
≡ 0 In addition to MD, we also perform simulations using
stresses are significant and Eq. (1) underestimates the
a nonlinear conjugate gradient (CG) routine [34], which
stress needed to deform the material. This crossover
keepsthesystematalocalminimumofthepotentialen-
strain vanishes under quasistatic shear (γ˙ 0, filled
squares). (ii) Above a strain γ† the materia0l s→oftens and ergylandscape,whichitselfchangesasthesystemunder-
goes shearing. The dynamics are therefore quasistatic,
Hooke’s law overestimates the stress. This crossover is
i.e. the particle trajectories correspond to the limit of
rate-independent,consistentwithplasticeffects. (iii)For
vanishing strain rate.
strain rates above a scale γ˙† (triangles), Eq. (1) is never
Bubble packings consist of N = 128 to 2048 disks in
accurate and there is no strain interval where the mate-
a 50:50 bidisperse mixture with a 1.4:1 diameter ratio.
rial responds as a linear elastic solid.
Shear is implemented via Lees-Edwards “sliding brick”
boundary conditions. The stress tensor is given by
I. SOFT SPHERES: MODEL AND 1 (cid:88) 1 (cid:88)
σ = f r m v v , (5)
BACKGROUND αβ 2V ij,α ij,β − V i i,α i,β
ij i
We first introduce the soft sphere model and summa- where V is the volume (area in two dimensions) of the
rizepriorresultsregardinglinearelasticitynearjamming. packing, f(cid:126) is the sum of elastic and viscous contact
ij
forces acting on particle i due to particle j, and(cid:126)v is the
i
velocity of particle i. Greek indices label components
A. Model along the Cartesian coordinates x and y. The confin-
ing pressure is p = (1/D)(σ +σ ), where D = 2 is
xx yy
−
We perform numerical simulations of the Durian bub- the spatial dimension, while the shear stress is σ = σ .
xy
ble model [4], a mesoscopic model for wet foams and The second term on the righthand side of Eq. (5) is a ki-
emulsions. The model treats bubbles/droplets as non- netic stress, which is always negligible in the parameter
3
ranges investigated here. Initial conditions are isotropic
withatargetedpressurep,preparedusingCGand“shear
a)
stabilized” in the sense of Dagois-Bohy et al. [35], which 101 b)
guaranteesthattheinitialslopeofthestress-straincurve c)
is positive. Stresses and times are reported in dimen-
100
sionless units constructed from k , τ , and the average
0 0
particle diameter. γ0
σ/ 10−1 γ0
=
Gr
B. Distance to jamming 10−2
γ0=4×10−7
φwdihs−Wtearφeenccuφesaentodatnhjeadexmcczmeosnisrnfiegmnf.eienraTgtnhopecrtoeehnsxestcuaerrcseetssppnveoucalmtusimbvaeeermvf∆raealazucsteui=sorenazto∆−fjφatmzh=ce-, 1100−−1430−1 1γγγ00000===444×××111100001−−−654 10N2p==11002−1440.53 104 105 106
c c t/τ
0
ming,arealsofrequentlyusedforthispurpose[10,36,37].
These three alternative order parameters are related via
FIG. 2. The ensemble-averaged relaxation modulus G at
r
p ∆φ ∆z2. (6) pressurep=10−4.5 forfourvaluesofthestrainamplitudeγ0.
k ∼ ∼ Inallfourcases,Gr displaysaninitialplateaucorresponding
to affine particle motion (inset a), followed by a power law
Here k should be understood as a typical value of the decayastheparticledisplacementsbecomeincreasinglynon-
contact stiffness in Eq. (3). The harmonic case (α = 2) affine (b). At long times the stress is fully relaxed and the
is straightforward because the contact stiffness is a con- final particle displacements are strongly non-affine (c).
stant. For other values of α, however, k depends on
the pressure. As the typical force trivially reflects its
bulk counterpart, f p, the contact stiffness scales as strain amplitude is small enough that the packing nei-
k f/δ p(α−2)/(α−∼1). Inthefollowing,allscalingrela- ther forms new contacts, nor breaks existing ones. They
∼ ∼
tionswillspecifytheirdependenceonkandthetimescale further found that the typical strain at the first contact
τ . In the present work τ is independent of the overlap change depends on pressure and system size as
0 0
betweenparticles(asintheviscoelasticHertziancontact
problem[38]),butweincludeτ becauseonecouldimag- (p/k)1/2
0 γ(1) . (8)
ine a damping coefficient kτ0 with more general overlap cc ∼ N
dependence than the form treated here.
Similar to the findings of Schreck et al. [25], this scale
vanishesinthelargesystemlimit,evenatfinitepressure.
C. Shear modulus and the role of contact changes
II. STRESS RELAXATION
In large systems the linear elastic shear modulus G
0
vanishes continuously with pressure,
We will characterize mechanical response in jammed
G /k (p/k)µ, (7) solids using stress relaxation and flow start-up tests, two
0
∼
standard rheometric tests. In the linear regime they are
with µ = 1/2. Hence jammed solids’ shear stiffness can equivalent to each other and to other common tests, in-
be arbitrarily weak. The scaling of G0 has been deter- cluding creep response and oscillatory rheology, as com-
mined multiple times, both numerically: [36, 39, 40] and pleteknowledgeoftheresultsofonetestpermitscalcula-
theoretically: [15, 41, 42]; it is verified for our own pack- tionoftheothers[2]. Thisequivalencebreaksdownonce
ings in Fig. 3a, as discussed in Section II. the response becomes nonlinear.
TherearetwostandardapproachestodeterminingG0. We employ stress relaxation tests to access the time
The first, which we employ, is to numerically impose a scaleτ∗ overwhichviscouseffectsaresignificant,andwe
small shear strain and relax the packing to its new en- use flow start-up tests to determine the strain scale γ†
ergyminimum[36,39]. Inthesecondapproachonewrites beyond which the stress-strain curve becomes nonlinear.
down the DN equations of motion and linearizes them We consider stress relaxation first.
about a reference state, which results in a matrix equa- In a stress relaxation test one measures the time-
tion that can be solved for the response to an infinitesi- dependent stress σ(t,γ ) that develops in a response to
0
mallyweakshear[15,35,40,42–44]. Thislatterapproach a sudden shear strain with amplitude γ , i.e.
0
allows access to the zero strain limit, but it is blind to
(cid:26)
the influence of contact changes. Van Deen et al. [30] 0 t<0
γ(t)= (9)
verifiedthatthetwoapproachesagree,providedthatthe γ t 0.
0
≥
4
The relaxation modulus is
N=128 N=256 N=512 N=1024 N=2048
G (t,γ ) σ(t,γ0). (10) 102 a) 10−1 b)
r 0
≡ γ0 10−2
We determine the relaxation modulus by employing the µ2N 101 µ λ2N10−3 −λ
0 /
shear protocol of Hatano [7]. A packing’s particles and G τ∗10 4
−
simulation cell are affinely displaced in accordance with
asimpleshearwithamplitudeγ0. E.g.forasimpleshear 10100 1100 101 102 103 104 10−150 1100 101 102 103 104
in the xˆ-direction, the position of a particle i initially at − −
pN2 pN2
(x ,y )instantaneouslybecomes(x +γ y ,y ),whilethe
i i i 0 i i
Lees-Edwards boundary conditions are shifted by γˆ0Ly, 102
where Ly is the height of the simulation cell. Then the c)
particles are allowed to relax to a new mechanical equi- γ0=4×10−6
librium while the Lees-Edwards offset is held fixed.
N=128...2048
moTdhuelimofaianspinagnleelpoafckFinigg. e2quililluibstrraatteedsaftouprrersesluarxeatpio=n G0101 p=10−5...10−2
/
10−4.5 and then sheared with strain amplitudes vary- Gr
ing over three decades. All four undergo a relaxation
from an initial plateau at short times to a final, lower θ
−
plateau at long times. The character of the particle mo-
100
tionschangesasrelaxationprogressesintime. Whilethe
particle motions immediately after the deformation are 10 4 10 3 10 2 10 1 100 101 102 103 104
− − − −
affine (Fig. 2a), they become increasingly non-affine as t/τ
∗
the stresses relax to a new static equilibrium (Fig. 2b,c).
This non-affine motion is a consequence of slowly relax-
FIG. 3. (a) Finite size scaling collapse of the linear shear
ing eigenmodes of the packing that become increasingly
modulus G in harmonic packings with exponent µ = 1/2.
0
abundant on approach to jamming [15]. These modes
(b) Finite size scaling collapse of the relaxation time τ∗ with
favor sliding motion between contacting particles [40],
exponent λ 1.13. (c) The relaxation modulus G collapses
r
reminiscent of zero energy floppy modes [45], and play to a master≈curve when G and t are rescaled with G and
r 0
an important role in theoretical descriptions of mechan- τ∗,respectively,asdeterminedin(a)and(b). Atshorttimes
ical response near jamming [15, 41, 42, 44, 46]. the master curve decays as a power law with exponent θ =
Forsufficientlysmallstrainamplitudes,linearresponse µ/λ 0.44 (dashed line).
≈
is obtained and any dependence of the relaxation mod-
ulus on γ is sub-dominant. The near-perfect overlap of
0
the moduli for the two smaller strain amplitudes Fig. 2 We showed in Fig. 2 that a packing relaxes in three
indicatesthattheyresideinthelinearregime. Thelong- stages. Theshort-timeplateauistrivial,inthesensethat
time plateau is then equal to the linear elastic modulus viscousforcespreventtheparticlesfromrelaxingatrates
G0. In practice there is a crossover time scale τ∗ such faster than 1/τ0; hence particles have not had time to
that for longer times t τ∗ viscous damping is negli- depart significantly from the imposed affine deformation
(cid:29)
gible and the relaxation modulus is well approximated and the relaxation modulus reflects the contact stiffness,
wbcryiollsitsdsoevateesrrymtmimipnteeoittseh,τeG∗sc≈ra(cid:39)l1in0Gg4τ0o0.f.τFIn∗orwthtihethefodplalroetwsasiuinnrgeF.Siegc.ti2oantwhee Gtimr e∼skc.alWesett(cid:29)herτe0f.ore focus hereafter on the response on
TodemonstratedynamiccriticalscalinginG ,wefirst
r
determinethescalingofitslong-timeasymptoteG . We
0
then identify the time scale τ∗ on which G significantly
r
A. Scaling in the relaxation modulus deviates from G . Finally, we show that rescaling with
0
these two parameters collapses the relaxation moduli for
We now characterize stress relaxation in linear re- a range of pressures to a single master curve. While we
sponse by measuring the relaxation modulus, ensemble- addressvariationswithstraininsubsequentSections,the
averaged over ensembles of packings prepared at varying strainamplitudehereisfixedtoavalueγ0 =10−5.5. We
pressure. We will show that G collapses to a critical have verified that this strain amplitude is in the linear
r
scalingfunctiongovernedbythedistancetothejamming regime for all of the data presented in this Section.
point, consistent with recent theoretical predictions by As noted above, at long times the relaxation modu-
Tighe [15]. Our main focus is on numerically measuring lus approaches the linear quasistatic modulus, G (t
r
→
the time scale beyond which viscous effects fade and the ) G . We verify the scaling for G from Eq. (7) in
0 0
∞ (cid:39)
responsebecomesquasistatic, whichispredictedtoscale ourharmonicpackingsbyrepeatingthefinitesizescaling
as τ∗ kτ /p. analysis of Goodrich et al. [47], who showed that finite
0
∼
5
size effects become important when a packing has O(1) 10−3
contacts in excess of isostaticity, or equivalently when
p/k 1/N2 – c.f. Eq. (6). Consistent with their results, 10−4
wefi∼ndthat G N2µ forvaryingN andpcollapsesto
0
a master curGve≡when plotted versus x pN2, as shown 10−5
≡
in Fig. 3a. The scaling of Eq. (7) is verified by this data
collapse together with the requirement for the modulus σ 10−6
s
to be an intensive property of large systems. To see this, es
note that G is intensive only if xµ for large x. str 10−7
0
G ∼
Again referring to Fig. 2, there is clearly some time
scale τ∗ such that for t < τ∗ the relaxation modulus 10−8 p=10−5.0 p=10−3.0
deviates significantly from the quasistatic modulus. To p=10−4.5 p=10−2.5
determine the scaling of τ∗ with p, we perform the fi- 10−9 p p=10−4.0 p=10−2.0
nite size scaling analysis presented in Fig. 3b. The re- p=10−3.5 slope1
laxation time is determined from the point where Gr, 10−1100−7 10−6 10−5 10−4 10−3 10−2 10−1 100
strainγ
averaged over an ensemble of at least 100 packings per
condition, has decayed to within a fraction ∆ of its fi-
nal value, G (t = τ∗) = (1+∆)G . We present data FIG.4. Averagedstress-straincurvesunderquasistaticshear
r 0
for ∆=1/e, but similar scaling results for a range of ∆ atvaryingpressurep. Solidanddashedcurveswerecalculated
[33]. Werequiretherescaledpressuretoremainx=pN2 using different strain protocols. Dashed curves: fixed strain
and collapse the data by rescaling the relaxation time as steps of 10−3, sheared to a final strain of unity. Solid curves:
τ∗/N2λ for a positive exponent λ. It follows that τ∗ di- logarithmicallyincreasingstrainsteps,beginningat10−9 and
verges in large systems near jamming as reaching a total strain of 10−2 after 600 steps.
(cid:18) (cid:19)λ
k
τ∗ τ as N . (11)
∼ p 0 →∞ III. FINITE STRAIN
We find the best data collapse for λ = 1.13, close to
When does linear elasticity break down under increas-
but somewhat higher than the value λ = 1 predicted by
ing strain, and what lies beyond? To answer these ques-
theory[15],althoughourcurrentnumericalresultsdonot
tions, we now probe shear response at finite strain using
exclude this possibility.
flow start-up tests.
We now use the linear quasistatic modulus G and
0
the characteristic time scale τ∗ to collapse the relax-
ation modulus to a master curve (s). Fig. 3c plots
G /G versus s t/τ∗ for aRrange of pressures A. Flow start-up
r 0
R ≡ ≡
and system sizes; data from the trivial affine regime at
times t < 10τ0 have been excluded. The resulting data Inaflowstart-uptest,strain-controlledboundarycon-
collapse is excellent, and the master curve it reveals has ditions are used to “turn on” a flow with constant strain
two scaling regimes: 1 for s 1, and s−θ rate γ˙ at time t=0, i.e.
R (cid:39) (cid:29) R ∼ 0
for s 1. The plateau at large s occurs by construc-
tion a(cid:28)nd corresponds to the quasistatic scaling Gr G0. (cid:26) 0 t<0
(cid:39) γ(t)= (13)
Thepowerlawrelaxationatshortertimescorrespondsto γ˙ t t 0
0
G G (t/τ∗)−θ for some exponent θ. By considering a ≥
r 0
∼
marginal solid prepared at the jamming point, one finds Toimplementflowstart-upinMD,attimet=0apack-
that the prefactor of t−θ cannot depend on the pressure. ing’sparticlesandsimulationcellareinstantaneouslyas-
InvokingthepressurescalingofG0 andτ∗ inthelargeN signed an affine velocity profile (cid:126)vi = (γ˙0yi,0)T in ac-
limit, identified above, we conclude that θ =µ/λ. Hence cordance with a simple shear with strain rate γ˙ ; the
0
in large systems the relaxation modulus scales as Lees-Edwards images of the simulation cell are assigned
Gr(t) (cid:26)(τ0/t)θ 1(cid:28)t/τ0 (cid:28)(k/p)λ (12) atocoemvomlveensaucrcaotredvinegloctiotyN. Tewhteonnt’shelapwasrtiwclheisleartehaellLoweeesd-
k ∼ (p/k)µ (k/p)λ (cid:28)t/τ0. Edwards boundary conditions maintain constant veloc-
with µ=1/2, λ 1.13, and θ =µ/λ 0.44. ity, so that the total strain γ(t) grows linearly in time.
≈ ≈
Anomalous stress relaxation with exponent θ 1/2 We also perform quasistatic shear simulations using
≈
was first observed in simulations below jamming [7] and nonlinearCGminimizationtorealizethelimitofvanish-
is also found in disordered spring networks [48, 49]. It is ing strain rate. Particle positions are evolved by giving
relatedviaFouriertransformtotheanomalousscalingof the Lees-Edwards boundary conditions a series of small
the frequency dependent complex shear modulus G∗ strainincrementsandequilibratingtoanewminimumof
(ıω)1−θ found in viscoelastic solids near jamming [15∼]. theelasticpotentialenergy. Thestressσisthenreported
WerevisitthescalingrelationofEq.(12)inSectionIIIF. as a function of the accumulated strain. For some runs
6
we use a variable step size in order to more accurately
101
determine the response at small strain.
p=10−5.0 p=10−3.0
Fig. 1 illustrates the output of both the finite strain p=10−4.5 p=10−2.5
rate and quasistatic protocols. p=10−4.0 p=10−2.0 N=1024
p=10−3.5
100
B. Quasistatic stress-strain curves
/Gγ0 1100−−21
To avoid complications due to rate-dependence, we σ 10−3
consider the limit of vanishing strain rate first. 10−1 γ†1100−−54 p
Fig. 4 plots the ensemble-averaged stress-strain curve 10−6
σ(γ) for harmonic packings at varying pressure. Pack- 1100−−87 ∆∆∆===000...765 ∆∆slo==pe00..=431
ings contain N = 1024 particles, and each data point is 10−190−6 10−5 10−4 10−3 10−2
averaged over at least 600 configurations. Several fea- p
tures of the stress-strain curves stand out. First, there 10−120−3 10−2 10−1 100 101 102
γ/p
is indeed a window of initially linear growth. Second,
beyond a strain of approximately 5 - 10% the system
achievessteadyplasticflowandthestress-straincurveis FIG.5. (mainpanel)DatafromFig.4,expressedasadimen-
sionlesseffectiveshearmodulusσ/G γandplottedversusthe
flat. Finally,theendoflinearelasticityandthebeginning 0
of steady plastic flow do not generally coincide; instead rescaledstrainγ/p. (inset)Thecrossoverstrainγ† wherethe
effective shear modulus has decayed by an amount ∆ in a
there is an interval in which the stress-strain curve has
system of N =1024 particles.
a complex nonlinear form. We shall refer to the end of
thelinearelasticregimeas“softening”becausethestress
initially dips below the extrapolation of Hooke’s law. (In
range 0.87 to 1.06, depending on the value of ∆. Bear-
the plasticity literature the same phenomenon would be
ingtheabovesubtletyinmind, weneverthelessconclude
denoted “strain hardening”.) Moreover, for sufficiently
that an effective power law with ν = 1 provides a rea-
low pressures there is a strain interval over which the
sonable description of the softening strain. Section IIA
stress increases faster than linearly. This surprising be-
presents further evidence to support this conclusion.
haviorisworthyoffurtherattention,butthefocusofthe
present work will be on the end of linear elasticity and
the onset of softening. This occurs on a strain scale γ†
D. Hertzian packings
that clearly depends on pressure.
In the previous section the pressure-dependence of γ†
was determined for harmonic packings. We now gener-
C. Onset of softening
alize this result to other pair potentials, with numerical
verification for the case of Hertzian packings (α=5/2).
Wenowdeterminethepressureandsystemsizedepen-
Recall that the natural units of stress are set by the
dence of the softening (or nonlinear) strain scale γ†.
contactstiffnessk,whichitselfvarieswithpressurewhen
Fig. 5 replots the quasistatic shear data from Fig. 4 α = 2. Based on the linear scaling of γ† in harmonic
(solid curves), now with the linear elastic trend G γ (cid:54)
0 packings, we anticipate
scaledout. Therescalingcollapsesdataforvaryingpres-
sures in the linear regime and renders the linear regime γ† p p1/(α−1), (14)
flat. The strain axis in Fig. 5b is also rescaled with the ∼ k ∼
pressure, a choice that will be justified below. The onset
of softening occurs near unity in the rescaled strain co- which becomes γ† p2/3 in the Hertzian case. To test
∼
ordinate for all pressures, which suggests that γ† scales thisrelation,werepeattheanalysisoftheprecedingSec-
linearly with p in harmonic packings (α=2). tion; results are shown in Fig. 6. We again find a finite
Unlike the linear relaxation modulus in Fig. 3c, the linear elastic window that gives way to softening. Soft-
quasistatic shear data in Fig. 5 do not collapse to a mas- ening onset can again be described with a ∆-dependent
ter curve; instead the slope immediately after softening exponent(seeinset). Itsvaluehasanarrowspreadabout
steepens (in a log-log plot) as the pressure decreases. As 2/3; power law fits give slopes between 0.63 and 0.74.
a result, it is not possible to unambiguously identify a
correlation γ† pν between the crossover strain and the
pressure. To c∼larify this point, the inset of Fig. 5 plots E. Relating softening and contact changes
the strain where σ/G γ has decayed by an amount ∆
0
from its plateau value, denoted γ†(∆). This strain scale Whydoesthelinearelasticwindowclosewhenitdoes?
is indeed approximately linear in the pressure p (dashed We now seek to relate softening with contact changes on
curves), but a power law fit gives an exponent ν in the the particle scale [21–25, 30]. Specifically, we identify a
7
100 100
101
pppp====11110000−−−−7665....0505 ppp===111000−−−544...050 N=1024 ncc1100−−21 a) p=10−4 ncc111000−−−321 bp) pp==1100−−54.5
/Gγ0 100 1100−−21 1100−−1430 610 510 410 3NNNN10====1251251028622140 1 111000−−−16540 N610=5110024410 31ppppp=====01111100000−−−−−2433221..550 1
σ 10−3 − − − − − − − − − − − −
10−1 γ†111000−−−654 p 103 γ γ
1100−−87 ∆∆∆===000...654 ∆slo=pe0.=32/3 102 c)
10−190−8 10−7 10−6 10−5 10−4 10−3
p 101
10−120−3 10−2 10−1 100 101 102
γ/p2/3 τp 100
/
FIG. 6. (main panel) The dimensionless shear modulus of ncc10−1 τ=0.5
quasistatically sheared Hertzian packings plotted versus the 10 2 p=10−5...10−2
rescaled strain γ/p2/3. (inset) Pressure-dependence of the − N=128...1024
crossover strain γ†. 10−3
10 4
−10 4 10 3 10 2 10 1 100 101 102 103 104
− − − −
correlation between the softening strain γ†, the cumula- γ/p
tive number of contact changes, and the distance to the
isostatic contact number z . In so doing we will answer
c
FIG. 7. The contact change density shown for (a) varying
the question first posed by Schreck and co-workers [25],
system size and (b) varying pressure. (c) Data collapse for
who asked how many contact changes a packing can ac- pressures p = 10−2...10−5 in half decade steps and system
cumulate while still displaying linear elastic response.
sizesN =128...1024inmultiplesof2. Dashedlinesindicate
We begin by investigating the ensemble-averaged con- slopes of 1 and 1/2.
tact change density n (γ) [N (γ)+N (γ)]/N,
cc make break
≡
whereN andN arethenumberofmadeandbro-
make break
kencontacts,respectively,accumulatedduringastrainγ. tweencontactchangesatagivenstrain. Hencetheinitial
Contact changes are identified by comparing the contact slope of n is fixed by γ(1):
cc cc
network at strain γ to the network at zero strain.
(cid:32) (cid:33)
InFig.7aweplotn forpackingsofharmonicparticles 1 γ
cc n (γ) as γ 0. (15)
at pressure p=10−4 and varying system size. The data cc (cid:39) N γ(1) →
cc
collapse to a single curve, indicating that n is indeed
cc
an intensive quantity. The effect of varying pressure is From Fig. 7 it is apparent that ncc remains linear in γ
shown in Fig. 7b. There are two qualitatively distinct up to the crossover strain γ†. We conclude that γ(1)
cc
regimes in ncc, with a crossover governed by pressure. describes the strain between successive contact changes
To better understand these features, we seek to col- over the entire interval 0 γ < γ†. In the softening
≤
lapse the n data to a master curve. By plotting regime the strain between contact changes increases; it
cc
N ≡ncc/pτ versus y ≡γ/p, we obtain excellent collapse scales as ncc ∼γ1/2 (see Fig. 7c).
for τ = 1/2, as shown in Fig. 7b for the same pressures Let us now re-interpret the softening crossover strain
as in Fig. 7a and system sizes N = 128...1024. The γ† ∆z2 (c.f. Eq. (6)) in terms of the coordination of
∼
scaling function y for small y, while yτ for the contact network. We recall that ∆z = z zc is the
y > 1. The rescaNled∼strain y provides furthNer∼evidence difference between the initial contact number−z and the
for∼a crossover scale γ† p/k, now apparent at the mi- isostatic value zc, which corresponds to the minimum
∼
croscale. Moreover,thefactthatdataforvaryingsystem number of contacts per particle needed for rigidity. The
sizes all collapse to the same master curve is an impor- excess coordination ∆z is therefore an important char-
tant indicator that γ† is an intensive strain scale that acterization of the contact network. The contact change
remains finite in the large system size limit. density at the softening crossover, n† , can be related to
cc
The scaling collapse in Fig. 7c generalizes the re- ∆z via Eq. (15), while making use of Eq. (6),
sults of Van Deen et al. [30], who determined the strain
n† n (γ†) ∆z. (16)
scale γ(1) (p/k)1/2/N associated with the first con- cc ≡ cc ∼
cc
∼
tact change. To see this, note that the inverse slope Hence we have empirically identified a topological crite-
(dγ/dn )/N represents the average strain interval be- rionfortheonsetofsoftening: aninitiallyisotropicpack-
cc
8
strains, however, the effective shear modulus is stiffer
102 than the quasistatic curve and decays as σ/γ t−θ
γ˙0 γγ˙˙00==1100−−1110 (see inset). This is rate-dependence: for a given∼strain
γ˙0=10−9 amplitude, the modulus increases with increasing strain
101 γ˙0=10−8 rate. Correspondingly, thecharacteristicstrainγ∗ where
γ˙0=10−7 curves in the main panel of Fig. 8 reach the linear elas-
QS
tic plateau (σ/G γ 1) grows with γ˙ . For sufficiently
0 0
γ ≈
G0 100 highstrainratesthereisnolinearelasticplateau;forthe
σ/ 102 data in Fig. 8 this occurs for γ˙0 10−8. Hence there is
101 a characteristic strain rate, γ˙†, b≈eyond which the linear
10−1 σ/γ101−010 −−−−−θθθθθ eγ˙l†asatriecawlwinadyoswrahtaes-dcelposeendd:enptaacknidn/gosrsshteraarinedsofaftsetnerintgh.an
10−1200101102103104105106 Np==1100−244 strTaoinsu,nwdeerrsetvainsdittthhee rrealtaex-daetipoenndmeondturleusspdoentseermatinsemdainll
10−120−9 10−8 10−7t 10−6 10−5 10−4 10−3 10−2 SectionII.Inlinearresponsethestressafterflowstart-up
γ depends only on the elapsed time t=γ/γ˙ ,
0
FIG. 8. The effective shear modulus during flow start-up for σ 1 (cid:90) t
= G (t(cid:48))dt(cid:48). (17)
packings of N =1024 particles at pressure p=10−4, plotted γ t r
0
versus strain for varying strain rates γ˙ . (inset) The same
0
datacollapsesforearlytimeswhenplottedversust,decaying Employing the scaling relations of Eq. (12), one finds
as a power law with exponent θ=µ/λ 0.44 (dashed line).
≈
σ (cid:16)τ (cid:17)θ
k 0 , τ <t<τ∗, (18)
0
γ ∼ t
ing softens when it has undergone an extensive number
of contact changes that is comparable to the number of as verified in Fig. 8 (inset). Linear elasticity σ/γ G0
contacts it initially had in excess of isostaticity. (This is only established at longer times, when γ > γ˙0(cid:39)τ∗
does not mean the packing is isostatic at the softening (k/p)λγ˙0τ0. Hence the relaxation time τ∗ plays a∼n
crossover,asn countsbothmadeandbrokencontacts.) important role: it governs the crossover from rate-
cc
dependent to quasistatic linear response. The system
requiresatimeτ∗ torelaxafteraperturbation. Whenit
is driven at a faster rate, it cannot relax fully and hence
F. Rate-dependence
its response depends on the driving rate.
We can now identify the characteristic strain rate γ˙†
To this point we have considered nonlinear response
where the linear elastic window closes. This rate is
exclusively in the limit of quasistatic shearing. A mate- reachedwhentheboundonquasistaticity, γ >γ˙ τ∗, col-
0
rial accumulates strain quasistatically when the imposed lides with the bound on linearity, γ <γ†, giving
strain rate is slower than the longest relaxation time
in the system. Because relaxation times near jamming (p/k)1+λ
are long and deformations in the lab always occur at fi- γ˙† , (19)
∼ τ
0
nite rate, we can anticipate that quasistatic response is
difficult to achieve and that rate-dependence generically with 1 + λ 2.1. This strain rate vanishes rapidly
≈
playsasignificantrole. Henceitisimportanttoconsider nearjamming,andpackingsmustbeshearedincreasingly
shear at finite strain and finite strain rate. We now con- slowlytoobserveastress-straincurvethatobeysHooke’s
sider flow start-up experiments in which a finite strain law. As a practical consequence, experiments near jam-
rate γ˙ is imposed at time t=0, cf. Eq. (13). ming are unlikely to access the linear elastic regime.
0
Fig. 8 displays the mechanical response to flow start-
upforvaryingstrainrates. Tofacilitatecomparisonwith
thequasistaticdataoftheprevioussection,flowstart-up IV. DISCUSSION
data are plotted in terms of the dimensionless quantity
σ(t;γ˙0)/G0γ,whichweshallrefertoastheeffectiveshear Usingacombinationofstressrelaxationandflowstart-
modulus. The data are for systems of N = 1024 parti- upexperiments,wehaveshownthatsoftsolidsnearjam-
cles, averaged over an ensemble of around 100 realiza- ming are easily driven out of the linear elastic regime.
tions each. Here we plot data for the pressure p=10−4; There is, however, a narrow linear elastic window that
results are qualitatively similar for other pressures. For survivestheaccumulationofanextensivenumberofcon-
comparison, we also plot the result of quasistatic shear tact changes. This window is bounded from below by
(solid circles) applied to the same ensemble of packings. rate-dependent viscoelasticity and bounded from above
Packings sheared sufficiently slowly follow the qua- by the onset of strain softening. Close to the transition
sistatic curve; see e.g. data for γ˙ = 10−11. For smaller these two bounds collide and the linear elastic window
0
9
˙γ areforfrequencieshigherthanγ˙†, whereviscousstresses
rate dominate. There are also qualitative differences between
n
strai ellainsteiacrity thequasistaticshearmodulus,whichcannotbecollapsed
to a master curve (Fig. 5), and the storage modulus in
k γ† oscillatory shear, which can [14, 33]. We speculate that
/
p there are corresponding microstructural differences be-
re γ∗ tween packings in steady state and transient shear [20],
ssuγ˙†1/2 similar to those which produce memory effects [50].
e
r
p Soft sphere packings near jamming approach the iso-
static state, which also governs the rigidity of closely
related materials such as biopolymer and fiber net-
works [51–54]. It is therefore remarkable to note that,
strainγ whereas sphere packings soften under strain, quasistat-
ically sheared amorphous networks are strain stiffening
FIG. 9. In a flow start-up test, quasistatic linear response beyond a crossover strain that scales as ∆z [55], which
(G G0) occupies a strain window γ∗ <γ <γ† (shaded re- vanishes more slowly than γ† ∆z2 in packings. Hence
≈ ∼
gions). For smaller strains the response is rate-dependent, nonlinearity sets in later and with opposite effect in net-
with a crossover strain γ∗ that depends on both pressure works [56]. We expect that this difference is attributable
and strain rate. Softening sets in for higher strains, with
to contact changes, which are absent or controlled by
a crossover γ† that depends only on the pressure. The inter- slow binding/unbinding processes in networks.
sectionoftherate-dependentandsofteningcrossoversdefines
We have demonstrated that the onset of softening oc-
a strain rate γ˙† above which there is no quasistatic linear
response, i.e. the shaded region closes. curs when the system has accumulated a finite number
of contact changes correlated with the system’s initial
distance from the isostatic state. This establishes an im-
closes. Finally, weakly jammed materials are generally portantlinkbetweenmicroscopicandbulkresponse. Yet
rate-dependentand/orstrainsofteningonscalesrelevant further work investigating the relationship between mi-
to the laboratory, because the strains and strain rates croscopicirreversibility,softening,andyieldingisneeded.
bounding the linear elastic window vanish rapidly near The inter-cycle diffusivity in oscillatory shear, for exam-
jamming. Fig. 9 provides a qualitative summary of our ple, jumps at yielding [21, 24], but its pressure depen-
results. dence has not been studied. Shear reversal tests could
Whileoursimulationsareintwodimensions,weexpect alsoprovideinsightintotheconnectionbetweenjamming
thescalingrelationswehaveidentifiedtoholdforD >2. and plasticity.
To the best of our knowledge, all scaling exponents near While the onset of softening can be probed with qua-
jammingthathavebeenmeasuredinboth2Dand3Dare sistatic simulation methods, rate dependent effects such
the same. There is also numerical evidence that D = 2 as the strain scale γ∗ should be sensitive to the manner
is the transition’s upper critical dimension [32, 47]. in which energy is dissipated. The dissipative contact
Our work provides a bridge between linear elasticity forces considered here are most appropriate as a model
near jamming, viscoelasticity at finite strain rate, and for foams and emulsions. Hence useful extensions to the
nonlinearityatfinitestrainamplitude. Themeasuredre- present work might consider systems with, e.g., lubrica-
laxationmodulusG isingoodagreementwiththelinear tion forces or a thermostat.
r
viscoelasticity predicted by Tighe [15]. Consistent with
the granular experiments of Coulais et al., we identify a
crossover to nonlinear strain softening. Their crossover
V. ACKNOWLEDGMENTS
scales differently with the distance to jamming, possibly
due to the presence of static friction. The emulsions of
Knowltonetal.alsosoften[21]. Theydisplayacrossover WethankP.Boukany,D.J.Koeze,M.vanHecke,and
strain that is roughly linear in ∆φ, consistent with both S. Vasudevan for valuable discussions. JB, DV and BPT
ourㆠandtheresultsofOtsukiandHayakawa[14], who were supported by the Dutch Organization for Scientific
simulated large amplitude oscillatory shear at finite fre- Research(NWO).ESwassupportedbytheJ´anosBolyai
quency. The agreement between the crossover strains in Research Scholarship of the Hungarian Academy of Sci-
ourquasistaticsimulationsandtheoscillatoryshearsim- ences. This work was carried out on the Dutch national
ulations of Ref. [14] is surprising, as most of their results e-infrastructure with the support of SURF Cooperative.
[1] L. D. Landau and E. M. Lifshitz, Theory of Elasticity [2] H. A. Barnes and J. F. Hutton, An Introduction to Rhe-
(Butterworth-Heineman, Oxford, 1997). ology (Elsevier, 1989).
10
[3] F. Bolton and D. Weaire, Phys. Rev. Lett. 65, 3449 (2014).
(1990). [31] C. P. Goodrich, A. J. Liu, and S. R. Nagel, Phys. Rev.
[4] D. J. Durian, Phys. Rev. Lett. 75, 4780 (1995). Lett. 112, 049801 (2014).
[5] B.P.Tighe,E.Woldhuis,J.J.C.Remmers,W.vanSaar- [32] C. P. Goodrich, S. Dagois-Bohy, B. P. Tighe, M. van
loos, and M. van Hecke, Phys. Rev. Lett. 105, 088303 Hecke, A. J. Liu, and S. R. Nagel, Phys. Rev. E 90,
(2010). 022138 (2014).
[6] G.Katgert,B.P.Tighe, andM.vanHecke,SoftMatter [33] S. Dagois-Bohy, E. Somfai, B. P. Tighe, and M. van
9, 9739 (2013). Hecke, (in preparation) (2014).
[7] T. Hatano, Phys. Rev. E 79, 050301 (2009). [34] D. V˚agberg, P. Olsson, and S. Teitel, Phys. Rev. E 83,
[8] J. R. Seth, L. Mohan, C. Locatelli-Champagne, 031307 (2011).
M. Cloitre, and R. T. Bonnecaze, Nat Mater 10, 838 [35] S. Dagois-Bohy, B. P. Tighe, J. Simon, S. Henkes, and
(2011). M. van Hecke, Phys. Rev. Lett. 109, 095703 (2012).
[9] S. V. Franklin and M. D. Schattuck, eds., Handbook of [36] C. S. O’Hern, L. E. Silbert, A. J. Liu, and S. R. Nagel,
Granular Materials (CRC Press, 2015). Phys. Rev. E 68, 011306 (2003).
[10] M. van Hecke, J. Phys. Cond. Matt. 22, 033101 (2010). [37] G. Katgert and M. van Hecke, EPL 92, 34002 (2010).
[11] A.J.LiuandS.R.Nagel,Ann.Rev.Cond.Matt.Phys. [38] R.Ram´ırez,T.P¨oschel,N.V.Brilliantov, andT.Schwa-
1, 347 (2010). ger, Phys. Rev. E 60, 4465 (1999).
[12] D. A. Head, Phys. Rev. Lett. 102, 138001 (2009). [39] H. P. Zhang and H. A. Makse, Phys. Rev. E 72, 011301
[13] C.Coulais,A.Seguin, andO.Dauchot,Phys.Rev.Lett. (2005).
113, 198001 (2014). [40] W.G.Ellenbroek,E.Somfai,M.vanHecke, andW.van
[14] M. Otsuki and H. Hayakawa, Phys. Rev. E 90, 042202 Saarloos, Phys. Rev. Lett. 97, 258001 (2006).
(2014). [41] M. Wyart, Annales de Physique 30, 1 (2005).
[15] B. P. Tighe, Phys. Rev. Lett. 107, 158303 (2011). [42] A. Zaccone and E. Scossa-Romano, Phys. Rev. B 83,
[16] L.R.Go´mez,A.M.Turner,M.vanHecke, andV.Vitelli, 184205 (2011).
Phys. Rev. Lett. 108, 058001 (2012). [43] L.E.Silbert,A.J.Liu, andS.R.Nagel,Phys.Rev.Lett.
[17] S. Ulrich, N. Upadhyaya, B. van Opheusden, and 95, 098301 (2005).
V. Vitelli, PNAS 110, 20929 (2013). [44] M. Wyart, S. R. Nagel, and T. A. Witten, Euro-
[18] S. van den Wildenberg, R. van Loo, and M. van Hecke, phys. Lett. 72, 486 (2005).
Phys. Rev. Lett. 111, 218003 (2013). [45] S. Alexander, Phys. Rep 296, 65 (1998).
[19] M. Lundberg, K. Krishan, N. Xu, C. S. O’Hern, and [46] C. Maloney, Phys. Rev. Lett. 97, 035503 (2006).
M. Dennin, Phys. Rev. E 77, 041505 (2008). [47] C. P. Goodrich, A. J. Liu, and S. R. Nagel, Phys. Rev.
[20] I. Regev, T. Lookman, and C. Reichhardt, Phys. Rev. Lett. 109, 095704 (2012).
E 88, 062401 (2013). [48] B. P. Tighe, Phys. Rev. Lett. 109, 168303 (2012).
[21] E.D.Knowlton,D.J.Pine, andL.Cipelletti,SoftMat- [49] M. Sheinman, C. P. Broedersz, and F. C. MacKintosh,
ter 10, 6931 (2014). Phys. Rev. E 85, 021801 (2012).
[22] N. C. Keim and P. E. Arratia, Phys. Rev. Lett. 112, [50] N.C.Keim,J.D.Paulsen, andS.R.Nagel,Phys.Rev.
028302 (2014). E 88, 032306 (2013).
[23] N. C. Keim and P. E. Arratia, Soft Matter 11, 1539 [51] C. Heussinger and E. Frey, Phys. Rev. Lett. 96, 017802
(2015). (2006).
[24] T. Kawasaki and L. Berthier, arXiv:1507.04120 (2015). [52] C. Heussinger and E. Frey, Phys. Rev. Lett. 97, 105501
[25] C.F.Schreck,T.Bertrand,C.S.O’Hern, andM.Shat- (2006).
tuck, Phys. Rev. Lett. 107, 078301 (2011). [53] C. P. Broedersz, X. Mao, T. C. Lubensky, and F. C.
[26] C. F. Schreck, T. Bertrand, C. S. O’Hern, and M. D. MacKintosh, Nat. Phys. 7, 983 (2011).
Shattuck, arxiv:1306.1961 (2013). [54] M. Das, D. Quint, and J. Schwarz, PloS One 7, e35939
[27] C. F. Schreck, R. S. Hoy, M. D. Shattuck, and C. S. (2012).
O’Hern, Phys. Rev. E 88, 052205 (2013). [55] M.Wyart,H.Liang,A.Kabla, andL.Mahadevan,Phys.
[28] C.Schreck,C.O’Hern, andM.Shattuck,GranularMat- Rev. Lett. 101, 215501 (2008).
ter 16, 209 (2014). [56] B. P. Tighe, Granular Matter 16, 203 (2014).
[29] T. Bertrand, C. F. Schreck, C. S. O’Hern, and M. D.
Shattuck, Phys. Rev. E 89, 062203 (2014).
[30] M. S. van Deen, J. Simon, Z. Zeravcic, S. Dagois-Bohy,
B.P.Tighe, andM.vanHecke,Phys.Rev.E90,020202