Calculate the radius of the overlaped area between two adjacent particles.

Asked by Lin 1999

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')
]
#############

Question information

Language:
English Edit question
Status:
Solved
For:
Yade Edit question
Assignee:
No assignee Edit question
Solved by:
Jan Stránský
Solved:
Last query:
Last reply:
Revision history for this message
Best Jan Stránský (honzik) said :
#1

Hello,

please try to make your code a MWE [1].
M = minimal
two (or just a few) spheres are sufficient to target your problem. Moreover, with fixed spheres, you may already know the expected result.

Knowing sphere centers and radii, it can be computed with simple math
###
for i in O.interactions:
    b1 = O.bodies[i.id1]
    center1 = b1.state.pos
    radius1 = b.shape.radius
    # same for id2
    radius = getRadiusFrom(center1,radius1,center2,radius2) # some simple math
###

Or you can try to use a builtin method. You are using MindlinPhys. Searching MindlinPhys for something about radius gives [2], is it what you are looking for?
###
for i in O.interactions:
    contactRadius = i.phys.radius
###

Cheers
Jan

[1] https://www.yade-dem.org/wiki/Howtoask
[2] https://yade-dem.org/doc/yade.wrapper.html#yade.wrapper.MindlinPhys.radius

Revision history for this message
Lin 1999 (sarcasticli) said :
#2

Sorry. I did not noticed the code is too complex. I will pay attention to the description next time. I ignored it can be solved with math. This solve my difficulty. Thanks!

Cheers
lin

Revision history for this message
Lin 1999 (sarcasticli) said :
#3

Thanks Jan Stránský, that solved my question.