ConstitutiveLawHowto

How to create new constitutive law

Every interaction has 2 components: InteractionGeometry and InteractionPhysics. InteractionGeometry gives you information about various spatial characteristics of that interaction (such as normal and shear strain). InteractionPhysics has all other parameters and internal variables.

Constitutive law declaration

Constitutive law is a special "functor" that receives InteractionGeometry and InteractionPhysics based on their type. Every constitutive law declares the types it acts upon in its declaration.

Let's take the class Law2_Dem3Dof_Elastic_Elastic as an example: <source lang="cpp"> class Law2_Dem3Dof_Elastic_Elastic: public ConstitutiveLaw{ public: virtual void go(shared_ptr<InteractionGeometry>& _geom, shared_ptr<InteractionPhysics>& _phys, Interaction* I, MetaBody* rootBody); FUNCTOR2D(Dem3DofGeom,ElasticContactInteraction); REGISTER_CLASS_AND_BASE(Law2_Dem3Dof_Elastic_Elastic,ConstitutiveLaw); REGISTER_ATTRIBUTES(ConstitutiveLaw,/*nothing here*/); }; REGISTER_SERIALIZABLE(Law2_Dem3Dof_Elastic_Elastic); </source> What you see:

1. Every constitutive law derives from ConstitutiveLaw class
2. the virtual void go(...) line is common to all constitutive laws
3. FUNCTOR2D(Dem3DofGeom,ElasticContactInteraction) declares types for which it will be dispatched. Therefore, in our case, the InteractionGeometry will be of type Dem3DofGeom and InteractionPhysics will be of type ElasticContactInteraction.
4. REGISTER_CLASS_AND_BASE, REGISTER_ATTRIBUTES and REGISTER_SERIALIZABLE are macros related to type information and serialization. They are all mandatory.
5. The name Law2_Dem3Dof_Elastic_Elastic should means "functor dispatching based on 2 arguments", which are of type "Dem3Dof" (geometry) and Elastic(=ElasticContactLaw) for physics, with the epitet of "Elastic".

Constitutive law implementation

We are going to describe very simple constitutive law, which is linear elastic in normal sense and has Mohr-Coulomb friction in the shear sense (linear plasticity function).

Let's have a look at the source for Law2_Dem3Dof_Elastic_Elastic: <source lang="cpp"> void Law2_Dem3Dof_Elastic_Elastic::go(shared_ptr<InteractionGeometry>& ig, shared_ptr<InteractionPhysics>& ip, Interaction* contact, MetaBody* rootBody){ Dem3DofGeom* geom=static_cast<Dem3DofGeom*>(ig.get()); ElasticContactInteraction* phys=static_cast<ElasticContactInteraction*>(ip.get()); Real displN=geom->displacementN(); if(displN>0){rootBody->interactions->requestErase(contact->getId1(),contact->getId2()); return; } phys->normalForce=phys->kn*displN*geom->normal; Real maxFsSq=phys->normalForce.SquaredLength()*pow(phys->tangensOfFrictionAngle,2); Vector3r trialFs=phys->ks*geom->displacementT(); if(trialFs.SquaredLength()>maxFsSq){ geom->slipToDisplacementTMax(sqrt(maxFsSq)); trialFs*=maxFsSq/(trialFs.SquaredLength());} applyForceAtContactPoint(phys->normalForce+trialFs,geom->contactPoint,contact->getId1(),geom->se31.position,contact->getId2(),geom->se32.position,rootBody); } </source>

What happens here? Let's comment the code by pieces.

• At the beginning, we cast received arguments to the types declared by FUNCTOR2D. We are sure that they are of that type (taken care of by ConstitutiveLawDispatcher, resp. InteractionDispatchers), there we can use static_cast. It is not necessary to case to shared_ptr, as we will not keep reference: .get() gets raw pointer from the shared pointer, which we cast afterwards.

<source lang="cpp"> Dem3DofGeom* geom=static_cast<Dem3DofGeom*>(ig.get()); ElasticContactInteraction* phys=static_cast<ElasticContactInteraction*>(ip.get()); </source> (if you are curious: the virtual method go(...) has its arguments declared in the top class (ConstitutiveLaw), therefore cannot receive arguments directly cast to their types)

