Facing problem in stabilizing the DEM simulation sample when using Hertz-mindlin contact model

Asked by Sarath Chandra Reddy on 2021-05-30

Dear Friends,
I am facing problem in stabilizing the DEM model sample and reaching the isotropic compression stage when I am using Hertz-mindlin contact model. The following is my code. Please run the code once and give your suggestions on where I need to correct the code.

####################
from yade import pack

############################################
### DEFINING VARIABLES AND MATERIALS ###
############################################
# The following 5 lines will be used later for batch execution
nRead=readParamsFromTable(
 num_spheres=10000,# number of spheres
 compFricDegree = 5, # contact friction during the confining phase
 key='_triax_base_', # put you simulation's name here
 unknownOk=True
)
from yade.params import table

num_spheres=table.num_spheres# number of spheres
key=table.key
targetPorosity = 0.368 #the porosity we want for the packing
compFricDegree = table.compFricDegree # initial contact friction during the confining phase (will be decreased during the REFD compaction process)
finalFricDegree = 30 # contact friction during the deviatoric loading
rate=-0.02 # loading rate (strain rate)
damp=0.2 # damping coefficient
stabilityThreshold=0.01 # we test unbalancedForce against this value in different loops (see below)
young=5e7 # contact stiffness
mn,mx=Vector3(0,0,0),Vector3(.02,.02,.02) # corners of the initial packing

## create materials for spheres and plates
O.materials.append(FrictMat(young=young,poisson=0.2,frictionAngle=radians(compFricDegree),density=2600,label='spheres'))
O.materials.append(FrictMat(young=young,poisson=0.2,frictionAngle=0,density=3000,label='walls'))

## create walls around the packing
walls=aabbWalls([mn,mx],thickness=0,material='walls')
wallIds=O.bodies.append(walls)

## use a SpherePack object to generate a random loose particles packing
sp=pack.SpherePack()

clumps=False #turn this true for the same example with clumps
if clumps:
 ## approximate mean rad of the futur dense packing for latter use
 volume = (mx[0]-mn[0])*(mx[1]-mn[1])*(mx[2]-mn[2])
 mean_rad = pow(0.09*volume/num_spheres,0.3333)
 ## define a unique clump type (we could have many, see clumpCloud documentation)
 c1=pack.SpherePack([((-0.2*mean_rad,0,0),0.5*mean_rad),((0.2*mean_rad,0,0),0.5*mean_rad)])
 ## generate positions and input them in the simulation
 sp.makeClumpCloud(mn,mx,[c1],periodic=False)
 sp.toSimulation(material='spheres')
 O.bodies.updateClumpProperties()#get more accurate clump masses/volumes/inertia
else:
 sp.makeCloud(mn,mx,num=num_spheres,rMean=0.0005,periodic=False, seed=1)
 #sp.makeCloud(mn,mx,-1,0.3333,num_spheres,False, 0.95,seed=1) #"seed" make the "random" generation always the same
 O.bodies.append([sphere(center,rad,material='spheres') for center,rad in sp])
 print(len(sp),' particles generated.')
 #or alternatively (higher level function doing exactly the same):
 #sp.toSimulation(material='spheres')

############################
### DEFINING ENGINES ###
############################

triax=TriaxialStressController(
 ## TriaxialStressController will be used to control stress and strain. It controls particles size and plates positions.
 ## this control of boundary conditions was used for instance in http://dx.doi.org/10.1016/j.ijengsci.2008.07.002
 maxMultiplier=1.+9e4/young, # spheres growing factor (fast growth)
 finalMaxMultiplier=1.+2e3/young, # spheres growing factor (slow growth)
 thickness = 0,
 stressMask = 7,
 internalCompaction=False, # If true the confining pressure is generated by growing particles
)

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()],
  #[Law2_ScGeom_FrictPhys_CundallStrack()]
  [Ip2_FrictMat_FrictMat_MindlinPhys( betan =0.2 , betas = 0.2, krot = 0.2 , eta = 0.2 )],
  [Law2_ScGeom_MindlinPhys_Mindlin()]
 ),
 GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.25),
 triax,
 TriaxialStateRecorder(iterPeriod=100,file='WallStresses'+table.key),
 newton
]

#Display spheres with 2 colors for seeing rotations better
#Gl1_Sphere.stripes=0
#if nRead==0: yade.qt.Controller(), yade.qt.View()
triax.max_vel=1
#the value of (isotropic) confining stress defines the target stress to be applied in all three directions
triax.goal1=triax.goal2=triax.goal3=-100000

