Simulation compared to real tests

Asked by Amedeo Galletti on 2020-11-20

Goodmoring,

using this script [1] I've made some simulation with different value of SN parameter and PI (25000, 50000 an 75000 for SN and respectively 2500, 5000 and 7500 for PI). By comparing those result with real shear test performed on a Bristol Sand ( confining stress was set on 25, 50 and 75 kPa) I've observed that the simulated stress-strain graph is exactly the half of the stress-strain curve at each confining stress. I've tried to double SN and PI parameters on the simulation, and the results o a stress-strain curve match the real stress-strain curve at each confining stress. How is that possible? Is There something I don't understand that makes this possible, like a "trick" on the script?

Thanks for your help.

[1]
###

from yade import pack,plot,export
import math

sp=pack.SpherePack()
O.periodic=True

# dimensions of sample (fixed by particle size such as L/D~15)
RADIUS=0.001
length=0.06
height=length/3
width=length
thickness=RADIUS

# friction angles
compFRIC=15. # during compaction (controls porosity)
FRIC=40. # during shear

# boundary conditions
PI=5e3 # sample preparation: pressure applied for the isotropic compaction
SN=50e3 # normal stress
RATE=0.1 # shear velocity (top plate)

# simulation control
DAMPSHEAR=0
ITER=8e4
VTK=20
OUT='TEST'

####

O.cell.hSize=Matrix3(length,0,0,0,3*height,0,0,0,width)

O.materials.append(CohFrictMat(isCohesive=True,density=2500,young=1e8,poisson=0.5,frictionAngle=radians(0.),normalCohesion=1e100,shearCohesion=1e100,label='boxMat'))
O.materials.append(CohFrictMat(isCohesive=True,density=1600,young=230e6,poisson=0.3,frictionAngle=radians(compFRIC),normalCohesion=0,shearCohesion=0,label='sphereMat'))

upBox = utils.box(center=(0,2*height+thickness/2.0,0),orientation=Quaternion(1,0,0,0),extents=(2*length,thickness/2.,2*width),fixed=1,wire=False,color=(1,0,0),material='boxMat')
lowBox = utils.box(center=(0,height-thickness/2.0,0),orientation=Quaternion(1,0,0,0),extents=(2*length,thickness/2.,2*width),fixed=1,wire=False,color=(1,0,0),material='boxMat')
O.bodies.append([upBox,lowBox])

sp.makeCloud((0,height+1.5*RADIUS,0),(length,2*height-1.5*RADIUS,width),rMean=RADIUS,rRelFuzz=0.001,periodic=True)
O.bodies.append([utils.sphere(s[0],s[1],color=(0,0,1),material='sphereMat') for s in sp])

effCellVol=(O.bodies[0].state.pos[1]-O.bodies[1].state.pos[1])*O.cell.hSize[0,0]*O.cell.hSize[2,2]
volRatio=(O.cell.hSize[0,0]*O.cell.hSize[1,1]*O.cell.hSize[2,2])/effCellVol

#print 'volRatio=',volRatio

O.engines=[
 ForceResetter()
 ,InsertionSortCollider([Bo1_Box_Aabb(),Bo1_Sphere_Aabb()],verletDist=-0.1,allowBiggerThanPeriod=True)
 ,InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom6D(),Ig2_Box_Sphere_ScGeom6D()],
  [Ip2_CohFrictMat_CohFrictMat_CohFrictPhys()],
  [Law2_ScGeom6D_CohFrictPhys_CohesionMoment()]
 )
 ,GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8,defaultDt=utils.PWaveTimeStep())
 ,PeriTriaxController(dynCell=True,maxUnbalanced=1e-3,relStressTol=1e-4,stressMask=7,goal=(-PI/volRatio,-PI/volRatio,-PI/volRatio),globUpdate=1,maxStrainRate=(1,1,1),doneHook='triaxDone()',label='triax')
 ,NewtonIntegrator(damping=0.3,label='newton')
 ,PyRunner(command='dataRecorder()',iterPeriod=100,label='recData',dead=False)
 ,VTKRecorder(fileName=OUT+'.',iterPeriod=1,skipNondynamic=1,recorders=['spheres','colors','velocity','bstresses','stress','boxes','intr','boxes'],label='saveSolid',dead=True)
]

