yade-users team mailing list archive

[Bruno Chareyre <bruno.chareyre@xxxxxxxxxxx>]

> 2. current "trial" shear strain εT is computed*, with corresponding
> trial stress σT(εT); for admissible state f(σN,|σT(εT)|)<=0, there is no
> dissipation and the contact is done for this timestep. In case of
> non-admissible f(σN,|σT|)>0, shear strain is modified** to have new
> value εT2 so that f(σN,|σT(εT2)|=0; dissipated energy increment is |
> εT-εT2|*σT(εT2), which is added to cummulative dissipated energy. In the
> Cpm model, plastic flow rule is non-associated and σN doesn't have to be
> recomputed again, after changing σT.***
>   
> > Worth than that : unload the normal component with null shear velocity. 
> > It will meet the plasticity criterion at some point (assuming purely 
> > frictional law). Some elastic energy will be transfered to plasticity 
> > with total dUs=0. I see this special case only now... I'll have to 
> > tackle this even for FrictPhys_basic.
> >     
>
> That works just fine in the scenario described above, as far as I see.
>
>   
Interesting. For forces, I used exactly the same formulae that you write 
for strain/stress, and indeed it should work just fine in the special 
case described.

if( shearForce.SquaredLength() > maxFs ){
       Real ratio = Mathr::Sqrt(maxFs) / shearForce.Length();
       if (traceEnergy) {//define the plastic work input and increment 
the total plastic energy dissipated
           plasticDissipation +=
               
(shearDisp+(1/currentContactPhysics->ks)*(shearForce-prevForce))//plastic 
disp.
               .Dot(shearForce);//active force
       }


For elastic energy, I don't see any problem(*) in my case :

Real Law2_ScGeom_FrictPhys_Basic::elasticEnergy()
{
   Real energy=0;
   FOREACH(const shared_ptr<Interaction>& I, *scene->interactions){
       if(!I->isReal()) continue;
       FrictPhys* phys = 
dynamic_cast<FrictPhys*>(I->interactionPhysics.get());<---------------------------------- 
see (*)
       if(phys) {
           energy += 0.5*(phys->normalForce.SquaredLength()/phys->kn + 
phys->shearForce.SquaredLength()/phys->ks);}
   }
   return energy;
}

I'm compiling and commiting in minutes. I changed "void updateShear" to 
"Vector3r updateShear", to avoid computing shear displacement twice.
(*) I used to use classIndex to check types, instead of dynamic_casts. 
Is it still possible? What is the best way in your opinion?

Bruno