← Back to team overview

yade-users team mailing list archive

Re: How to get shearstress of particles.

 

Hi Liquing,

I am sorry to not be able to give you a precise answer right now being
busy with externalities; I will get back to you next week (let ne know
if I forget)

 In short, the code should have computed some stress measure on a
particle given interactions with normal and shear forces over them;
after a few minutes of looking at it today, I was not able to see how
was it supposed to work, so I refer you to my next message.

As far as I know, this function was/is only used in VTKRecorder (by
myself) and results looked meaningful, though I give no guarantee. I
will work it out mathematically and document the code (and perhaps fix
it, if necessary) the next week and let you know, as said.

Thanks for asking, you pointed out a weak spot.

Cheers, Vaclav


> I tried to understand the partial codes in shop.cpp. I can't
> understand them. If you know the algorithms, could you tell me? Maybe
> you know where are the related algorithms from, could you tell me the
> information of the related papers. thank you so much. 
> The confusing codes are as follow:
> void Shop::getStressForEachBody(vector<Shop::bodyState>& bodyStates){
>     const shared_ptr<Scene>& scene=Omega::instance().getScene();
>     bodyStates.resize(scene->bodies->size());
>     FOREACH(const shared_ptr<Interaction>& I, *scene->interactions){
>         if(!I->isReal()) continue;
>         const NormShearPhys* phys =
> YADE_CAST<NormShearPhys*>(I->interactionPhysics.get());
>         Dem3DofGeom*
> geom=YADE_CAST<Dem3DofGeom*>(I->interactionGeometry.get());     //For
> the moment only for Dem3DofGeom!!!
>         if(!phys) continue;
>         const body_id_t id1=I->getId1(), id2=I->getId2();
>         
>         Real
> minRad=(geom->refR1<=0?geom->refR2:(geom->refR2<=0?geom->refR1:min(geom->refR1,geom->refR2)));
>         Real crossSection=Mathr::PI*pow(minRad,2);
>         
>         Vector3r
> normalStress=((1./crossSection)*geom->normal.dot(phys->normalForce))*geom->normal;
>         Vector3r shearStress;
>         for(int i=0; i<3; i++){
>             int ix1=(i+1)%3,ix2=(i+2)%3;
>             shearStress[i]=geom-&g
> t;normal[ix1]*phys->shearForce[ix1]+geom->normal[ix2]*phys->shearForce[ix2];
>             shearStress[i]/=crossSection;
>         }
>         
>         bodyStates[id1].normStress+=normalStress;
>         bodyStates[id2].normStress+=normalStress;
>         
>         bodyStates[id1].shearStress+=shearStress;
>         bodyStates[id2].shearStress+=shearStress;
>     }
> }
> How to get the shearstress? I know maybe it is from the t=s-sigma(n),
> but I can't get the coherent expression. I think it should be "
> shearStress[i]=phys->shearForce[i]-(phys->shearForce[i]*geom->normal[i]+geom->normal[ix1]*phys->shearForce[ix1]+geom->normal[ix2]*phys->shea rForce[ix2])*geom->normal[i];". If you have the free time, could you help me to check it? thank you.
> 
> Regards
> Liqing





References