def dataRecorder():
 h=vol=vol_s=nb_s=0.
 h=O.bodies[0].state.pos[1]-O.bodies[1].state.pos[1]
 vol=h*O.cell.hSize[0,0]*O.cell.hSize[2,2]
 contactStress=getStress(vol)
 for o in O.bodies:
   if isinstance(o.shape,Sphere) and o.shape.color[0]!=1:
    nb_s+=1
    vol_s += 4.*pi/3.*(o.shape.radius)**3
 n = 1-vol_s/vol
 nbFrictCont=0.
 for i in O.interactions:
  if i.isReal and i.phys.cohesionBroken:
   nbFrictCont+=1
 plot.addData(
  iter=O.iter
  ,contactStress01=((contactStress[0,1])/1000),contactStress11=(contactStress[1,1])
  ,xW=(O.bodies[0].state.pos[0]*1000)
  ,height=h
  ,volume=vol
  ,porosity=n
 )

def triaxDone():
 global phase
 volRatio=(O.cell.hSize[0,0]*O.cell.hSize[1,1]*O.cell.hSize[2,2])/((O.bodies[0].state.pos[1]-O.bodies[1].state.pos[1])*O.cell.hSize[0,0]*O.cell.hSize[2,2])
 h=O.bodies[0].state.pos[1]-O.bodies[1].state.pos[1]
 vol=h*O.cell.hSize[0,0]*O.cell.hSize[2,2]
 contactStress=getStress(vol)
 vol_s=Rmean=Rmax=nbSph=0
 Rmin=1e6
 for o in O.bodies:
  if isinstance(o.shape,Sphere):
   nbSph+=1
   Rmean+=o.shape.radius
   if o.shape.radius>Rmax: Rmax=o.shape.radius
   if o.shape.radius<Rmin: Rmin=o.shape.radius
   vol_s += 4.*pi/3.*(o.shape.radius)**3
 Rmean=Rmean/nbSph
 n = 1-vol_s/vol
 print 'DONE! iter=',O.iter,'| sample generated: nb spheres',nbSph,', Rmean=',Rmean,', Rratio=',Rmax/Rmin,', porosity=',n
 print 'Changing contact properties now'
 tt=TriaxialCompressionEngine()
 tt.setContactProperties(FRIC)
 triax.dead=True
 O.pause()

#### Initialization
print 'SAMPLE PREPARATION!'

recData.dead=True # uncomment to record what is happening during stress initialization
O.run(600000,1)
saveSolid.dead=False
#saveSolid.fileName=OUT+'_isoConfined.'
O.step()
saveSolid.dead=True
O.save(OUT+'_initialState.yade')

print 'Normal stress (platen) = ',O.forces.f(0)[1]/(O.cell.hSize[0,0]*O.cell.hSize[2,2])
print 'Normal stress (contacts) = ',getStress((O.bodies[0].state.pos[1]-O.bodies[1].state.pos[1])*O.cell.hSize[0,0]*O.cell.hSize[2,2])[1,1]

#### Applying normal stress
print 'NORMAL LOADING! iter=',O.iter

recData.dead=True
stage=0
stiff=fnPlaten=currentSN=0.
def servo():
 global stage,stiff,fnPlaten,currentSN
 if stage==0:
  currentSN=O.forces.f(0)[1]/(O.cell.hSize[0,0]*O.cell.hSize[2,2])
  unbF=unbalancedForce()
  #print 'SN=',SN,'| current SN = ',currentSN,' | unbF=',unbF
  boundaryVel=copysign(min(0.1,abs(0.5*(currentSN-SN))),currentSN-SN)
  O.bodies[0].state.vel[1]=boundaryVel
  if ( (abs(currentSN-SN)/SN)<0.001 and unbF<0.001 ):
   stage+=1
   fnPlaten=O.forces.f(0)[1]
   print 'Normal stress =',currentSN,' | unbF=',unbF
   ## the following computes the stiffness of the plate (used for stress control of the top plate)
   for i in O.interactions.withBody(O.bodies[0].id):
    stiff+=i.phys.kn
   print 'DONE! iter=',O.iter
   O.pause()
 if stage==1:
  fnDesired=SN*(O.cell.hSize[0,0]*O.cell.hSize[2,2])
  #boundaryVel=copysign(abs(0.5*(O.forces.f(0)[1]-fnDesired)/stiff/O.dt),O.forces.f(0)[1]-fnDesired)
  boundaryVel=copysign(min(RATE,abs(0.333*(O.forces.f(0)[1]-fnDesired)/stiff/O.dt)),O.forces.f(0)[1]-fnDesired)
  O.bodies[0].state.vel[1]=boundaryVel

