# effect of rolling stiffness on peak friction angle of sample

Asked by Elham Hosseinkhani on 2020-02-02

Hello,
I want to model a sample that is in bellow paper:

Simulation of triaxial response of granular materials by modified DEM
WANG XiaoLiang & LI JiaChun
doi: 10.1007/s11433-014-5605-z
December 2014 Vol. 57 No. 12: 2297–2308

I'm trying to model effects of rolling stiffness on macroscopic mechanical parameters like peak friction angle. I used Ip2_CohFrictMat_CohFrictMat_CohFrictPhys and Law2_ScGeom6D_CohFrictPhys_CohesionMoment and i wrote bellow code.The basic parameters like final friction degree, confining pressure,poisson, psd, young modulus,strain rate and porosity are the same as parameters in that paper. I didn't see any effect and changes on shear strength and peak friction angle after setting rolling stiffness(alphaKr). Actullay i couldn't simulate graphs that are in the paper with different alfaKr (for example for alphaKr =0.01, 1).

1) Is there any problem in my code that i can't simulate graphs in that paper?
2) Should i expect that different alphaKr make a difference in peak friction angle of sample?

## encoding: utf-8
import matplotlib; matplotlib.rc('axes',grid=True)
import matplotlib.pyplot as plt
import pylab
import numpy as np
from numpy import *
import math
import sys
import time
from math import *

## Define Parameters
num_spheres=10000
compFricDegree=20
finalFricDegree=30
confiningS=-100000 # [Pa]
RhoW=1000
graindensity=2600
poissonRatio=0.4
youngModulus=270e6 # [Pa]

AxialstrainMatrix1=[]
DeviatoricStressMatrix1=[]
VolumetricStrainMatrix1=[]
PrincipalStressRatioMatrix1=[]
qPMatrix1=[]
AxialstrainMatrix2=[]
DeviatoricStressMatrix2=[]
VolumetricStrainMatrix2=[]
PrincipalStressRatioMatrix2=[]
qPMatrix2=[]

alphaKrMatrix=[0.01,1]
for i in range(0,len(alphaKrMatrix),1):
psdSizes=[0.107,0.12,0.132,0.145,0.158,0.171,0.183,0.196,0.209,0.221,0.234] #(mm)
psdCumm=[0.001,0.029,0.069,0.121,0.187,0.268,0.367,0.49,0.63,0.8,1.0] #cumulative
psdSizesArray=np.array(psdSizes)
psdSizesMeter=psdSizesArray*0.001 #Convert the size of particles to meter
sp=pack.SpherePack()
mn,mx=Vector3(0,0,0),Vector3(0.0045,0.0075,0.0045) #initial box size
sp.makeCloud(minCorner=mn,maxCorner=mx,num=num_spheres,psdSizes=psdSizesMeter,psdCumm=psdCumm,distributeMass=True,seed=True)
sp.psd(bins=200,mass=True)

O.materials.append(FrictMat(young=youngModulus,poisson=poissonRatio,frictionAngle=0,density=0,label='frictionless'))

walls=aabbWalls((mn,mx),thickness=0,material='frictionless')
wallIds=O.bodies.append(walls)