while 1:
  O.run(1000, True)
  ##the global unbalanced force on dynamic bodies, thus excluding boundaries, which are not at equilibrium
  unb=unbalancedForce()
  print ('unbalanced force:',unb,' mean stress: ',triax.meanStress)
  if unb<stabilityThreshold and abs(-100000-triax.meanStress)/100000<0.001:
    break
print ("Particle Distribution is Stable")
O.save('confinedState'+key+'.yade.gz')
print ("### Isotropic state saved ###")

###################################################
### REACHING A SPECIFIED POROSITY PRECISELY ###
###################################################

import sys #this is only for the flush() below
while triax.porosity>targetPorosity:
 ## we decrease friction value and apply it to all the bodies and contacts
 compFricDegree = 0.95*compFricDegree
 setContactFriction(radians(compFricDegree))
 print ("\r Friction: ",compFricDegree," porosity:",triax.porosity)
 sys.stdout.flush()
 ## while we run steps, triax will tend to grow particles as the packing
 ## keeps shrinking as a consequence of decreasing friction. Consequently
 ## porosity will decrease
 O.run(500,1)

#O.save('compactedState'+key+'.yade.gz')
#print "### Compacted state saved ###"

##############################
### DEVIATORIC LOADING ###
##############################

##We move to deviatoric loading, let us turn internal compaction off to keep particles sizes constant
triax.internalCompaction=False

## Change contact friction (remember that decreasing it would generate instantaneous instabilities)
setContactFriction(radians(finalFricDegree))

##set stress control on x and z, we will impose strain rate on y
triax.stressMask = 5
##now goal2 is the target strain rate
triax.goal2=rate
## we define the lateral stresses during the test, here the same 10kPa as for the initial confinement.
triax.goal1=-100000
triax.goal3=-100000

##we can change damping here. What is the effect in your opinion?
newton.damping=0.1

##Save temporary state in live memory. This state will be reloaded from the interface with the "reload" button.
O.saveTmp()

#####################################################
### Example of how to record and plot data ###
#####################################################

from yade import plot

### a function saving variables
def history():
 plot.addData(e11=-triax.strain[0], e22=-triax.strain[1], e33=-triax.strain[2],
   ev=-triax.strain[0]-triax.strain[1]-triax.strain[2],
   s11=-triax.stress(triax.wall_right_id)[0],
   s22=-triax.stress(triax.wall_top_id)[1],
   s33=-triax.stress(triax.wall_front_id)[2],
   i=O.iter)
 plot.saveDataTxt('calibration-hz.txt')

if 1:
  ## include a periodic engine calling that function in the simulation loop
  O.engines=O.engines[0:5]+[PyRunner(iterPeriod=20,command='history()',label='recorder')]+O.engines[5:7]
  ##O.engines.insert(4,PyRunner(iterPeriod=20,command='history()',label='recorder'))
else:
  O.engines[4]=PyRunner(iterPeriod=20,command='history()',label='recorder')

O.run(100000,True)

#plot.plots={'e22':('s11','s22','s33',None,'ev')}

## display on the screen (doesn't work on VMware image it seems)
#plot.plot()

Thanking you in advance,
Sarath

Question information

Language:
English Edit question
Status:
Open
For:
Yade Edit question
Assignee:
No assignee Edit question
Last query:
Last reply:
Revision history for this message
Launchpad Janitor (janitor) said :
#1

This question was expired because it remained in the 'Open' state without activity for the last 15 days.

Revision history for this message
Jan Stránský (honzik) said :
#2

Hello,

> I am facing problem

please be more specific [1]. Do you get errors? Does it take too much time? The target cannot be reached? ... ?

> when I am using Hertz-mindlin contact model

with other model(s) it is ok or not? (like why you mention explicitly "using Hertz-Mindlin model").

> Please run the code

please describe the problem such that running the code is not the primary source of information [1].

Cheers
Jan

[1] https://www.yade-dem.org/wiki/Howtoask

Revision history for this message
Sarath Chandra Reddy (sarath4321cityu) said :
#3

Thank you Jan for the reply.
>please be more specific [1]. Do you get errors? Does it take too much time? The target cannot be reached? ... ?

I am not getting any error but the unbalanced force and mean stress are not reaching equilibrium even after long time (say even after few days).

>when I am using Hertz-mindlin contact model - with other model(s) it is ok or not? (like why you mention explicitly "using Hertz-Mindlin model").

I think this especially when I am using Hertz-Mindlin model as I only changed small things in original code of Dr.Bruno (mainly about the contact model).

Please give your advice/suggestions on any measures to be taken to attain stability quickly.

best wishes
Sarath

Can you help with this problem?

Provide an answer of your own, or ask Sarath Chandra Reddy for more information if necessary.

To post a message you must log in.