O.engines = O.engines[:4]+[PyRunner(command='servo()',iterPeriod=1,label='servo')]+O.engines[4:]

O.run(35000,1)
print 'STABILIZING! iter=',O.iter
O.run(10000,1)
# coloring particles to see vertical stripes in material
dxi=4*(2*RADIUS) # can be a function of cell length dxi=O.cell.hSize[0,0]/5.
n=int(length/dxi)
xmin=1e6
for i in range(0,n):
 for o in O.bodies:
  if o.id>1:
   if o.state.pos[0]<xmin: xmin=o.state.pos[0]
   if (o.state.pos[0]>=i*dxi) and (o.state.pos[0]<((i+0.5)*dxi)):
    o.shape.color[1]=1
for o in O.bodies:
 if (o.state.pos[0]>(xmin+0.5*O.cell.hSize[0,0]-RADIUS)) and (o.state.pos[0]<(xmin+0.5*O.cell.hSize[0,0]+3*RADIUS)):
  o.shape.color[2]=0

saveSolid.dead=False
#saveSolid.fileName=OUT+'_normalLoaded.'
O.step()
saveSolid.dead=True
O.save(OUT+'_normallyLoaded.yade')
if recData.dead==False: plot.saveDataTxt(OUT)

print 'Normal stress (platen) = ',O.forces.f(0)[1]/(O.cell.hSize[0,0]*O.cell.hSize[2,2])
print 'Normal stress (contacts) = ',getStress((O.bodies[0].state.pos[1]-O.bodies[1].state.pos[1])*O.cell.hSize[0,0]*O.cell.hSize[2,2])[1,1]

#### preparing for shearing
Gl1_Sphere.stripes=1
print 'Gluing spheres to boundary walls'
gluedSpheres=[]

## gluing particles in contact with the walls
for i in O.interactions:
  if i.isReal:
   if isinstance(O.bodies[i.id1].shape,Box):
    O.bodies[i.id2].shape.color[0]=1
    gluedSpheres += [O.bodies[i.id2]]
   if isinstance(O.bodies[i.id2].shape,Box):
    O.bodies[i.id1].shape.color[0]=1
    gluedSpheres += [O.bodies[i.id1]]

for i in O.interactions:
  if i.isReal and ( O.bodies[i.id1].shape.color[0]==1 and O.bodies[i.id2].shape.color[0]==1 ):
   O.bodies[i.id1].mat.normalCohesion=O.bodies[i.id1].mat.normalCohesion
   O.bodies[i.id2].mat.normalCohesion=O.bodies[i.id1].mat.normalCohesion
   O.bodies[i.id1].mat.shearCohesion=O.bodies[i.id1].mat.shearCohesion
   O.bodies[i.id2].mat.shearCohesion=O.bodies[i.id1].mat.shearCohesion
   i.phys.initCohesion=True

print 'nb of glued spheres=',len(gluedSpheres)

#### shearing
print 'SHEARING! iter=',O.iter

saveSolid.dead=False
O.step()
saveSolid.dead=True
O.save(OUT+'_shearInit.yade')

recData.dead=False
newton.damping=DAMPSHEAR
saveSolid.dead=False
saveSolid.iterPeriod=int(ITER/VTK)
shearVel=0
iterShear=O.iter
while 1:
 O.run(int(10),1)
 if shearVel<RATE:
  shearVel+=(RATE/100.)
 # only top wall moves
 O.bodies[0].state.vel[0]=shearVel
 ## top and bottom walls move
 #O.bodies[0].state.vel[0]=0.5*shearVel
 #O.bodies[1].state.vel[0]=-0.5*shearVel
 if ( O.iter >= (int(iterShear+ITER)) ):
  print 'iter=',O.iter,' -> FINISHED!'
  plot.saveDataTxt(OUT)
  O.save(OUT+'_sheared.yade')
  sys.exit(0)

