← Back to team overview

yade-users team mailing list archive

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

 

New question #699779 on Yade:
https://answers.launchpad.net/yade/+question/699779

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 *

utils.readParamsFromTable(seed=1,num_spheres=10000,compFricDegree = 25.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
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=radians(compFricDegree),density=2600,label='spheres'))
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)


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)

O.engines=[
	ForceResetter(),
	InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
	InteractionLoop(
		[Ig2_Sphere_Sphere_ScGeom6D(),Ig2_Box_Sphere_ScGeom()],
		[Ip2_FrictMat_FrictMat_MindlinPhys()],
		[Law2_ScGeom_MindlinPhys_Mindlin(includeAdhesion=False,includeMoment=False)]
	),
	GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
	VTKRecorder(iterPeriod=1000,fileName="./export",recorders=['spheres','intr','color'],dead=1,label='VTK'),
	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 
setContactFriction(radians(finalFricDegree))
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.stressMask=5
triax.goal1=triax.goal3=confiningS
triax.goal2=-rate
triax.wall_bottom_activated=0
#O.engines[2].lawDispatcher.functors[1].always_use_moment_law=True
VTK.dead=1
#cohesiveLaw.always_use_moment_law=False
#recorder.dead=0


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()



-- 
You received this question notification because your team yade-users is
an answer contact for Yade.