← Back to team overview

yade-users team mailing list archive

Re: Granular ratchetting explained

 

On 14 June 2010 09:22, Bruno Chareyre <bruno.chareyre@xxxxxxxxxxx> wrote:

> Hello Chiara,
>
> Interesting stuff, I agree to put it in trunk when ready.
> 1) I don't know how jumpSe3 and rotation engine work. I modified your
> script to get correct results using NewtonIntegrator and assigning constant
> velocity to the moving body. This script will not work for you unless you
> compile with this modified version of Newton (making non-dynamic bodies move
> at constant speed). At the end of the loop, the residual value of shear
> force is of the order 1e-17*fs.
>
> 2) You cannot compare Law2_ScGeom_FrictPhys_Basic with and without
> ratcheting since ratcheting is always disabled (I anounced that in the list
> before). This line of the script is doing nothing since the functor has no
> attribute "preventGranularRatcheting" (you are just creating a new
> attribute) :
>
Oh I see now! Sorry I was not updated.

>
> contact.preventGranularRatcheting=True
>
> This line has an effect OTOH, but the associated code is not tested (by
> me), I would use False :
>
I tested it by myself, it is giving the right result (I compared it with the
other function rotateAndGetShear since the latter, as you suggested, can be
used also to get total shear displacement, not just total shear force). That
is why I decided to use it.

>
> contact.useShear=True
>
> If you still want to compare ratcheting in Law2_ScGeom_FrictPhys_Basic, a
> quick hack is to change "true"->"false" at line ElasticContactLaw.cpp:77
> Please tell me the result. :)
>
Sure, I let you know.

thank you, chiara