triax=TriaxialStressController(
internalCompaction=False,
goal1=confiningS,
goal2=confiningS,
goal3=confiningS,
label="triax"
)
O.engines=[
ForceResetter(),
InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
InteractionLoop(
[Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
[Ip2_FrictMat_FrictMat_FrictPhys(), Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(label="cohesiveIp")],
[Law2_ScGeom_FrictPhys_CundallStrack(traceEnergy=True),Law2_ScGeom6D_CohFrictPhys_CohesionMoment(useIncrementalForm=True,always_use_moment_law=False,label='cohesiveLaw')]
),
GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
triax,
NewtonIntegrator(damping=0.3,label="newton"),
]

O.dt=utils.PWaveTimeStep()
O.dynDt=False

while 1:
O.run(1000,True)
unb=unbalancedForce()
e=(triax.porosity)/(1-triax.porosity)
print ('unb:',np.round(unb,2),'e: ',np.round(e,3),'meanstress:', np.round(-triax.meanStress,3),' SigmaY: ',np.round(-triax.stress(3),3),' SigmaX: ',np.round(-triax.stress(1),3),' SigmaZ: ',np.round(-triax.stress(5),3))
if unb<0.1 and abs(confiningS-triax.meanStress)/abs(confiningS)<0.01 and e<0.634:
break

##########################
## Start test ###
##########################
# triaxial conditions
triax.internalCompaction=False
triax.goal1=triax.goal3=confiningS

triax.wall_bottom_activated=False
cohesiveLaw.always_use_moment_law=True
triax.goal2=-5
if i==0:
while 1:
O.run(1000,True)
AxialstrainMatrix1.append((-triax.strain)*100)
qPMatrix1.append((-triax.stress(triax.wall_top_id)-(-triax.stress(triax.wall_right_id)-triax.stress(triax.wall_front_id))/2.0)/(-triax.meanStress))
print ('eyy=',triax.strain,'Sxx=', triax.stress(0),'Syy=',triax.stress(3), 'Szz=',triax.stress(4),'e=', (triax.porosity)/(1-triax.porosity))
if triax.strain <-0.3:
break
if i==1:
while 1:
O.run(1000,True)
AxialstrainMatrix2.append((-triax.strain)*100)
qPMatrix2.append((-triax.stress(triax.wall_top_id)-(-triax.stress(triax.wall_right_id)-triax.stress(triax.wall_front_id))/2.0)/(-triax.meanStress))
print ('eyy=',triax.strain,'Sxx=', triax.stress(0),'Syy=',triax.stress(3), 'Szz=',triax.stress(4),'e=', (triax.porosity)/(1-triax.porosity))
if triax.strain <-0.3:
break
O.reset()

plt.figure()
plt.subplot(224)
plt.xlabel('Axial Strain (%)')
plt.ylabel('q/P')
plt.grid(True)
plt.plot(AxialstrainMatrix1,qPMatrix1,linestyle='-',color='b',linewidth=1)
plt.plot(AxialstrainMatrix2,qPMatrix2,linestyle='-',color='r',linewidth=1)
plt.show()

## Question information

Language:
English Edit question
Status:
Solved
For:
Assignee:
No assignee Edit question
Solved by:
Elham Hosseinkhani
Solved:
2020-02-18
Last query:
2020-02-18
2020-02-12
 Robert Caulk (rcaulk) said on 2020-02-02: #1

>>I want to model a sample that is in bellow paper:
>> Is there any problem in my code that i can't simulate graphs in that paper?

Have you tried contacting the corresponding authors for assistance with this ? It is their responsibility to make sure their plots are reproducible. They may be using different contact laws than the ones you've selected.

 LI JiaChun, email: <email address hidden>; WANG XiaoLiang, email: wangxiaoLiang52086@126.com

 Elham Hosseinkhani (ehosseinkhani) said on 2020-02-02: #2

Thanks Robert Caulk for your help,

>> Have you tried contacting the corresponding authors for assistance with this?

Actually, based on relations of contact laws that are in the paper, i found that i need these two contact laws for simulating rolling stiffness in my sample. But i'll send an email for paper authors and ask for more help.

 Luc Sibille (luc-sibille) said on 2020-02-02: #3

Dear Elham,

Please note that the two papers  and  (below) are based on the
contact law you want to use.
It has been shown in paper  that the peak friction angle does not
depend on alphaKr when alphaKr is sufficiently high. Maybe you are in
this latter case...?

Regards,
Luc

 Discrete numerical modeling of loose soil with spherical particles
and interparticle rolling friction
RA Hosn, L Sibille, N Benahmed, B Chareyre
Granular matter 19 (1), 4
http://hal.univ-grenoble-alpes.fr/hal-01951842/document

 Quantitative prediction of discrete element models on complex
Luc Sibille, Pascal Villard, Félix Darve, Rodaina Aboul Hosn
International Journal for Numerical and Analytical Methods in Geomechanics
43(5): 858-887
https://hal.archives-ouvertes.fr/hal-02059431/document

--
Luc Sibille
Université Grenoble Alpes / IUT1 de Grenoble
Laboratoire 3SR: Sols, Solides, Structures, Risques

Tel lab.: +33 (0)4 76 82 63 48
Tel IUT: +33 (0)4 76 82 53 36

 Elham Hosseinkhani (ehosseinkhani) said on 2020-02-03: #4

Hello Luc Sibille, and thanks for your comment.

I've studied paper  before and I've observed that point you said (peak friction angle does not depend on alphaKr when alphaKr is sufficiently high). But my problem is about the amount of peak friction angle before considering rolling stiffness (I mean without considering CohesionMomment) and after that. when I use just FrictPhys_CundallStrack as the contact law, with the Inter-Particle friction angle of 30 degree, I can't get peak friction angle of experimental data that is about 40 degree. Because of this problem, I tried to model a sample with more stiffness until getting more peak friction angle. After using CohesionMoment and setting alpaKr (for example alphaKr=0.5 or 1), I couldn't see any changes.
I should say that I know about using more Inter-Particle friction angle( higher than 30 Degree like 50, 60,...) and getting more peak friction angle in the macroscopic response. I've tried this dependency an observed exact results. But I want to see effects of rolling parameters like stiffness and strength in macroscopic parameters.
Thanks Luc Sibille for your suggestion for more studying. I also will study paper .

>> Discrete numerical modeling of loose soil with spherical particles
and interparticle rolling friction
RA Hosn, L Sibille, N Benahmed, B Chareyre

>> Quantitative prediction of discrete element models on complex
Luc Sibille, Pascal Villard, Félix Darve, Rodaina Aboul Hosn

 Luc Sibille (luc-sibille) said on 2020-02-03: #5

> But my problem is about the amount of peak friction angle before considering rolling stiffness (I mean without considering CohesionMomment) and after that. when I use just FrictPhys_CundallStrack as the contact law, with the Inter-Particle friction angle of 30 degree, I can't get peak friction angle of experimental data that is about 40 degree.
This is a classical result: with spherical particles and without rolling
resistance, internal peak friction angle relatively high (typically 40
deg) cannot be reproduced.
> Because of this problem, I tried to model a sample with more stiffness until getting more peak friction angle. After using CohesionMoment and setting alpaKr (for example alphaKr=0.5 or 1), I couldn't see any changes.
Again, stiffnesses, as normal, tangential or rolling stiffnesses, if
their are chosen sufficient high should not affect the failure
properties as the peak friction angle. So it may be normal that peak
friction angle is not affected by such stiffnesses. Then peak friction
angle is more strongly affected by the contact friction angle and
plastic rolling moment coefficient (or rolling friction coefficient)
etaRoll (and not alphaKr).

Besides, be sure that the rolling resistance is activated. It is not
because you use Law2_ScGeom6D_CohFrictPhys_CohesionMoment that rolling
resistance is necessarily activated: rolling resistance can be enable or
disable depending on some arguments in
Law2_ScGeom6D_CohFrictPhys_CohesionMoment and/or CohFrictMat

Luc

--
Luc Sibille
Université Grenoble Alpes / IUT1 de Grenoble
Laboratoire 3SR: Sols, Solides, Structures, Risques

Tel lab.: +33 (0)4 76 82 63 48
Tel IUT: +33 (0)4 76 82 53 36

 Elham Hosseinkhani (ehosseinkhani) said on 2020-02-09: #6

Hello,
Unfortunately,I couldn't get the results that I expected based on the papers. I changed my script and modify the script in YADE examples. I added CohesionMomentLaw and necessary parameters like alphaKr and etaRoll. I get PrincipalStressRatio about 2.4 with this parameters in this script. But I expect more PrincipalStressRatio after using CohesionMomentLaw.

Could you please run this script and check PrincipalStressRatio?
Is there any problem that i couldn't get changes in PrincipalStressRatio and peak friction angle?

# -*- coding: utf-8 -*-

num_spheres=10000,
young = 500e6,
poisson = 0.4,
compFricDegree = 30,
alphaKr = 1.25,
etaRoll = 0.55,
finalFricDegree = 30,
targetPorosity = 0.4,
confining =-300000,
unknownOk=True
)
num_spheres=table.num_spheres
young = table.young
poisson = table.poisson
compFricDegree = table.compFricDegree
alphaKr = table.alphaKr
etaRoll = table.etaRoll
finalFricDegree = table.finalFricDegree
targetPorosity = table.targetPorosity
confining = table.confining
rate=-0.02
damp=0.2