• In the second part, we compute normal force.
  • First, we ask for the normal displacement (relative to the equilibrium − initial − position); as this involves some (small) computation and we may need the number later, store it in local var displN.

<source lang="cpp"> Real displN=geom->displacementN(); </source>

• • If normal displacement is positive, the two bodies don't touch anymore. As this particular constitutive law is not cohesive, we may ask to erase this contact. As there will be no force resulting, we mya return from the call immediately.

<source lang="cpp"> if(displN>0){rootBody->interactions->requestErase(contact->getId1(),contact->getId2()); return; } </source>

• • Set NormalInteraction::normalForce to the desired value; phys->kn × displN (rigidity × displacement) is the force magnitude, multiplied by geom->normal (unit vector from body A (I->getId1()) to body B (I->getId2()), we get it as 3d vector. The orientation of the force is as it will apply on A, but we will reverse it on B later.

<source lang="cpp"> phys->normalForce=phys->kn*displN*geom->normal; </source>

• Now, regarding the shear force.
  • We get maxFsSq (maximum shear force Fs for our, already computed, normal force; actually radius of plastic surface for this Fn; squared), and "trial force" (like trial stress, but already multiplied by length and area). geom->displacementT() computes the vector from spatial positions of spheres; the Vector3r it returns is in the contact plane, therefore perpendicular to geom->normal.

<source lang="cpp"> Real maxFsSq=phys->normalForce.SquaredLength()*pow(phys->tangensOfFrictionAngle,2); Vector3r trialFs=phys->ks*geom->displacementT(); </source> (It is better here to get square, since square root is quite expensive operation and we need the value only for comparison now.)

• • If the "trial force" is beyond the plasticity surface, we perform plastic slip. This is done by (i) doing the slip itself by asking the Dem3DofGeom class, providing magnitude of the shear force requested, square-rooted now; it will not change the direction. (ii) We make the trialFs smaller by the factor sqrt(maxFsSq/trialFs.SquaredLength()) so that it returns to the plasticity surface. From this point on, trialFs is no longer "trial", but the real force.

<source lang="cpp"> if(trialFs.SquaredLength()>maxFsSq){ geom->slipToDisplacementTMax(sqrt(maxFsSq)); trialFs*=sqrt(maxFsSq/(trialFs.SquaredLength()));} </source>

• Forces are computed, now they have to be applied to the bodies. For this, there is utility function ConstitutiveLaw::applyForceAtContactPoint, which needs to receive lots of data. The first argument (normalForce+trialFs) is applied on body A, and applied in reverse to body B. Torques resulting from the fact that the force is located at the contact point (and not at body centroids) are computed automatically as well.

<source lang="cpp"> applyForceAtContactPoint(phys->normalForce+trialFs,geom->contactPoint,contact->getId1(),geom->se31.position,contact->getId2(),geom->se32.position,rootBody); </source>

Going further

Custom InteractionPhysics class

For most your application, you will need specific InteractionPhysics to hold your interaction parameters and contact internal variables. To have that, you need (again) a functor that will derive the interaction parameters from Body::physicalParameters. For an example, see the CpmPhys class. The functor that creates CpmPhys is called Ip2_CpmMat_CpmMat_CpmPhys (2d functor creating InteractionPhysics (=Ip2) taking CpmMat and CpmMat of 2 bodies, returning type CpmPhys).

In this case, I needed even special PhysicalParameters class (CpmMat, standing for "material") to hold a few parameters specific to each body.

There is a high-level overview of the CPM model in its header file, which can get you some idea about the more complicated cases.

InteractionGeometry

The Dem3DofGeom class is abstract class for computing geometry (displacements) on 2 bodies. It has the same interface for both Sphere-Sphere and Sphere-Facet contacts. the 3 degrees of freedom in the name mean normal and shear displacements.

If you need torsion and bending in your code, there is, as of now unimplemented, Dem6DofGeom class that sould additionaly feature bending and torsion as well. Ask on the mailing list to have it implemented, it would be quite easy to do.

Personal note: do not use SpheresContactGeometry, it lacks any intelligence on its own and makes you put lots of code into the constitutive law. Avoid that. Keep your classes simple and readable.