Yade

Hertz Mindlin Contact Law including Rolling

Asked by Hossein

Hi every body

I need to rolling resistance in order to increase the stiffness of the sample during a triaxial test. I have already probe with below parameters such as young modules and Poisson ratio and the result is the same with a simulation including rolling stiffness (krot) and coefficent of bending moment (eta) was performed. So my question is that :

1) How to should I assign rolling and twisting in a simulation with considering Hertz Mindlin contact law ?

Also

2) I want to active this contact law after the sample achieve desire void ratio, since the computation take a long time to achieve desire void ratio, So How can I do that?

Here in below I mention my code where the contact law must be activated:

O.materials.append(FrictMat(young=9e10,poisson=.25,frictionAngle=radians(compFricDegree),density=2600,label='spheres'))
O.materials.append(FrictMat(young=9e10,poisson=.25,frictionAngle=0,density=0,label='frictionless'))

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

O.bodies.append([utils.sphere(center,rad,material='spheres') for center,rad in sp])

triax=TriaxialStressController(
 internalCompaction=False,
 goal1=confiningS,
 goal2=confiningS,
 goal3=confiningS,
 max_vel=10,
 label="triax"
)

newton=NewtonIntegrator(damping=0.7)

O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom6D(),Ig2_Box_Sphere_ScGeom()],
  [Ip2_FrictMat_FrictMat_MindlinPhys(krot=.4,eta=.2)],
  [Law2_ScGeom_MindlinPhys_Mindlin(includeAdhesion=False,includeMoment=True)]
 ),
 GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
 VTKRecorder(iterPeriod=1000,fileName="./export",recorders=['spheres','intr','color'],dead=1,label='VTK'),
 triax,
 newton
]

Question information

Language:
English Edit question
Status:
Solved
For:
Yade Edit question
Assignee:
No assignee Edit question
Solved by:
Karol Brzezinski
Solved:
Last query:
Last reply:
Revision history for this message
Best Karol Brzezinski (kbrzezinski) said : 		#1

Hi,

MWE (Minimum Working Example)[1] would be helpful with emphasis on "working". When running your code, the first problem I encounter is an undefined variable (compFricDegree).
I cannot help with twisting, but rolling seems to be properly included in your code.

If you want to start without this, and activate it later maybe you should add a label to your engine, and change its properties by referring to this label. Not tested, but this workflow would look like this:

#### in engines set a label to Law2
Law2_ScGeom_MindlinPhys_Mindlin(includeMoment=False, label = 'Mindlin_Law')

(... run the simulation until some condition is reached ...)
#change the engine
Mindlin_Law.includeMoment=True
O.interactions.clear()# I would reset all the interactions just in case

(... run the rest of the simulation ...)

Cheers,
Karol

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

Revision history for this message
Hossein (hossein75) said (last edit ): 		#2

 Thanks Karol
your comment was really practical

I have another problem that when I wanna use Twophase flow engine with adhesion (gamma parameters which I must use) this engine work not very well. So, I have to have two different engine, one of them (Cundall Contact Law) must be used by the Twophase and another one (Hertz Mindlin )one must be used terms of loading. The way that came on my mind is that I should turn down an engine while another one must be replaced. I mention my code in below and also mark a line that I think the new engine
(Hertz Mindlin Include Adhesion ) must be add.

Actually, I do need to use Hertz Mindlin (Include Adhesion with defining surface energy (gamma)) in Stage Two in below cod.

###################################
from mpl_toolkits.axes_grid1 import host_subplot
import matplotlib; matplotlib.rc('axes',grid=True)
import matplotlib.pyplot as plt
from yade import pack
import pylab
from numpy import *
import numpy as np
from math import *

utils.readParamsFromTable(seed=1,num_spheres=10000,compFricDegree =8.0)
from yade.params import table

seed=table.seed
num_spheres=table.num_spheres# number of spheres
compFricDegree = table.compFricDegree # initial contact friction during the confining phase (will be decreased during the REFD compaction process)
confiningS=-1e5

young1=10e7
psdSizes,psdCumm=[0.00509,0.00704,0.0101,0.0127,0.0166,0.0222,0.0338],[.322,.376,.452,.533,.632,.726,1]
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.0001,0.0001,0.0001)
sp.makeCloud(minCorner=mn,maxCorner=mx,psdSizes=psdSizesMeter,psdCumm=psdCumm,distributeMass=True,num=num_spheres,seed=seed)
sp.psd(bins=4000,mass=True)
key='_Youngmadules_'+str(young1)

O.materials.append(FrictMat(young=70e7,poisson=.2,frictionAngle=radians(compFricDegree),density=2600,label='spheres'))
O.materials.append(FrictMat(young=70e7,poisson=.2,frictionAngle=0,density=0,label='frictionless'))

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

O.bodies.append([utils.sphere(center,rad,material='spheres') for center,rad in sp])

triax=TriaxialStressController(
 internalCompaction=False,
 goal1=confiningS,
 goal2=confiningS,
 goal3=confiningS,
 max_vel=10,
 label="triax"
)
newton=NewtonIntegrator(damping=0.4)

O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom6D(),Ig2_Box_Sphere_ScGeom()],
  [Ip2_FrictMat_FrictMat_FrictPhys(),Ip2_FrictMat_FrictMat_MindlinPhys(betan=0.04,betas=0.04,gamma=7)],
  [Law2_ScGeom_FrictPhys_CundallStrack(label = 'Cundall_Law'),Law2_ScGeom_MindlinPhys_Mindlin(includeAdhesion=True,includeMoment=False,label
  = 'Mindlin_Law'),

  ]
 ),
 GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
 VTKRecorder(iterPeriod=1000,fileName="./export",recorders=['spheres','intr','color'],dead=1,label='VTK'),
 triax,
 TwoPhaseFlowEngine(label="unsat"),
 newton
]

