# yade-users team mailing list archive

## [Question #699779]: Using Hertz - Mindlin Contact Law

 Thread Previous • Date Previous • Date Next • Thread Next

```New question #699779 on Yade:

Hi everybody
I want to carry out triaxial simulation on sand with considering a non-linear contact law (Hertz-Mindlin). I mention the code  in below which is working very well with Cohesion moment law but in terms of Hertz- MIndlin there are some problems in the case of creating a pack  such as an unacceptable unbalanced (almost 0.67) force for desire void ratio which is around 0.83.
Do you know what should I do to well-create a pack with this non-linear contact law ?
Also, I need to get to know the Hertz-Mindlin equations that are using on YADE, do you know the references of this contact law?

# encoding: utf-8

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 *

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
rate=5

psdSizes,psdCumm=[0.00012,0.00014,0.00018,0.00021,],[0.0030,0.039,0.21,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=10,mass=True)

## create material #0, which will be used as default
O.materials.append(FrictMat(young=5e9,poisson=.3,frictionAngle=0,density=0,label='frictionless'))

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

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

newton=NewtonIntegrator(damping=0)

O.engines=[
ForceResetter(),
InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
InteractionLoop(
[Ig2_Sphere_Sphere_ScGeom6D(),Ig2_Box_Sphere_ScGeom()],
[Ip2_FrictMat_FrictMat_MindlinPhys()],
),
GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
triax,
newton
]

#s11=[];e22=[]
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), 'Size=', (triax.width))
if unb<0.01 and abs(triax.goal1-triax.meanStress)/abs(triax.goal1)<0.01 and (triax.porosity)/(1-triax.porosity)<0.8399999999 :
break
print ("Particle.Distribution.is.stable")
########################################################
##Turn into triaxil test
########################################################

finalFricDegree = 18.0
triax.internalCompaction=False
while 1:
O.run(1000,True)
unb=unbalancedForce()
if unb<0.01 and abs(triax.goal1-triax.meanStress)/abs(triax.goal1)<0.01:
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]

triax.depth0=triax.depth
triax.height0=triax.height
triax.width0=triax.width
triax.goal1=triax.goal3=confiningS
triax.goal2=-rate
triax.wall_bottom_activated=0
#O.engines[2].lawDispatcher.functors[1].always_use_moment_law=True
#cohesiveLaw.always_use_moment_law=False

Exaxil=[5.98e-04,6.11e-04,0.00125,0.0013,0.00249,0.0031,0.00491,0.00729,0.010,0.017,0.025]
Exdev=[8.69,16.30,51.09,81.52,92.40,110.88,142.40,170.67,196.77,218.54,229.44]

ee=[];e22=[];s22=[]
file=open('Results'+'.txt',"w")

while 1:
O.run(1000,True)
ee.append((triax.porosity)/(1-triax.porosity))
s22.append(((-triax.stress(triax.wall_top_id)[1])-(-triax.stress(triax.wall_right_id)[0]-triax.stress(triax.wall_front_id)[2])/2.0)*1e-3)
e22.append((-triax.strain[1]))
print('ax=',(-triax.strain[1]), 'deviatoric=',(((-triax.stress(triax.wall_top_id)[1])-(-triax.stress(triax.wall_right_id)[0]-triax.stress(triax.wall_front_id)[2])/2.0)*1e-3), 'ee=', (triax.porosity)/(1-triax.porosity), 's11=',((-triax.stress(triax.wall_top_id)[1])*1e-3), 's33=',(((-triax.stress(triax.wall_right_id)[0])-(triax.stress(triax.wall_front_id)[2]))/2.0)*1e-3)
if abs(-triax.strain[1])>0.1 :
break
print(" Triaxil test is finished ")

file.write(str(s22)+" "+str(e22)+"\n")
file.write(str(s22)+" "+str(ee)+"\n")
file.close()

plt.figure()
plt.xlabel('Axial strain')
plt.ylabel('deviatoric')
plt.grid(True)
plt.plot(Exaxil,Exdev,linestyle='-',marker="*", markersize=8,color='g',linewidth=2)
plt.plot(e22,s22,linestyle='-',marker="^", markersize=8,color='r',linewidth=2)
plt.legend(('Exprimental data','simulation data'),loc='upper right', shadow=True)

plt.show()

--

```

 Thread Previous • Date Previous • Date Next • Thread Next