key='_P='+str(confining)
stabilityThreshold=0.01
mn,mx=Vector3(0,0,0),Vector3(1,2,1)

O.materials.append(FrictMat(young=young,poisson=poisson,frictionAngle=0,density=0,label='walls'))

walls=aabbWalls([mn,mx],thickness=0,material='walls')
wallIds=O.bodies.append(walls)

sp=pack.SpherePack()

clumps=False
if clumps:

volume = (mx-mn)*(mx-mn)*(mx-mn)
sp.makeClumpCloud(mn,mx,[c1],periodic=False)
sp.toSimulation(material='spheres')
O.bodies.updateClumpProperties()
else:
sp.makeCloud(mn,mx,-1,0.3333,num_spheres,False, 0.95,seed=1) #"seed" make the "random" generation always the same

triax=TriaxialStressController(
maxMultiplier=1.+2e5/young,
finalMaxMultiplier=1.+2e4/young,
thickness = 0,
internalCompaction=True,
)

newton=NewtonIntegrator(damping=damp)

O.engines=[
ForceResetter(),
InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
InteractionLoop(
[Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
[Ip2_FrictMat_FrictMat_FrictPhys(),Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(label="cohesiveIp")],
[Law2_ScGeom_FrictPhys_CundallStrack(),Law2_ScGeom6D_CohFrictPhys_CohesionMoment(
useIncrementalForm=True,
always_use_moment_law=False,
label='cohesiveLaw')]
),
GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
triax,
TriaxialStateRecorder(iterPeriod=100,file='WallStresses'),
newton
]

Gl1_Sphere.stripes=0

triax.goal1=triax.goal2=triax.goal3=confining

while 1:
O.run(1000, True)
unb=unbalancedForce()
print ('unbalanced force:',unb,' mean stress: ',triax.meanStress)
if unb<stabilityThreshold and abs(confining-triax.meanStress)/abs(confining)<0.01:
break

import sys
while triax.porosity>targetPorosity:
compFricDegree = 0.95*compFricDegree
print ("\r Friction: ",compFricDegree," porosity:",triax.porosity)
sys.stdout.flush()
O.run(500,1)

#############################
#############################

triax.internalCompaction=False
triax.goal2=rate
triax.goal1=confining
triax.goal3=confining
O.engines.lawDispatcher.functors.always_use_moment_law = True
#cohesiveLaw.always_use_moment_law=True
triax.wall_bottom_activated=False
newton.damping=0.1
O.saveTmp()

def history():
ev=-triax.strain-triax.strain-triax.strain,
s11=-triax.stress(triax.wall_right_id),
s22=-triax.stress(triax.wall_top_id),
s33=-triax.stress(triax.wall_front_id),
PrincipalStressRatio=abs((-triax.stress(3))/((-triax.stress(triax.wall_right_id)-triax.stress(triax.wall_front_id))/2)),
q=-triax.stress(triax.wall_top_id)-(-triax.stress(triax.wall_right_id)-triax.stress(triax.wall_front_id))/2.0,
i=O.iter)

