Table Of ContentGranular Matter manuscript No.
(will be inserted by the editor)
A contact model for the yielding of caked granular materials
L. Brendel · J. To¨ro¨k · R. Kirsch · U. Bro¨ckel
1
1
0
2
n
a
J
7 Received: date /Accepted: date
2
] Abstract We present a visco-elastic coupling model tractionor liquidmenisci) andthose which vanishirre-
h
between caked spheres, suitable for DEM simulations, versiblywhenbeingovercome(likesinternecksorother
c
e which incorporates the different loading mechanisms solidbridges). In this work,we focus solely onthis last
m
(tension,shear,bending,torsion)inacombinedmanner type.
- and allows for a derivation of elastic and failure prop-
t Such cemented but breakable contacts form (delib-
a erties on a common basis. In pull, shear, and torsion
t erately or accidentally) in agglomerates of particles.
s failuretestswithagglomeratesofupto10000particles,
The failure behavior of these agglomerates is of prac-
.
t wecomparethefailurecriteriontodifferentapproxima-
a tical interest in both cases, but also from the point of
tive variants of it, with respect to accuracy and com-
m view of Distinct Element Methods (DEMs), modelling
putational cost. The failure of the agglomerates,which
- the hindering of all six relative degrees of freedom of
d behaveaccordingtoelasticparametersderivedfromthe
two particles in contact is much less studied than the
n contactelasticity,givesalsoinsightintotherelativerel-
wellestablishedfrictionalcontactmodelslikethelinear
o
evance of the different load modes.
c spring-dash pot, the Hertz contact and its extensions
[ (cf. e.g. [12,8]).
Keywords DEM simulations cemented contacts
1 contact failure model · · For modelling purposes two aspects of the contacts
v
have to be taken into account: its elastic response to
0
3 deformations and the maximal load it can bear. The
3 1 Introduction elastic properties are usually expressed as set of (often
5
decoupled) springs associated to each relative degree
.
1 The research of granular media can be separated into of motion (cf. e.g. [4,1,15]). Also for the maximal load
0
tworatherdistinctfields:Therearesystemswithpurely it is a common practice to employ separate thresholds.
1
1 repulsive granular interaction (“dry granular media”, Exceptions(inthecontextofgranularassemblies)toan
: e.g. [5]) and media with non-negligible additional at- independenttreatmentofjusttensileandshearloadare
v
i tractions(“wetgranularmedia”,e.g.[9]).Whileforthe [3,4,15]whoemployempiricalorplausiblecombination
X former the further distinction between granular gases functions.
ar (e.g. [11]) and dense systems can be made, the latter It is thus legitimate to use a simple but self consis-
is more adequately separated into interactions of per- tent method to model the cemented contact by a flat,
manent or recurringnature (like e.g. van der Waals at- elasticcylindersubjectedtotheTrescafailurecriterion
(section 2). The full Tresca criterion is only solvable
L. Brendel,J. T¨oro¨k
with iterations therefore we compare the full criterion
FacultyofPhysics,UniversityofDuisburg-Essen,47048Duis-
to its two simplified variants in DEM simulations of
burg, Germany
E-mail: [email protected] the quasi-brittlefailure ofacylinder made upof200to
10000particles(section3).Weconcludewithdiscussing
R. Kirsch,U. Bro¨ckel
Fachhochschule Trier, Umwelt-Campus Birkenfeld 55761 the relative relevance of the different loading modes in
Birkenfeld,Germany section 4.4.
2 L. Brendel etal.
2 Modelling breakable cemented contacts From these fields, we can immediately read off the de-
formation geometry as
2.1 Contact model 2d(1+ν)
X =e u (0,0,d)= F (5)
x· shear Eπa2 x
2d(1+ν)
Y =e u (0,0,d)= F (6)
The principle modelling idea for the cemented contact y · shear Eπa2 y
between two particles is the dominance of the bridge d
Z =e u (0,0,d)= F (7)
withrespecttotheelasticandfailureproperties.Hence, z · push/pull Eπa2 z
weenvisagethebridgeasaflatcylinder(heightdsmaller 4d(1+ν)
φ=e u (a,0,d)/a= T (8)
thandiameter2a)whiletreatingtherestofthetwopar- y · torsion Eπa4 z
ticles as solid bodies. 4d
θ = ∂ (e u (0,0,d))= T . (9)
− z y· bend Eπa4 x
We set up a local coordinate frame, where the bot-
Summing up (1) to (4), calculating the correspond-
tom “plate” of the cylinder is fixed in the x/y-plane
ingstrainandapplyingHook’slawyieldsastresstensor
with its center in the origin, while the upper one is lo-
of
cated (in the unstrained state) parallel at z = d (i.e.
the original contact normal is in the z direction). We
σ =σ =σ =σ =0 (10)
xx xy xy yy
call the torsion about the z-axis “twist” with angle φ,
a2F 2T y
whilewechoosethex-axisofthelocalframeastheaxis σ =σ = x− z (11)
of the other rotational load, the “tilt” θ of the upper xz zx πa4
a2F +2T x
plate. These rotations (being, due to the brittle limit, y z
σ =σ = (12)
smallastoallowforcommutation)areaccompaniedby yz zy πa4
the translational displacements, the pull/push Z and a2Fx+4Txy
σ = (13)
the shear X, Y. In total, we describe the load geomet- zz πa4
rically by these five parameters.This is independent of whichobeysvanishingdivergenceandvanishingnormal
the actual representation of particle orientation e.g. in stress at the free surface. Integration over the lateral
a DEM code. surfacesconfirmstheforceF=F e +F e +F e and
x x y y z z
reveals a torque (width respect to the point e d/2) of
z
Next,weneedtoknowtheelasticresponseinterms
T=(T +F d/2)e +F d/2e +T e onthebottomone
x x x y y z z
of pull/push force F , shear forces F ,F , and torques
z x y (and thus F and T on the top one, of course). The
T ,T opposingthetwistandtilt,respectively.Weavoid − −
z x additional torque contributions are indeed necessary,
introducing five arbitrary spring constants, and aim to
sincewedefinedtheshearmodeaskeepingbottomand
derive them from the elastic behavior of the cylinder,
top plate parallel.
instead.
With the definition of the spring constant k
n
≡
Eπa2/d (which in principle can be measured in exper-
The stress inside a deformed cylinder (possessing
iments with two particles [6]), the relation between de-
elastic modulus E and Poisson ratio ν) is in principle
formation and reaction force finally reads
a text-book problem (cf. e.g. [13]). We directly provide
thedisplacementfieldassociatedtoeachloadingmode:
F 1 0 0 0 0
x 2+2ν X
Fy 0 2+12ν 0 0 0 Y
F 0 0 1 0 0
F TTTxyz///zaaa=−kn2+1200ν2da 2+1200ν2da 000 0041 4+0014νaaZφθ .
upush/pull(x,y,z)= Eπza2 (−νxex−νxey +zez) (1) (14)
4T (1+ν)
u (x,y,z)= z ( yze +xze ) (2) When considering the cylinder (the geometry of which
torsion Eπa4 − x y doesnotexplicitlyenterintothecontactgeometryany-
2T
u (x,y,z)= x ( 2νxye +2νyze way)as being very flat (i.e. d/a 0), the stiffness ma-
bend Eπa4 − x z trix looses the T -line and becom→es diagonal. This de-
y
(ν(y2 x2)+z2)ey (3) coupling is very convenient for DEM modelling. More-
− −
2+2ν (cid:1) over,everyforce/torquecomponentcanbesupplemented
u (x,y,z)= (F ze +F ze ) (4)
shear Eπa2 x x y y by a viscous damping involving the time derivative of
Acontact model for theyieldingof caked granular materials 3
thedeformationandwithcoefficientsforcriticaldamp- where Φ2 and 4Φ¯2 differ only by an irrelevant x¯/y¯-
ing: independent term. Obviously, Φ¯2 is just a quadratic
form in x¯ and y¯ the center and steepness of which de-
Fx =2 ktm¯ X˙ +ktX (15) pend on the load. Unfortunately, when maximizing Φ¯2
−
F =2pk m¯ Y˙ +k Y (16) subjected to the constraint x¯2 + y¯2 1 (which, due
− y p t t to the absence of local maxima, can≤be replaced by
F =2 k m¯ Z˙ +k Z (17)
− z n n x¯2+y¯2 =1) one is facedwith a quartic,the radicalsof
T =apk I¯ θ˙+ kna2 θ (18) which exceed any sensible formula length.
x n x
− 4 Forseveralspecialcases,thecalculationcanbeper-
p
T =a 2k I¯ φ˙ + kta2 φ (19) formed, though: In the simplest case of zero torques,
z t z
− 2 σ becomes independent of x and y and thus
p Tresca
Here, k k /(2+2ν) while m¯, I¯ , and I¯ are the re-
t n x z
≡ πa2σ F = F2+F2+F2/4 (26)
duced mass and moments of inertia of the two-particle Tresca,max ≡ Tresca q x y z
system, respectively. Critical damping is used for con-
venience when the deformation rate is anyway slow and can be compared directly to the threshold F∗
enough as to neglect velocity effects. πa2σ∗. A less simple case is Fy = 0 which lets Φ¯2 tak≡e
on its maximum at y¯=sgn(F T 2F T ), producing
z x x z
−
the value
2.2 Failure criterion
F2 F T 2F T T2+T2
Having worked out the elastic response of the bridge, FTresca =rFx2+ 4z +2| z x−a x z| +4 xa2 z .
we now turn to its failure. The choice of failure criteria
(27)
for condensed matter to be found in the literature is
vast.We will use the wellknownTrescacriterion[14] of
Yet another illuminating combination, namely F =
x
maximal shear stress. From the stress tensor (10)-(13)
F =0, yields
z
the principle stresses are readily found to be
(F +2T¯ x¯)2
πa4σ =a2F /2+2T y a2√Φ2 (20) Φ¯2 = y z +(T¯2+T¯2)y¯2 (28)
1 z x − 4 x z
πa4σ =0 (21)
2
πa4σ =a2F /2+2T y+a2√Φ2 (22) which possesses a local maximum at x¯ = FyTz/(2Tx2),
3 z x
y¯2 =1 x¯2 (thesignofy¯beingirrelevant).Thisisvalid
only fo−r F T¯ <2T¯2, otherwise x¯=sgn(F T¯ ), y¯=0
where | y z| x y z
has to be chosen. In the latter case, the value taken on
a4Φ2 =a4(F2+F2+F2/4)+4(T2y2+T2(x2+y2)) is
x y z x z
+2a2(FzTxy+2Tz(Fyx−Fxy)) . (23) F = F2+4|FyTz| +4Tz2 (29)
Tresca r y a a2
Since σ < σ < σ 1, the maximal shear stress (with
1 2 3
respect to orientation) is
For the general case (occurring naturally in simu-
lations), we resorted to locating the maximum of (25)
σ σ √Φ2
σ (x,y)= 3− 1 = (24) numerically (by means of a Lagrangian multiplier, the
Tresca 2 πa2
Newton-Raphsonmethodandbeingcarefulaboutpoles).
andhas yetto be maximized with respectto x2+y2 Let us term this and the subsequent usage of the max-
≤
a2 and then compared to a material dependent critical imized Φ in (24) the application of the full (Tresca-
stressσ∗.That means,we areleft with a maximization based) criterion. In section 4, we will come back to it
of Φ2. Introducing forces T¯ = T /a, T¯ = T /a and and compare the results to the ones obtained by em-
x x z z
dimensionless position x¯=x/a, y¯=y/a yields ploying the following rather crude estimate.
To obtain an approximateclosedformula, we disre-
F 2 F T¯ 2F T¯ 2 gard the exact constraint x¯2 +y¯2 = 1 and replace it
Φ¯2 =T¯z2(cid:18)x¯+ 2T¯y (cid:19) +(T¯x2+T¯z2)(cid:18)y¯+ 4z(Tx¯2−+T¯x2)z(cid:19) by x¯2 = y¯2 = 1. Then again, the lack of local maxima
z x z in Φ¯2 enforces its maximum to be locatedat one of the
(25)
fourcornerpoints(x¯,y¯)=( 1, 1).Thisinturn,given
± ±
1 We usenegative signfor compressivestress. the form of (25), shows that the signs of x¯ and y¯ have
4 L. Brendel etal.
to be those of F /T¯ and F T¯ 2F T¯ , respectively.
y z z x x z
−
Plugging that back into (24) yields
F˜2 =F2+F2+F2/4+4(T2+2T2)/a2
Tresca x y z x z
+(2F T 2F T +4F T )/a . (30)
z x x z y z
| − | | |
Comparing it to the special cases (26), (27), and (29),
we recognize that the estimation can be improved by
just dropping the factor 2 in front of T2, i.e. we have
z
to check
F2 T2+T2
F2 =F2+F2+ z +4 x z
Tresca x y 4 a2
2F T 2F T +4F T
+ | z x− x z| | y z| (31) (a) (b) (c)
a
Fig. 1 (a) A typical testsetup with aspect ratio s≃3 and
againstthe thresholdF∗2. Let us call this the simplified N =200.Thesmooth,blueparticlesarethedrivingparticles,
criterion. thestripedbrownonesarethebulkparticles.(b)Thecontact
networkofthesamesystem.(c)Thebrokencontactnetwork
For comparison purposes, we introduce the decou-
after strain ε = 0.00017. Note, that particle displacement
pledcriterion aswell,whereallloadingmodesarechecked wouldnot bevisible.
individually:
F =max F2+F2, F /2,2T /a,2T /a whered¯istheaverageparticlediameter,N isthe num-
Tresca nq x y | z| | x| | z| o ber of particles and s is the desired height width sam-
(32)
ple ratio: s=H/(2R), for which values in the range of
0.5 10 were chosen. The value of f describes the aver-
Whatever criterion (i.e. definition of FTresca) used, age−volume taken up by a particle as 2d¯3/f and thus,
in the case of exceeding the threshold (i.e. F >
Tresca
for a random close packing with a volume fraction of
F∗),thecontactisirreversiblytransformedintoausual,
around 0.65, amounts to f 2.5.
frictionallinearspring-dashpotonewiththesamevisco- ≈
elastic parameters. After allparticleswerepouredinto the cylinder the
samplewasuni-axiallycompressedbytwoplatesinab-
sence of gravity. The compression force employed was
3 Simulation small enough to keep the average deformations of the
friction-less,linearspring-dashpotcontactsbelow10−5
3.1 Simulation setup particleradii.Forthesubsequentfailuretests,thestiff-
ness was reduced by one order of magnitude and the
In this part, we describe in detail the procedure which
pairwise distance was rescaled to avoid pre-stressing.
was used to create standard yield tests on cylindrical
The tests weredone inabsenceofgravityandwith-
shaped specimens using the LAMMPS [10] molecular
out the compressing plates. There are fewer contacts
dynamicssimulationcode,extendedtotreatcakedcon-
fromthebulktoaplatethantoalayerofparticles,thus
tacts according to the model in section 2.
ifplateswereusedfordrivingitwouldbeexpectedthat
In order to create a homogeneous, random close
thecontactsbetweenaplateandthetouchingparticles
packing,first particles were poured into a half cylinder
would fail first. So we chose the last layers of particles
placedhorizontallyinaverticalgravityfield.Thecylin-
(the ones that arecloser than1.1 particle radius to the
derhadthesamediameterasthedesiredtestspecimen
compressing plate) to move as a rigid body and drive
but was necessarily much longer.
the system. This allows the failure zone to be every-
Particleshadanarrowsizedistributionthatislarge
where in the sample. A typical setup and test result is
enough to prohibit crystallizationon the cylinder wall.
shown on Fig. 1.
The number of particles varied from 200 to 10000, the
Two different driving methods were used: (i) the
friction coefficient in the preparation as well as in the
driving layer had infinite mass, (ii) the driving layer
tests was set to 0.
had infinite mass only in the driving direction but fi-
The radius of the cylinder was chosen to be
niteintheperpendiculardirections.Sincenosignificant
N 1/3 differencewasobserved,the resultsreportedinthispa-
R= d¯, (33)
(cid:18)fπs(cid:19) per are shown with driving particles of infinite mass.
Acontact model for theyieldingof caked granular materials 5
The driving velocity was always kept low enough as to whereN isthenumberofcontactsandV isthevolume
c
avoid kinetic effects. of the sample. The latter can be calculated using (33):
After thepreparationthesamplewastested.Atthe
z¯f(2k +3k )
n t
beginning of the tests all existing (i.e. force transmit- E = 12d¯(4k +k ), (36)
n t
ting)contactswere“glued”,allsubsequentcontacts(in-
where z¯is the averagecoordination number which was
cluding broken ones at later times) were handled with
measured to be z¯ 5.06 in all cases. This value is less
the standard Hooke contact law including static fric-
≃
than the isostatic limit for homogeneous packing due
tion. The simulation was run long enough such that
to the boundary effects. Using the values f=2.5 and
the upper and lower part are completely separated as
k =k /2,wegetanestimateforthemacroscopicYoung’s
shown on Fig. 1(c). t n
modulus:
The different tests were realized by giving the driv-
k
ing particles on both sides of the sample initially a ve- n
E =0.81 (37)
locity(orangularvelocity)intherequireddirectionbut d¯
with opposing sign. ThemacroscopicPoissonratiocanalsobe expressedas
function of the contact stiffness [7,2]:
k k
n t
ν = − =0.11 (38)
3.2 Parameters 4kn+kt
Theactualsimulationparameters(innaturalunits)
In order to interpret the results of numerical simula- chosen for the current study were taken to allow reli-
tions we have to take a closer look at the dimensions able and fast simulation of the tests and yet fulfill the
of the parameters. Due to the quasi static nature of requirements above:
the tests (up to the failure point) and the absence of
d¯=1 a=0.1
gravity,quantities liketime andmassareunimportant,
insteadforce and distance units determine the scalesof kn =1256 kt =kn/2
the numerical results. σ∗ =1.88 (39)
Twoparametersofdimensiondistance enter:theav-
erageparticlediameterd¯andthecontactgluediameter FromthosewecancalculatethemacroscopicYoung’s
2a. Since in all simulations the ratio a=0.1d¯was used, modulus(36)andPoissonratio(38)aswellasthelimit
d¯is taken as the unit of length. strain (34):
The stress or force scale is set by the stress thresh- E =1017 , ν =0.11 , ε∗ =9.4 10−5 (40)
×
old σ∗ (cf. (24)) and F∗ = πa2σ∗, respectively. Due to
technicalreasons,17timesthelatterwaschosenasunit
4 Yield tests
of force.
Thus, from now on all quantities will be given in
4.1 Poisson ratio
this natural units (d¯for length and 17F∗ for force).
The normal contact stiffness kn (the tangential one The easiest quantity that we can measure in our sim-
is set to kt kn/2, corresponding to the limit ν = 0, ulations is the Poisson ratio. In pull experiments we
≡
cf.(14))relatesforceandlengthandmustfulfillkn 1 recordedthedisplacementofthedrivingparticlesinthe
≫
to assure the brittle limit. This can also be expressed z directionandtheradiusofthesampleatz=H/2.The
as the failure strain scale latter was done by averaging the outer radius over all
surfaceparticlesthathaveatleasttwocontacts(remove
F∗
ε∗ = knd¯ (34) rmaetntltinisg,vdearyngslminagllp(aorfttichleeso)r.dSeirncoef1th0e−6r)elwateivceandiuspselatchee-
linearized expression for the Poissonratio:
being very small.
∆R
The preparation was done with the aim to achieve ν (41)
asampleashomogeneousandisotropicaspossible.For ≃−∆H
such a random configuration, the macroscopic Young’s We averaged for different samples and system size to
modulus can be calculated analytically [2,7] from the get the result:
stiffness parameters of the contacts:
ν =0.108 0.01 (42)
±
N d¯(2k +3k ) This value is in perfect agreement with the prediction
c n t
E = (35)
of (40).
3V 4k +k
n t
6 L. Brendel etal.
4.2 Stress/strainbehavior The shear deformation is a bit more involved than
the previous ones, since its behavior is fundamentally
Threedifferenttestsweredone:pull,shearandtwist.In differentwhethertheaspectratioislargeorsmallcom-
the tests the driving was realized by giving the driving paredtounity(withacrossoverinbetween).Theshort
particles a constant (angular) velocity. The net force cylinder limit was covered in Sec. 2.1. The long cylin-
(torque), the number of intact contacts, the deforma- der limit is more complicated, but consists of a combi-
tion of the sample was recorded as function of the dis- nation of bending and shearing[13] as well. We define
placement of the driving particles. here the strain as the most straightforward interpola-
The stress/strain evolution of the samples in the tion of the two cases that recovers the short cylinder
homogeneous and isotropic case are described by three resultfor s 0 andthe long cylinder form in the limit
parameters: Young’s modulus2 E, Poisson ratio ν and s : →
∗ →∞
maximalstrainε .Sotheevolutionofalltestshouldbe
X (fπ)1/3X
described by one single curve. Thus the measured dis- ε = =
s,max H(1+ H ) 2s2/3N1/3d¯(1+s/4)
placement and force (torque) data has to be converted 8R
4ε
to a stress/strainrelation. t,max
ε¯ = , (47)
s
Itismoreinformativeifthisconversionisdonewith 3π
the dimensionless quantities of the aspect ratio s and Applying Hooke’s law to the whole strain field yields
the number of particles N since it allows immediately the interpolating stress
the scaling of different tests in a transparent manner.
In general the stress and the strain are not homo- Fx 2ν+2+ 3HR22
σ = (cid:16) (cid:17) =
geneous,so we focus on their two important character- s,max πR2(1+ H )
8R
istics: the maximum and the average. The maximum
stress governsthe point where the first contact breaks, = fs 2/3 (cid:16)2ν+2+ 43s2(cid:17)Fx
the average value indicates the importance of the in- (cid:18)√πN(cid:19) (1+s/4) d¯2
homogeneities of the stress and we compare it to the 4σ
s,max
σ¯ = . (48)
maxima of the stress/strain curve. s
3π
Thepulldeformationishomogeneousthusthemax-
Figure 2 (a) shows the average stress/strain curves
imum is the same as the average:
for different experiments, system size, sample length.
Z (fπ)1/3Z The measured Young’s modulus varies around a mean
ε¯ =ε = = , (43)
p p,max H 2s2/3N1/3d¯ value E 880, which is 15% less than the theoretical
≃
valueof(37).Thereisaslightscatterofthelinearslope
whereZ isthe verticaldisplacementofthedrivingpar-
ofdifferentsystemsbutinspiteofthemoderatesystem
ticles andf is the volumefactorof(33).The stresscan
sizestillonecanobserveaconvincingdatacollapsewith
be taken directly from section 2.1:
different system size, aspect ratio and test type.
2/3
F fs F InFig.2(b)themaximaoftheaveragestress/strain
z z
σ¯ =σ = = (44)
p p,max πR2 (cid:18)√πN(cid:19) d¯2 curves are shown as function of the average strain. It
is evident that the elastic behavior is very close to the
Thetwistdeformationhasaradiusdependence,but
Hooke’s law and can be described by a single Young’s
apart from this is very similar to the previous case.
modulus.Theaveragemacroscopicstressandstrainare
Rα α alwayssmallerthanthescale(40),asmustbeexpected
ε = =
t,max
H 2s due to spatial disorder.
ε¯ =2ε /3, (45) The maximal stress/strain values at the moment
t t,max
of the first broken contact are plotted in Fig. 2 (c).
where α is the angle covered by the driving layer. The
The values are scattered around ε∗ of (40). We note
form of the stress again is identical to the microscopic
two important features: (i) points which are above ε∗
one in 2.1:
are mainly for system size N=200 and for shear and
4T(1+ν) 4f(1+ν)s T twist tests. This is the result of the inhomogeneities of
σ = =
t,max πR3 N d¯3 our system. Due to the preparation procedure the vol-
σ¯ =2σ /3 (46) ume fraction at the lateral surface of the cylinder is
t t,max
larger than in the bulk. This was tested by repeating
2 Justlikeintheprevioussection,weusethesymbolsEand
a few tests with a sample cut out from a bigger one,
νforthemacroscopicquantitiesinsteadoftheonesatcontact
and this feature disappeared. (ii) Maximal stresses for
level as in section 2. The same applies to the deformations
X, Z, to forces/torque Fx, Fz, T and stresses. long cylinders in shear tests are somewhat larger than
Acontact model for theyieldingof caked granular materials 7
0.1 0.16
σ
σ
max
0.14
0.08 (a)
0.12
0.06
0.1
0.04 0.08
0.06
0.02
0.04
4 6 8 10 12 14 16
0
0 2 4 6 8 10 105 ε
max
105 ε
Fig. 3 Thestressandthestrainpointswherethefirstcon-
0.1
σ tactwasbroken.Colorsrepresentexperimenttypepull:blue,
shear: green, twist: red. Symbols indicate the differentbreak
0.08 (b) criterions: Full +, Simplified: ×, Decoupled (cid:3). System size
N=200, 2000, 10000, are from bottom to top respectively
with artificial shiftingof 0.3. The aspect ratio was chosen to
0.06
bes=2
0.04
same initial configuration with the different break cri-
0.02 terions. Fig. 3 shows the summary of the results.
The results show hardly any difference between the
0 two coupled criterions (simplified (31) and full (25),
0 2 4 6 8 10
however the decoupled criterion allows the system to
105 ε evolve without failure beyond the point where the oth-
0.12
σ ers already broke.
max It is remarkable that the simple coupling of the
0.1
(c) different deformation modes leads to indistinguishable
0.08 failure points for large systems. For small systems the
only difference is in the twist deformation.
0.06
The elapsed time during simulations was also mea-
sured. We observed no measurable difference in simu-
0.04
lation duration of tests with decoupled and simplified
0.02 criterions.Simulationswiththefullcriteriononaverage
require 55% more time independently of system size.
0
0 2 4 6 8 10 12 14 Thisadvocatestheuseofthesimplifiedcriterionifcom-
105 ε putationaltimeiscriticalbutthisoverheadisrelatively
max
small in spite of the numerical maximization of (25)
Fig.2 Coloronline.(a)Stressstrainevolution(b)Themax-
which is in general affordable if precision is required.
imumoftheaveragestress,(c)Endofelasticity(firstbroken
contact) in the stress/strain plain. Colors represent experi-
ment type pull: blue, shear: green, twist: red. Symbol fill in-
tensityindicatessystemsizeN:200:half,2000:empty,10000: 4.4 Failure process
full.Symboltypesshowaspectratios:1:⋄,2:(cid:3),3:(cid:13),5:△,
7: ▽, 8: +, 10: × Inthissectionwestudytheprocessofthefailureinthe
tests.Oursimulationsallowthestudyofeachindividual
brokencontact.Thequestionwewanttoanswerhereis:
expected, these are the result of the same skin effect
Whichisthedeformationtypethatisthemostrelevant
during the bending of the sample.
for a give test.
We identify 6 different deformationtypes which are
the core of the simplified break criterion (31):
4.3 Yield criteria
tensile: F2/4
z
In the simulations we used all three criterions intro- shear: F2+F2
x y
duced in Sec. 2.2. We performed test using exactly the tilt: 4T2/a2
x
8 L. Brendel etal.
40 40
(a) (b)
F 2 /F2 F 2 /F2
i T i T
30 30
20 20
10 10
0 0
0 25 50 75 100 0 25 50 75 100
n n
b b
40 40
(c) (d)
F 2 /F2 F 2 /F2
i T i T
30 30
20 20
10 10
0 0
0 25 50 75 100 0 25 50 75 100
n n
b b
Fig.4 Coloronline.Average(sampleandtime)deformationforcesinjustbrokencontactbridgesscaledbythemaximalforce
versus N˜b. (a) Pull test with s=3 (b) Twist test with s=3 (c) Shear test with s=3, (d) Shear test with s=0.5. The symbols
representthe followingdeformation component: tensile: +, shear: ×, tilt:•, torsion: (cid:3), tilt+tensile: (cid:4), torsion+shear: (cid:13). On
therighthand sideof theplotsa typical exampleof force network plotis shownfor thegiven systemafter thetest.
torsion: 4T2/a2 small chains, namely a single particle having two con-
z
tilt+tensile 2F t 2F T /a tactsonetoeachsideofthesample.Ingeneral,sinceall
z x x z
| − |
torsion+shear 4F T /a, (49) “strong” contacts are already broken, its two contacts
y z
| |
are loaded mainly with a combination of torsion and
In order to study the relevance of the above defor-
tilting. A remarkable feature is that tilting gets stabi-
mationtypes,the followingprocedurewasused:Atthe
lized by shear and tensile forces, while torsion is not.
moment of a contact break, the above quantities are
This can be seen by the fact that tilting+tensile defor-
recorded.
mationis almostthe secondmostimportantloadingat
Since we have a random packing, the values show
latestages.Sincewehaveusedarathernarrowcontact
huge variations from contact to contact. An average
diameter a=0.1d¯, contacts break easily under torsional
is necessary to be able to observe the different trends.
load.
We averaged the recorded limit forces for each succes-
sive 5% of broken contacts. So if Nb(t) is the number Twist tests [Fig. 4 (b)] show a surprisingly similar
of broken contacts at time t then the x axes shows picture to pull experiment, except that at early stages
nb=Nb(t)/Nb(T), with T being the time at the end of it is the shear deformation which dominates. The late
the experiment. This representation allows us to aver- stage is also dominated by torsional deformation.
age also for different sample realizations (5-10) so that
The behavior of shear tests are fundamentally dif-
early and late stages of failure can be identified, as
ferentdependingontheaspectratio.Iftheaspectratio
shown in Fig. 4.
is small [s < 1, Fig. 4 (d)] the overall behavior is sim-
The beginning ofthe pull test[Fig.4(a)]is marked
ilar to the twist deformation: shear dominates at early
by contactsstrainedin the tensile direction.This dom-
stages and torsion at late ones.
inance is clearforthe first25%brokencontactsbut re-
mains the most important up to 60% broken contacts. It is common to all of the above cases that the
The end stage (after 75%ofthe contactswere broken), stress/strain behavior of the sample is determined by
is surprisingly characterized by contacts breaking due the contacts that broke in the early stages, moreover
to torsion. Our view is that this is due to breaking of the wholefractureis realizedveryfastinacascadelike
Acontact model for theyieldingof caked granular materials 9
5100 decoupling the loads leads to a rather large error of
N
c 25%, whereas the simplified criterion agreed with the
5000
full one within statistical errors. This applies to our
4900 tests,butinthe caseofthe fullcriterionturningoutto
be necessary in other situations, the found overheadof
4800
55% for the numerical maximization is affordable.
4700 As a benchmark, we have done pull, shear and tor-
siontestswithcylindricalspecimenscreatedfromglued
4600 particles. These exhibited, up to the failure point, an
elastic behavior with the same Young’s modulus for all
4500
2 3 4 5 6 7 8 9 10 experiments and test parameters. This enabled us to
105 ε collapse all different test results onto one single curve.
Thedeviationsofthetheoreticalpredictionforthespec-
Fig. 5 (Color online) The number of contacts as function
oftheaveragestrainfors=3andN=2000,forthetests:pull: imens’ elastic properties were below 15%.
blue,shear:green, twist: red During the tests we payed special attention to the
dominant loading of the contacts that are about to
break. We found that the first half of broken contacts
manner as shown on Fig. 5, as to be expected in the
in pull tests mostly break due to tensile stress,in twist
brittle limit.
testsduetoshearstress,andinsheartestsofveryshort
The only test where a different scenario can be ob-
samples due to shear stress as well. Long samples be-
served is the shear test with long cylinders [Fig. 4 (c)].
have differently in shear tests, though, due to bending
Atearlystagesthetensiledeformationdominates.This
the dominating load gets tensile instead of shear. Fi-
is due to the bending of the sample. The force network
nally, the late stages of the failure is mainly character-
showsthat the contacts startto break attwopoints si-
ized by contacts breaking under torsional load due to
multaneously.These points arethe ones where the ma-
break up of twisted weak chains connecting the almost
terialisstretchedthemost(rightatthemoving“front”
separated parts.
of the driving layer).The stress is very localized so de-
For the sake of a clearer study of this relative im-
spite the damage the sample has still a significant re-
portance of the loading types and the different failure
sistance which is manifested by the long oscillating tail
criteria, we used a constant diameter of the contact
ofthe stress/straincurve[Fig.2(a)]andthe slowstep-
cylinders. In future works, the stochastic nature of the
wise decreasing of the number of contacts.
geometryandstrengthofthe contact,possiblyobeying
a wide distribution[6], will be taken into account.
5 Conclusion
InthispaperweintroducedanewmodelforDEMsim- 6 Acknowledgements
ulations with breakable cemented contacts. Contacts
were considered as small, flat, elastic cylinders which Theauthorsacknowledgethefinancialsupportfromthe
keep the particles in contact until the shear stress in GermanResearchFoundation(DFG grantBR3729/1)
the cementing cylinder reaches a maximal value con- andthankSandiaNationalLaboratories(oftheUSDe-
forming to the Tresca criterion of failure. This allowed partment of Energy) for putting LAMMPS under the
for deriving the elastic and the failure behavior of the GNU Public License.
contact on a common basis. In the limit of very flat
cylinders, the stiffness matrix (14) gets decoupled, jus-
tifying this common a priori modelling Ansatz. References
The problem of finding the maximal shear stress
in an arbitrarily loaded cylinder turned out to possess 1. Antonyuk,S.,Khanal,M.,Tomas,J.,Heinrich,S.,M¨orl,
L.: Impact breakage of spherical granules: Experimen-
analytic solutions too complicated for practical pur-
tal study and dem simulation. Chemical Engineering
poses (but is numerically well accessible). Therefore and Processing: Process Intensification 45(10), 838–856
weconsideredtwoapproximations:Thedecoupledcrite- (2006)
rion coincideswiththeexactmaximizationonlyforsin- 2. Borghi,R.,Bonelli,S.:Towardsaconstitutivelawforthe
unsteady contact stress in granular media. Continuum
gle loading modes,while the simplified criterion agrees
Mech. Thermodyn. 19, 329–345 (2007)
with the latter in more general situations, due to tak-
3. Cundall, P., Strack, O.: A discrete numerical model for
ing into account load combinations. We showed that granular assemblies. Geotechnique29, 47–65 (1979)
10 L. Brendel etal.
4. Delenne,J.Y.,Youssoufi,M.S.E.,Cherblanc,F.:Mechan-
ical behaviour and failure of cohesive granular materi-
als. Int.J.Numer.Anal.Meth.Geomech.28,1577–1594
(2004)
5. Herrmann, H.J., Hovi, J.P., Luding, S. (eds.): The
Physics of Dry Granular Media, vol. 350. Kluwer, Dor-
drecht(1998)
6. Kirsch, R., Brendel, L., T¨oro¨k, J., Bro¨ckel, U.: Measur-
ing tensile, shear and torsional strength of solid bridges
betweenparticlesinthemillimeterregime. Submittedto
thisjournal
7. Kuhl, E., d’Addetta, G.A., Herrmann, H.J., Ramm, E.:
A comparison of discrete granular material models with
continuousmicroplaneformulations. GranularMatter2,
113–122 (2000)
8. Luding, S.: Collisions & contacts between two particles.
In: Herrmann etal. [5], pp.285–304
9. Mitarai, N., Nori, F.: Wet granular materials. Advances
in Physics55(1), 1–45 (2006)
10. Plimpton, S.: Fast parallel algorithms for short-range
molecular dynamics. J. Comp. Phys.117, 1–19 (1995)
11. Po¨schel, T., Luding, S. (eds.): Lecture Notes in Physics:
Granular Gases,vol. 564. Springer(2001)
12. Roux, S.: Quasi-static contacts. In: Herrmann et al. [5],
pp.267–284
13. Sokolnikoff, I.S.: Mathematical Theory of Elasticity.
Krieger, Malabar, Florida (1983)
14. Tresca, H.: M´emoire sur l’´ecoulement des corps solides
soumis `a de fortes pressions. C.R. Acad. Sci. Paris 59,
754 (1864)
15. Wang, Y., Alonso-Marroquin, F.: A finite deformation
method for discrete modeling: particle rotation and pa-
rametercalibration. GranularMatter11,331–343(2009)
