yade-users team mailing list archive

[Lin 1999 <question701222@xxxxxxxxxxxxxxxxxxxxx>]

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

Hi! 
I have a quetion with getting the contact radius of the contact particles.
A Hertz-Minglin contact model is used to describe sphere-sphere interactions when I prepared the sample using "TriaxialStressController()".
I want to calculate the radius of the overlaped area between two adjacent particles. But I do not know how to do. 
Are there some methods to achieve it? Thanks!

Here is my code to prepare the sample.
##############
from __future__ import print_function
from yade import pack,plot,os,timing
import matplotlib;matplotlib.rc('axes',grid=True)
import pylab

#Material constants
density=2650
poisson=0.31
young=29e9
damp=0.25
number=5000
frictionangle=0

#Wall constants
walldensity=0
wallfrictionangle=0
wallpoisson=0.5
wallyoung=30000e9
mn=Vector3(0.,0.,0.)
mx=Vector3(0.01,0.01,0.01)

compress=-5000
rate=0.1

O.materials.append(FrictMat(young=young,poisson=poisson,frictionAngle=radians(frictionangle),density=density,label='spheres')) 
O.materials.append(FrictMat(young=wallyoung,poisson=wallpoisson,frictionAngle=radians(wallfrictionangle),density=walldensity,label='walls'))
wallIds=O.bodies.append(aabbWalls([mn,mx],thickness=0.00001,material='walls')) 

# PSD of the particles
psdSizes,psdCumm=[0.00025,0.0004,0.0006,0.00075],[0,0.15,0.85,1.0]
pylab.plot(psdSizes,psdCumm,label='precribed mass PSD')

#Add particles
sp = pack.SpherePack()
sp.makeCloud(mn,mx,-1,0.33,number,False,psdSizes=psdSizes,psdCumm=psdCumm,distributeMass=True)
pylab.plot(*sp.psd(bins=30,mass=True),label='PSD of (free) %d random spheres'%len(sp))
pylab.legend()
sp.toSimulation(material='spheres')

#Color
for b in O.bodies:
    if isinstance(b.shape,Sphere):
       r=b.shape.radius
       if r>0.00025 and r<0.000375:
          b.shape.color=(1.0,0/255.,51/255.)
       if r>0.0 and r<0.00025:
          b.shape.color=(1.0,1.0,0)

#Energy
O.trackEnergy=True

#Defining triaxil engines
triax=TriaxialStressController(
   wall_bottom_id = wallIds[2],
   wall_top_id = wallIds[3],
   wall_left_id = wallIds[0],
   wall_right_id = wallIds[1],
   wall_back_id = wallIds[4],
   wall_front_id = wallIds[5],
   thickness=0.00001,
   internalCompaction=False, 
   stressMask=7, 
   goal1=compress,
   goal2=compress,
   goal3=compress,
)

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_MindlinPhys()],
            [Law2_ScGeom_MindlinPhys_Mindlin()]
  ),
GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8,defaultDt=4*utils.PWaveTimeStep()),
   triax,
   newton,
   TriaxialStateRecorder(iterPeriod=5000,file='WallStress'),
   PyRunner(realPeriod=1,command='checkUnbalanced()',label='check'),
   PyRunner(command='addPlotData()',iterPeriod=5000,label='record')
]
#############

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