yade-users team mailing list archive

[Raphaël Maurin <question429604@xxxxxxxxxxxxxxxxxxxxx>]

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

Raphaël Maurin gave more information on the question:
Hi,

Even though I am not quite sure to understand everything, I found  the
origin of the problem by comparing the Law2_ScGeom_ViscElPhys_Basic to
the Law2_ScGeom_FrictPhys_CundallStrack. In the following of the
message, I focus only on the elastic part of
Law2_ScGeom_ViscElPhys_Basic, considering a case with zero viscous
damping.

When it is not periodic, there are not much difference in the evaluation of the normal and shear elastic contact force between the two contact laws, even if sign conventions are different. There is only a difference in the application of the contact force to the interacting elements. 
In the Law2_ScGeom_ViscElPhys_Basic, force and torques are evaluated and directly added to the particles following:
force = phys.normalForce + shearForce + shearForceVisc;
torque1 = -c1x.cross(force);
torque2 = c2x.cross(force);

with c1x = (geom.contactPoint - de1.pos) defined between the position of the particle (de1.pos) and the contact point.
and then

addForce (id1,-force,scene);
addForce (id2, force,scene);
addTorque(id1, torque1,scene);
addTorque(id2, torque2,scene);


In the  Law2_ScGeom_FrictPhys_CundallStrack it is done as: 
	if (!scene->isPeriodic && !sphericalBodies) {
		State* de1 = Body::byId(id1,scene)->state.get();
		State* de2 = Body::byId(id2,scene)->state.get();
		applyForceAtContactPoint(-phys->normalForce-shearForce, geom->contactPoint, id1, de1->se3.position, id2, de2->se3.position);}
	else {//we need to use correct branches in the periodic case, the following apply for spheres only
		Vector3r force = -phys->normalForce-shearForce;
		scene->forces.addForce(id1,force);
		scene->forces.addForce(id2,-force);
		scene->forces.addTorque(id1,(geom->radius1-0.5*geom->penetrationDepth)* geom->normal.cross(force));
		scene->forces.addTorque(id2,(geom->radius2-0.5*geom->penetrationDepth)* geom->normal.cross(force));
	}

Different things regarding that:
- I didn't succeed to find the location of the function applyForceAtContactPoint, but I believe that it is doing exactly the same as in the viscous law where the force and torque are applied taking into account the branch vectors c1x and c2x defined with respect to the contact point. 
- We can see that when the scene is periodic (or when we consider spheres), the force is however not applied with respect to the contact point in the Law2_ScGeom_FrictPhys_CundallStrack, but with respect to the radius and the penetration depth. I think these two methods are completely equivalent for two spheres of the same size. 

If I am right with these two points, there are no differences between
the two contact laws regarding the elastic contribution in non-periodic
frameworks and when considering monodisperse particles.

However, in the periodic case, while there is no change in the
evaluation of the branch vector in Law2_ScGeom_ViscElPhys_Basic, it is
different in the Law2_ScGeom_FrictPhys_CundallStrack where it is assumed
that the second part of the if statement apply for all interactions.
While this is a priori not appropriate for all interactions (such as an
interaction between a sphere and a wall), this formulation seems to
avoid some troubles. Indeed, with the Law2_ScGeom_ViscElPhys_Basic not
modified formulation, the partices tend to spin a lot leading to
unphysical behavior.

Could someone explain me why this trick is used for periodic BC and how
it solves this spinning problem ?

Replacing the branch vector formulation based on the contact point by
the one based on the penetration depth such as done in
Law2_ScGeom_FrictPhys_CundallStrack, solves the whole problem, even when
re-activating the viscous damping.

Last remark: The way the contact force is evaluated in
Law2_ScGeom_ViscElPhys_Basic does not take into account the implemented
ratcheting effect correction, and I believe it would be better to base
it on the same formulation as the CundallStrack contact law to account
for that. I can do it, but in order to do that, I need to know how to
access the function getIncidentVel of ScGeom inside the contact law. Can
someone tell me how to do that ?

Thanks !
Cheers, 

Raphael

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