>
> Bruno
>
>
>
>
>
>  I wrote a simple script which might be helpful to test contact laws in
>> order to see if granular ratcheting has any effect, as you suggest below.
>> Please could you have a quick look at the attached script? Try to run it, I
>> just attempted to follow the way you suggest below (the four steps of
>> displacements). There must be something strange in my script because also in
>> the case of preventGranularRatcheting==True I have got a shear force
>> (numerical approximations?) and anyway if I try to change that bool from
>> True to False nothing appears to change.
>>
>> thanks, Chiara
>>
>>
>>
>> On 26 May 2010 08:51, Bruno Chareyre <bruno.chareyre@xxxxxxxxxxx <mailto:
>> bruno.chareyre@xxxxxxxxxxx>> wrote:
>>
>>    For those interested, I elaborated the comment a little in
>>    ScGeom.cpp (possible wiki paragraph in the future), as this
>>    "granular ratchetting" needed explanation.
>>    We could put a simple py script to simulate the cycle explained
>>    below, and test any law in Yade to see if it generates ratchetting.
>>
>>
>>          1. translation "dx" in the normal direction
>>          2. rotation "a"
>>          3. translation "-dx" (back to initial position)
>>          4. rotation "-a" (back to initial orientation)
>>
>>
>> ------------------------------------------------------------------------
>>
>>
>> _______________________________________________
>> Mailing list: https://launchpad.net/~yade-users<https://launchpad.net/%7Eyade-users>
>> Post to     : yade-users@xxxxxxxxxxxxxxxxxxx
>> Unsubscribe : https://launchpad.net/~yade-users<https://launchpad.net/%7Eyade-users>
>> More help   : https://help.launchpad.net/ListHelp
>>
>
>
> --
> _______________
> Bruno Chareyre
> Associate Professor
> Grenoble INP
> Lab. 3SR
> BP 53 - 38041, Grenoble cedex 9 - France
> Tél : 33 4 56 52 86 21
> Fax : 33 4 76 82 70 43
> ________________
>
>
> === modified file 'pkg/dem/Engine/GlobalEngine/NewtonIntegrator.cpp'
> --- pkg/dem/Engine/GlobalEngine/NewtonIntegrator.cpp    2010-06-02 14:19:37
> +0000
> +++ pkg/dem/Engine/GlobalEngine/NewtonIntegrator.cpp    2010-06-14 08:12:23
> +0000
> @@ -96,12 +96,13 @@
>                        State* state=b->state.get();
>                        const body_id_t& id=b->getId();
>                        // clump members are non-dynamic; we only get their
> velocities here
> -                       if (!b->isDynamic || b->isClumpMember()){
> +                       if (/*!b->isDynamic ||*/ b->isClumpMember()){
>                                saveMaximaVelocity(scene,id,state);
>                                continue;
>                        }
>
>                        if (b->isStandalone()){
> +                               if (b->isDynamic){//FIXME : do the same for
> clumps
>                                // translate equation
>                                const Vector3r&
> f=scene->forces.getForce(id);
>                                state->accel=f/state->mass;
> @@ -110,50 +111,57 @@
>                                else {//Apply damping on velocity
> fluctuations only rather than true velocity meanfield+fluctuation.
>                                        Vector3r
> velFluctuation(state->vel-prevVelGrad*state->pos);
>
>  cundallDamp(dt,f,velFluctuation,state->accel);}
> +                               }
>                                leapfrogTranslate(scene,state,id,dt);
>                                // rotate equation: exactAsphericalRot is
> disabled or the body is spherical
>                                if (!exactAsphericalRot ||
> (state->inertia[0]==state->inertia[1] &&
> state->inertia[1]==state->inertia[2])){
> +                                       if (b->isDynamic){//FIXME : do the
> same for clumps
>                                        const Vector3r&
> m=scene->forces.getTorque(id);
>
>  state->angAccel=diagDiv(m,state->inertia);
>
>  cundallDamp(dt,m,state->angVel,state->angAccel);
> +                                       }
>
>  leapfrogSphericalRotate(scene,state,id,dt);
>                                } else { // exactAsphericalRot enabled &
> aspherical body
>                                        const Vector3r&
> m=scene->forces.getTorque(id);
>
>  leapfrogAsphericalRotate(scene,state,id,dt,m);
>                                }
>                        } else if (b->isClump()){
> -                               // clump mass forces
> -                               const Vector3r&
> f=scene->forces.getForce(id);
> -                               Vector3r dLinAccel=f/state->mass;
> -                               cundallDamp(dt,f,state->vel,dLinAccel);
> -                               state->accel+=dLinAccel;
> +                               if (b->isDynamic){// clump mass forces
> +                                       const Vector3r&
> f=scene->forces.getForce(id);
> +                                       Vector3r dLinAccel=f/state->mass;
> +
> cundallDamp(dt,f,state->vel,dLinAccel);
> +                                       state->accel+=dLinAccel;}
>                                const Vector3r&
> m=scene->forces.getTorque(id);
> -                               Vector3r M(m);
> +
>                                // sum force on clump memebrs
>                                // exactAsphericalRot enabled and clump is
> aspherical
>                                if (exactAsphericalRot &&
> ((state->inertia[0]!=state->inertia[1] ||
> state->inertia[1]!=state->inertia[2]))){
> +                                       Vector3r M(m);
>                                        FOREACH(Clump::memberMap::value_type
> mm, static_cast<Clump*>(b.get())->members){
>                                                const body_id_t&
> memberId=mm.first;
>                                                const shared_ptr<Body>&
> member=Body::byId(memberId,scene); assert(member->isClumpMember());
>                                                State*
> memberState=member->state.get();
> -
> handleClumpMemberAccel(scene,memberId,memberState,state);
> -
> handleClumpMemberTorque(scene,memberId,memberState,state,M);
> +                                               if (b->isDynamic){
> +
> handleClumpMemberAccel(scene,memberId,memberState,state);
> +
> handleClumpMemberTorque(scene,memberId,memberState,state,M);}
>
>  saveMaximaVelocity(scene,memberId,memberState);
>                                        }
>                                        // motion
>
>  leapfrogTranslate(scene,state,id,dt);
>
>  leapfrogAsphericalRotate(scene,state,id,dt,M);
>                                } else { // exactAsphericalRot disabled or
> clump is spherical
> -                                       Vector3r
> dAngAccel=diagDiv(M,state->inertia);
> -
> cundallDamp(dt,M,state->angVel,dAngAccel);
> -                                       state->angAccel+=dAngAccel;
> +                                       if (b->isDynamic){
> +                                               Vector3r
> dAngAccel=diagDiv(m,state->inertia);
> +
> cundallDamp(dt,m,state->angVel,dAngAccel);
> +
> state->angAccel+=dAngAccel;}
>                                        FOREACH(Clump::memberMap::value_type
> mm, static_cast<Clump*>(b.get())->members){
>                                                const body_id_t&
> memberId=mm.first;
>                                                const shared_ptr<Body>&
> member=Body::byId(memberId,scene); assert(member->isClumpMember());
>                                                State*
> memberState=member->state.get();
> -
> handleClumpMemberAccel(scene,memberId,memberState,state);
> -
> handleClumpMemberAngAccel(scene,memberId,memberState,state);
> -
> saveMaximaVelocity(scene,memberId,memberState);
> +                                               if (b->isDynamic){
> +
> handleClumpMemberAccel(scene,memberId,memberState,state);
> +
> handleClumpMemberAngAccel(scene,memberId,memberState,state);
> +
> saveMaximaVelocity(scene,memberId,memberState);}
>                                        }
>                                        // motion
>
>  leapfrogTranslate(scene,state,id,dt);
> @@ -181,7 +189,7 @@
>        blockTranslateDOFs(state->blockedDOFs, state->accel);
>
>        if (scene->forces.getMoveRotUsed())
> state->pos+=scene->forces.getMove(id);
> -       if (homotheticCellResize>0) {
> +       if (scene->isPeriodic && homotheticCellResize>0) {
>                // update velocity reflecting changes in the macroscopic
> velocity field, making the problem homothetic.
>                //NOTE : if the velocity is updated before moving the body,
> it means the current velGrad (i.e. before integration in
> cell->integrateAndUpdate) will be effective for the current time-step. Is it
> correct? If not, this velocity update can be moved just after "state->pos +=
> state->vel*dt", meaning the current velocity impulse will be applied at next
> iteration, after the contact law. (All this assuming the ordering is
> resetForces->integrateAndUpdate->contactLaw->PeriCompressor->NewtonsLaw. Any
> other might fool us.)
>                //NOTE : dVel defined without wraping the coordinates means
> bodies out of the (0,0,0) period can move realy fast. It has to be
> compensated properly in the definition of relative velocities (see Ig2
> functors and contact laws).
>
>
> _______________________________________________
> Mailing list: https://launchpad.net/~yade-users<https://launchpad.net/%7Eyade-users>
> Post to     : yade-users@xxxxxxxxxxxxxxxxxxx
> Unsubscribe : https://launchpad.net/~yade-users<https://launchpad.net/%7Eyade-users>
> More help   : https://help.launchpad.net/ListHelp
>
>

Follow ups

References