Question information

Language:
English Edit question
Status:
Answered
For:
Yade Edit question
Assignee:
No assignee Edit question
Last query:
2020-11-20
Last reply:
2020-11-22
Luc Scholtès (luc) said : #1

Hello,

How did you define the interparticle properties in your simulation (interparticle stiffnesses, interparticle friction)?

Did you calibrate the model to make sure that the macroscopic emergent behavior is the behavior you expect?

My guess is that the emergent/macroscopic friction angle of your DEM model is exactly half the friction angle of the material you want to model.

Because of the spherical simplification (and other things), defining an interparticle friction angle of 40° will not result in a macroscopic/emergent friction angle of 40°. Generally, the emergent macroscopic friction angle is smaller than the interparticle friction angle.

Sorry if you already knew that but I just wanted to be sure before going any further.

Luc

Amedeo Galletti (amedeog) said : #2

Hello,

Thanks Luc for your quick answer.

The properties in the script are the ones that I used in the simulation.

Yes, I've calibrated the model.

Yes, I'm aware of the spherical simplification.

Amedeo Galletti (amedeog) said : #3

Hello,

Thanks Luc for your quick answer.

The properties in the script are the ones that I used in the simulation.

Yes, I've calibrated the model.

Yes, I'm aware of the spherical simplification.

How can I check if the emergent/macroscopic friction angle of my DEM model is exactly half of my real model?

Sorry for the double replay, I misclicked.
Regards,
Amedeo

Luc Scholtès (luc) said : #4

The fact that you asked about how to measure the emergent properties of your DEM model confuses me a bit since this is part of the calibration process which consists in finding the set of interparticle properties (stiffnesses, friction) that enables to obtain the right macroscopic properties. Usually, this is done by running either triaxial or direct shear tests to extract the emergent properties for a given set of interparticle properties.

Just to be sure: how did you calibrate your model?

Luc

Amedeo Galletti (amedeog) said : #5

I calibrated my model by changing the various parameters one by one to see what would happened to the results of the simulation. Then I choose the parameters that matches my soil. My soil is a Bristol Sand. I've the results of real direct shear test on that sand. The real friction angle is 38°.

Sorry for my poor English.
Thanks for you patience and your help.

Amedeo

Luc Scholtès (luc) said : #6

OK... I am still confused because you say that you chose the parameters that matches your soil but then you say that there is a discrepancy when you compare the stress strain curves...

If I understand correctly, your real soil has an internal friction angle of 38°. To obtain this value, you performed several direct shear tests at different normal stresses and then, based on the Mohr-Coulomb criterion, you used a linear regression to fit your data (peak or residual shear stresses and associated normal stresses).

My question is: what is the internal/emergent/macroscopic friction angle of your DEM model/numerical material?

I see that you defined a 40° in the following lines of your script:

# friction angles
compFRIC=15. # during compaction (controls porosity)
FRIC=40. # during shear

These lines defines the interparticle friction angle, the value that is used in CohFrictMat.

This value (40°) does not define the internal/emergent/macroscopic friction angle of your DEM model/material.

To estimate the internal/emergent/macroscopic friction angle of your DEM model/material, you need to do exactly what you have done with your soil: for a given value of the interparticle friction angle, perform 2 (at least) simulations with different normal stresses and then fit the results with the MC criterion. You will see that there is no direct correspondence between the interparticle friction angle defined in CohFrictMat anf the internal/emergent/macroscopic friction angle of the DEM material/material.

The calibration consists in adjusting the interparticle friction angle (the one in CohFrictMat) that produces an internal/emergent/macroscopic friction angle of 40°.

Sorry if all this was clear but, again, I want to make sure to understand the source of your issue.

Finally, you may want to have a look at this article: https://hal.univ-grenoble-alpes.fr/hal-01951842/document

Luc

Can you help with this problem?

Provide an answer of your own, or ask Amedeo Galletti for more information if necessary.

To post a message you must log in.