yade-users team mailing list archive

[chiara modenese <c.modenese@xxxxxxxxx>]

I am testing my script with the class you committed but there is something
wrong. In the contact law I used to "play" with the damping in the shear
direction I coded something as in the attached file. Please look at that (I
changed the terminology but it is clear what it does) so you can check if
you are doing well.
I do not commit my law simply because it is quite a mess and I was playing
with it to clarify this issue around the damping. Please update your law
following the procedure as I suggested.

Thanks, Chia

///////////////////////////////////////// SHEAR FORCE (as in Cundall)
///////////////////////////////////////// ROTATION OF SHEAR FORCE
	Vector3r axis;
        Real angle;
	Vector3r& shearElastic_cum = phys->shearElastic_cum; // alias for phys->shearElastic_cum
	axis = phys->prevNormal.Cross(scg->normal); // axis of rotation between previous and current normal
	shearElastic_cum -= shearElastic_cum.Cross(axis); // minus shearForce that is perpendicular (outward) to contact (not tangent)
	angle = dt*0.5*scg->normal.Dot(de1->angVel+de2->angVel);
	axis = angle*scg->normal;
	shearElastic_cum -= shearElastic_cum.Cross(axis); // minus shearForce that is perpendicular outward to contact (not tangent)
///////////////////////////////////////// END OF ROTATION - START SHEAR FORCE
	Vector3r x = scg->contactPoint;
	Vector3r c1x = (x - de1->se3.position); // distance of contact point from centre1
	Vector3r c2x = (x - de2->se3.position); // distance of contact point from centre2
	// The delta relative velocity as it performed is an INCREMENT. In fact we get it considering the relative motion "at contact"
	Vector3r deltaRelVelocity = (de1->vel + de1->angVel.Cross(c1x)) - (de2->vel + de2->angVel.Cross(c2x));
	Real deltaNormVelocity = scg->normal.Dot(deltaRelVelocity);
	Vector3r deltaShearVelocity = deltaRelVelocity - deltaNormVelocity*scg->normal; // shear component of the relative velocity
	Vector3r deltaShearDisplacement = deltaShearVelocity*dt; // in Dem3DofGeom displacementT() is cumulative
	Vector3r deltaShearElastic = phys->ks*deltaShearDisplacement; // incremental elastic shear force
	shearElastic_cum = shearElastic_cum + deltaShearElastic; // total elastic shear force
	Vector3r shearViscous = cs*deltaShearVelocity; // shear force due to viscous damping (same sign convention as for shear force since we oppose the motion)

	// Account for global damping in the normal direction (update normal force)
	phys->normalForce += (deltaNormVelocity*cn)*scg->normal; // here the sing is plus (but we are oppising the motion, in fact we copmute the reltive velocity as (v1-v2); if it was (v2-v1) we should have put the minus sign

 ///////////////////////////////////////// MOHR-COULOMB LAW
	Real maxFs = phys->normalForce.SquaredLength()*std::pow(phys->tangensOfFrictionAngle,2);
	if (shearElastic_cum.SquaredLength() > maxFs)
	{
	Real ratio = Mathr::Sqrt(maxFs)/shearElastic_cum.Length(); 
	shearElastic_cum *= ratio; 
	phys->shearForce = shearElastic_cum;
	}
	else
	{phys->shearForce = shearElastic_cum + shearViscous;} // add viscous force in the shear direction if we passed the failure criterion
	phys->shearElastic_cum = shearElastic_cum; // store value

///////////////////////////////////////// APPLY FORCES
	// Here it is necessary to follow the sign convention
	Vector3r f = phys->normalForce + phys->shearForce;
	ncb->forces.addForce(id1,-f);
	ncb->forces.addForce(id2,f); // force that body_1 applies to body_2
	ncb->forces.addTorque(id1,-c1x.Cross(f)); // about axis perpendicular to normal and force
	ncb->forces.addTorque(id2,(c2x).Cross(f));

	phys->prevNormal = scg->normal; // update normal direction
}