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

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

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

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

## Question information

Language:
English Edit question
Status:
Solved
For:
Assignee:
No assignee Edit question
Solved by:
Jan Stránský
Solved:
Last query:
 Revision history for this message Jan Stránský (honzik) said on 2022-04-07: #1

Hello,

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
# same for id2
###

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:
###

Cheers
Jan

 Revision history for this message Lin 1999 (sarcasticli) said on 2022-04-07: #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 on 2022-04-07: #3

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