← Back to team overview

yade-users team mailing list archive

Re: [Question #678364]: Normal stiffness for sphere-polyhedra interaction

 

Question #678364 on Yade changed:
https://answers.launchpad.net/yade/+question/678364

Description changed to:
Hello,

I have a question on how the normal stiffness between polyhedra and
sphere is calculated. I'm using Yade 2018.02b.

I checked the code files:

In Ig2_Sphere_Polyhedra_ScGeom,
geom->radius1 = geom->radius2 = radius;

In Ip2_FrictMat_PolyhedraMat_FrictPhys,
void Ip2_FrictMat_PolyhedraMat_FrictPhys::go(const shared_ptr<Material>& pp1, const shared_ptr<Material>& pp2, const shared_ptr<Interaction>& interaction){
	const shared_ptr<FrictMat>& mat1 = YADE_PTR_CAST<FrictMat>(pp1);
	const shared_ptr<PolyhedraMat>& mat2 = YADE_PTR_CAST<PolyhedraMat>(pp2);
	Ip2_FrictMat_FrictMat_FrictPhys().go(mat1,mat2,interaction);
}

And in Ip2_FrictMat_FrictMat_FrictPhys(),
	Real Ea 	= mat1->young;
	Real Eb 	= mat2->young;
	Real Va 	= mat1->poisson;
	Real Vb 	= mat2->poisson;

	//half the harmonic average of the two stiffnesses, when (2*Ri*Ei) is the stiffness of a contact point on sphere "i"
	Real Kn = 2*Ea*Ra*Eb*Rb/(Ea*Ra+Eb*Rb);
	//same for shear stiffness
	Real Ks = 2*Ea*Ra*Va*Eb*Rb*Vb/(Ea*Ra*Va+Eb*Rb*Vb);

I ran an example script, the sphere has radius 0.5, young 1e9, and polyhedra has young 1e10.
The kn between them in the simulation is 5e8, which is not the result by  Kn = 2*Ea*Ra*Eb*Rb/(Ea*Ra+Eb*Rb)

Wondering if there is other equation to calculate the kn for the sphere-
polyhedra interaction? Thank you!

Best regards,
Cong

The following is the example script.

from yade import polyhedra_utils
import random

polyMat  = PolyhedraMat(young=1e10,poisson=.05)
frictMat = FrictMat(young=1e9,poisson=.2)

O.materials.append((polyMat,frictMat))

poly = polyhedra_utils.polyhedra(polyMat,(1,2,3)); poly.wire=False
sph1 = sphere((2,2,2),.5,material=frictMat)
sph2 = sphere((-2,0,0),.5,material=frictMat)
sph3 = sphere((0,-2,0),.5,material=frictMat)
O.bodies.append((
	poly,
	sph1,
	sph2,
	sph3
))
	
O.engines=[
   ForceResetter(),
   InsertionSortCollider([Bo1_Polyhedra_Aabb(),Bo1_Sphere_Aabb()]),
   InteractionLoop(
      [Ig2_Sphere_Polyhedra_ScGeom()], 
      [Ip2_FrictMat_PolyhedraMat_FrictPhys()],
      [Law2_ScGeom_FrictPhys_CundallStrack()],
   ),
   NewtonIntegrator(),
]

O.dt = 1e-7
O.step()

poly.state.blockedDOFs = 'xyzXYZ'
for s in (sph1,sph2,sph3):
	r=random.random
	s.state.vel = -10*(s.state.pos+Vector3(r(),r(),r()))

from yade import qt
qt.View()

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