always_use_Cundall_Law=True
always_use_Mindlin_Law=False
while 1:
  O.run(1000,True)
  unb=unbalancedForce()
  print ('unb=', unb,'meanstress:',abs(triax.goal1-triax.meanStress)/abs(triax.goal1),'e=', triax.porosity/(1-triax.porosity), )
  if unb<0.5 and abs(triax.goal1-triax.meanStress)/abs(triax.goal1)<0.5 and (triax.porosity)/(1-triax.porosity)<0.979999999999999999 :
    break
print ("Particle.Distribution.is.stable")
##################################

##################################
finalFricDegree = 9.0
triax.internalCompaction=False
setContactFriction(radians(finalFricDegree))
while 1:
  O.run(1000,True)
  unb=unbalancedForce()
  if unb<0.5 and abs(triax.goal1-triax.meanStress)/abs(triax.goal1)<0.5:
    break

ei0=-triax.strain[0];ei1=-triax.strain[1];ei2=-triax.strain[2]
si0=-triax.stress(0)[0];si1=-triax.stress(2)[1];si2=-triax.stress(4)[2]

from yade import plot
O.engines=O.engines+[PyRunner(iterPeriod=20,command='history()',label='recorder')]
## a function saving variables
def history():
        plot.addData(e11=-triax.strain[0]-ei0, e22=-triax.strain[1]-ei1, e33=-triax.strain[2]-ei2,
  s11=-triax.stress(0)[0]-si0,
  s22=-triax.stress(2)[1]-si1,
  s33=-triax.stress(4)[2]-si2,
         ee=(((triax.porosity)/(1-triax.porosity))*1e-2),
         syy=((-triax.stress(3)[1])*1e-3),

         i=O.iter
 )
 #plot.addData( i=O.iter,ee=(triax.porosity)/(1-triax.porosity),s22=-triax.stress(3)[1]*1e-3)
 #plot.addData(e22=-triax.strain[1],t=O.time,s22=-triax.stress(2)[1],p=flow.MeasurePorePressure((0.5,0.5,0.5)))

##make nice animations:
#O.engines=O.engines+[PyRunner(iterPeriod=200,command='flow.saveVtk()')]

from yade import plot
plot.plots={'syy':('ee')}
plot.plot()
###################################
#Stage one ( Two phase flow )
###################################
print("Get saturation")

triax.stressMask=7
triax.goal1=triax.goal3=0 #confiningS
goalTop=triax.stress(3)[1]
triax.goal2=goalTop
triax.wall_bottom_activated=0

void=[];sy=[];Suction=[];Void2=[];valumstrain=[]

meanDiameter=(O.bodies[-1].shape.radius + O.bodies[6].shape.radius) / 2.

unsat.bndCondIsPressure=[0,0,1,1,0,0]
unsat.bndCondValue=[0,0,-1e8,0,0,0]
unsat.isPhaseTrapped=True
unsat.initialization()
unsat.surfaceTension = 0.0728

#file=open('Results'+key+'.txt',"w")
for pg in arange(1.0e-5,2,0.7):
  unsat.bndCondValue=[0,0,(-1.0)*pg*unsat.surfaceTension/meanDiameter,0,0,0]
  unsat.invasion()
  unsat.computeCapillaryForce()
  for b in O.bodies:
    O.forces.setPermF(b.id, unsat.fluidForce(b.id))
  while 1:
    O.run(1000,True)
    unb=unbalancedForce()
    if unb<0.5:
      break
  print ('PC:',((-unsat.bndCondValue[2])*1e-3), 'Sw:',(unsat.getSaturation(False)),'e:',(triax.porosity)/(1-triax.porosity), 'e33=', (triax.strain[2]), 'e22=', (-triax.strain[1]), 's22=', (-triax.stress(3)[1]*1e-3))

print ('----------LOADING--------------')
###############################################################
###Stage Two ( Loading )
###############################################################
O.engines[2].lawDispatcher.functors[1].always_use_Cundall_Law=False
O.engines[2].lawDispatcher.functors[1].always_use_Mindlin_Law=True

betan=0.04
betas=0.04
gamma=0.7
Mindlin_Law.includeMoment=True
Mindlin_Law.includeAdhesion=True

unsat.bndCondIsPressure=[0,0,0,1,0,0]
triax.stressMask=2
triax.goal1=triax.goal3=0
goalTop=triax.stress(3)[1]
triax.goal2=goalTop
triax.internalCompaction=False
triax.wall_bottom_activated=False
loadingMatrix=[-5e3,-2e5,-8.3e5,-4e5]

for i in arange(0,len(loadingMatrix),1):
  triax.goal2=loadingMatrix[i]
  O.run(1000,True)
  void.append((triax.porosity)/(1-triax.porosity))
  sy.append((-triax.stress(3)[1])*1e-3)
  valumstrain.append((-triax.strain[1]))
  print ('e=',((triax.porosity)/(1-triax.porosity)), 's22=',(-triax.stress(3)[1]*1e-3), 'Load:',(loadingMatrix[i]), 'Sw:',(unsat.getSaturation(False)), 'e22=',(-triax.strain[1]), 'PC=',((-unsat.bndCondValue[2])*1e-3))

plt.figure()

plt.subplot(221)
plt.grid(True)
plt.plot(sy,void,linestyle='-',marker="^", markersize=8,color='r',linewidth=2)
plt.legend(('Exprimental data','simulation data'),loc='upper right', shadow=True)

plt.show()

Revision history for this message
Hossein (hossein75) said : 		#3

Thanks Karol Brzezinski, that solved my question.