if 1:
O.engines=O.engines[0:5]+[PyRunner(iterPeriod=20,command='history()',label='recorder')]+O.engines[5:7]
else:
O.engines=PyRunner(iterPeriod=20,command='history()',label='recorder')

O.run(1000,True)
plot.plots={'e22':('PrincipalStressRatio',None,'ev')}
plot.plot()

 Luc Sibille (luc-sibille) said on 2020-02-12: #7

Did you check if there is actually a rolling moment not nil at the
contacts between spheres?
Luc

Le 09/02/2020 à 10:23, Elham Hosseinkhani a écrit :
> Question #688439 on Yade changed:
>
> Elham Hosseinkhani posted a new comment:
> Hello,
> Unfortunately,I couldn't get the results that I expected based on the papers. I changed my script and modify the script in YADE examples. I added CohesionMomentLaw and necessary parameters like alphaKr and etaRoll. I get PrincipalStressRatio about 2.4 with this parameters in this script. But I expect more PrincipalStressRatio after using CohesionMomentLaw.
>
> Could you please run this script and check PrincipalStressRatio?
> Is there any problem that i couldn't get changes in PrincipalStressRatio and peak friction angle?
>
> # -*- coding: utf-8 -*-
>
>
> num_spheres=10000,
> young = 500e6,
> poisson = 0.4,
> compFricDegree = 30,
> alphaKr = 1.25,
> etaRoll = 0.55,
> finalFricDegree = 30,
> targetPorosity = 0.4,
> confining =-300000,
> unknownOk=True
> )
> num_spheres=table.num_spheres
> young = table.young
> poisson = table.poisson
> compFricDegree = table.compFricDegree
> alphaKr = table.alphaKr
> etaRoll = table.etaRoll
> finalFricDegree = table.finalFricDegree
> targetPorosity = table.targetPorosity
> confining = table.confining
> rate=-0.02
> damp=0.2
>
> key='_P='+str(confining)
> stabilityThreshold=0.01
> mn,mx=Vector3(0,0,0),Vector3(1,2,1)
>
> O.materials.append(FrictMat(young=young,poisson=poisson,frictionAngle=0,density=0,label='walls'))
>
> walls=aabbWalls([mn,mx],thickness=0,material='walls')
> wallIds=O.bodies.append(walls)
>
> sp=pack.SpherePack()
>
> clumps=False
> if clumps:
>
> volume = (mx-mn)*(mx-mn)*(mx-mn)
> sp.makeClumpCloud(mn,mx,[c1],periodic=False)
> sp.toSimulation(material='spheres')
> O.bodies.updateClumpProperties()
> else:
> sp.makeCloud(mn,mx,-1,0.3333,num_spheres,False, 0.95,seed=1) #"seed" make the "random" generation always the same
>
> triax=TriaxialStressController(
> maxMultiplier=1.+2e5/young,
> finalMaxMultiplier=1.+2e4/young,
> thickness = 0,
> internalCompaction=True,
> )
>
> newton=NewtonIntegrator(damping=damp)
>
> O.engines=[
> ForceResetter(),
> InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
> InteractionLoop(
> [Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
> [Ip2_FrictMat_FrictMat_FrictPhys(),Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(label="cohesiveIp")],
> [Law2_ScGeom_FrictPhys_CundallStrack(),Law2_ScGeom6D_CohFrictPhys_CohesionMoment(
> useIncrementalForm=True,
> always_use_moment_law=False,
> label='cohesiveLaw')]
> ),
> GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
> triax,
> TriaxialStateRecorder(iterPeriod=100,file='WallStresses'),
> newton
> ]
>
> Gl1_Sphere.stripes=0
>
> triax.goal1=triax.goal2=triax.goal3=confining
>
> while 1:
> O.run(1000, True)
> unb=unbalancedForce()
> print ('unbalanced force:',unb,' mean stress: ',triax.meanStress)
> if unb<stabilityThreshold and abs(confining-triax.meanStress)/abs(confining)<0.01:
> break
>
> import sys
> while triax.porosity>targetPorosity:
> compFricDegree = 0.95*compFricDegree
> print ("\r Friction: ",compFricDegree," porosity:",triax.porosity)
> sys.stdout.flush()
> O.run(500,1)
>
> #############################
> #############################
>
> triax.internalCompaction=False
> triax.goal2=rate
> triax.goal1=confining
> triax.goal3=confining
> O.engines.lawDispatcher.functors.always_use_moment_law = True
> #cohesiveLaw.always_use_moment_law=True
> triax.wall_bottom_activated=False
> newton.damping=0.1
> O.saveTmp()
>
>
> def history():
> ev=-triax.strain-triax.strain-triax.strain,
> s11=-triax.stress(triax.wall_right_id),
> s22=-triax.stress(triax.wall_top_id),
> s33=-triax.stress(triax.wall_front_id),
> PrincipalStressRatio=abs((-triax.stress(3))/((-triax.stress(triax.wall_right_id)-triax.stress(triax.wall_front_id))/2)),
> q=-triax.stress(triax.wall_top_id)-(-triax.stress(triax.wall_right_id)-triax.stress(triax.wall_front_id))/2.0,
> i=O.iter)
>
> if 1:
> O.engines=O.engines[0:5]+[PyRunner(iterPeriod=20,command='history()',label='recorder')]+O.engines[5:7]
> else:
> O.engines=PyRunner(iterPeriod=20,command='history()',label='recorder')
>
> O.run(1000,True)
> plot.plots={'e22':('PrincipalStressRatio',None,'ev')}
> plot.plot()
>

--
Luc Sibille
Université Grenoble Alpes / IUT1 de Grenoble
Laboratoire 3SR: Sols, Solides, Structures, Risques

Tel lab.: +33 (0)4 76 82 63 48
Tel IUT: +33 (0)4 76 82 53 36

 Elham Hosseinkhani (ehosseinkhani) said on 2020-02-18: #8

Hello,
I found the problem of my code and now it works properly. I could reproduce some results in some papers.
I had to change the way of making pack. I used uniform distribution and activate this option; "always_use_moment_law=True", in contact law during making pack and deviatoric loading. I thank from Luc Sibille for his help.