yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #06523
Re: [Branch ~yade-dev/yade/trunk] Rev 2608: - r2608 again, sorry. Local tests are invalidated by last changes in setDynamic, but the commit i...
Hi Bruno,
sorry could not you spend two words on this last change? I am referring to
the "alpha" parameter you introduced when we avoid granular ratcheting. Many
people (I think) are using this formulation so it would be nice to know what
is causing this change.
Cheers, Chia
On 10 December 2010 21:23, <noreply@xxxxxxxxxxxxx> wrote:
> ------------------------------------------------------------
> revno: 2608
> committer: Bruno Chareyre <bruno.chareyre@xxxxxxxxxxx>
> branch nick: yade
> timestamp: Fri 2010-12-10 21:20:16 +0100
> message:
> - r2608 again, sorry. Local tests are invalidated by last changes in
> setDynamic, but the commit is correct.
> modified:
> pkg/common/Cylinder.cpp
> pkg/dem/ScGeom.cpp
> pkg/dem/ScGeom.hpp
> py/tests/pbc.py
>
>
> --
> lp:yade
> https://code.launchpad.net/~yade-dev/yade/trunk<https://code.launchpad.net/%7Eyade-dev/yade/trunk>
>
> Your team Yade developers is subscribed to branch lp:yade.
> To unsubscribe from this branch go to
> https://code.launchpad.net/~yade-dev/yade/trunk/+edit-subscription<https://code.launchpad.net/%7Eyade-dev/yade/trunk/+edit-subscription>
>
> === modified file 'pkg/common/Cylinder.cpp'
> --- pkg/common/Cylinder.cpp 2010-12-06 19:24:20 +0000
> +++ pkg/common/Cylinder.cpp 2010-12-10 20:20:16 +0000
> @@ -65,6 +65,12 @@
> direction = segment/length;
> dist = direction.dot(branch);
> if ((dist<-interactionDetectionFactor*cylinder->radius) &&
> isNew) return false;
> + if (cylinderSt->rank>0){//make sure there is no contact
> with the previous element in the chain, or consider it a duplicate and
> continue
> + const shared_ptr<Body> cylinderPrev =
> Body::byId(cylinderSt->chains[cylinderSt->chainNumber][cylinderSt->rank-1],scene);
> + Vector3r branchP =
> sphereSt->pos-cylinderPrev->state->pos;
> + Vector3r segmentP = cylinderSt->pos -
> cylinderPrev->state->pos;
> + if (segmentP.dot(branchP)<segmentP.dot(segmentP))
> return false;//will give false on existing interaction, as expected
> + }
> } else {//handle the last node with length=0
> segment = Vector3r(0,0,0);
> branch = sphereSt->pos-cylinderSt->pos-shift2;
>
> === modified file 'pkg/dem/ScGeom.cpp'
> --- pkg/dem/ScGeom.cpp 2010-12-06 19:24:20 +0000
> +++ pkg/dem/ScGeom.cpp 2010-12-10 20:20:16 +0000
> @@ -37,37 +37,30 @@
> }
>
> Vector3r ScGeom::getIncidentVel(const State* rbp1, const State* rbp2, Real
> dt, const Vector3r& shift2, const Vector3r& shiftVel, bool
> avoidGranularRatcheting){
> - Vector3r& x = contactPoint;
> - Vector3r c1x, c2x;
> if(avoidGranularRatcheting){
> //FIXME : put the long comment on the wiki and keep only a
> small abstract and link here.
> - /* B.C. Comment : The following definition of c1x and c2x
> is to avoid "granular ratcheting". This phenomenon has been introduced to me
> by S. McNamara in a seminar held in Paris, 2004 (GDR MiDi). The concept is
> also mentionned in many McNamara, Hermann, Lüding, and co-workers papers
> (see online : "Discrete element methods for the micro-mechanical
> investigation of granular ratcheting", R. García-Rojo, S. McNamara, H. J.
> Herrmann, Proceedings ECCOMAS 2004, @
> http://www.ica1.uni-stuttgart.de/publications/2004/GMH04/), where it
> refers to an accumulation of strain in cyclic loadings.
> -
> - Unfortunately, published papers tend to focus on the "good"
> ratcheting, i.e. irreversible deformations due to the intrinsic nature of
> frictional granular materials, while the talk of McNamara in Paris clearly
> mentioned a possible "bad" ratcheting, purely linked with the formulation of
> the contact laws in what he called "molecular dynamics" (i.e. Cundall's
> model, as opposed to "contact dynamics" from Moreau and Jean).
> -
> - Giving a short explanation :
> - The bad ratcheting is best understood considering this
> small elastic cycle at a contact between two grains : assuming b1 is fixed,
> impose this displacement to b2 :
> + /* B.C. Comment :
> + Giving a short explanation of what we want to avoid :
> + Numerical ratcheting is best understood considering a small
> elastic cycle at a contact between two grains : assuming b1 is fixed, impose
> this displacement to b2 :
> 1. translation "dx" in the normal direction
> 2. rotation "a"
> 3. translation "-dx" (back to initial position)
> 4. rotation "-a" (back to initial orientation)
>
> - If the branch vector used to define the relative shear in
> rotation*branch is not constant (typically if it is defined from the vector
> center->contactPoint like in the "else" below), then the shear displacement
> at the end of this cycle is not null : rotations a and -a are multiplied by
> branches of different lengths.
> - It results in a finite contact force at the end of the
> cycle even though the positions and orientations are unchanged, in total
> contradiction with the elastic nature of the problem. It could also be seen
> as an inconsistent energy creation or loss. It is BAD! And given the fact
> that DEM simulations tend to generate oscillations around equilibrium
> (damped mass-spring), it can have a significant impact on the evolution of
> the packings, resulting for instance in slow creep in iterations under
> constant load.
> - The solution to avoid that is quite simple : use a constant
> branch vector, like here radius_i*normal.
> + This problem has been analyzed in details by McNamara and
> co-workers. One will find interesting discussions in e.g. DOI
> 10.1103/PhysRevE.77.031304, even though solution it suggests is not fully
> applied here (equations of motion are not incorporating alpha, in
> contradiction with what is suggested by McNamara).
> */
> - // FIXME: For sphere-facet contact this will give an
> erroneous value of relative velocity...
> - c1x = radius1*normal;
> - c2x = -radius2*normal;
> + // For sphere-facet contact this will give an erroneous
> value of relative velocity...
> + Real alpha =
> (radius1+radius2)/(radius1+radius2-penetrationDepth);
> + Vector3r relativeVelocity = (rbp2->vel-rbp1->vel)*alpha +
> rbp2->angVel.cross(-radius2*normal) - rbp1->angVel.cross(radius1*normal);
> + relativeVelocity+=alpha*shiftVel;
> + return relativeVelocity;
> } else {
> - // FIXME: It is correct for sphere-sphere and sphere-facet
> contact
> - c1x = (x - rbp1->pos);
> - c2x = (x - rbp2->pos + shift2);
> - }
> - Vector3r relativeVelocity = (rbp2->vel+rbp2->angVel.cross(c2x)) -
> (rbp1->vel+rbp1->angVel.cross(c1x));
> - //update relative velocity for interactions across periods
> - relativeVelocity+=shiftVel;
> - return relativeVelocity;
> + // It is correct for sphere-sphere and sphere-facet contact
> + Vector3r c1x = (contactPoint - rbp1->pos);
> + Vector3r c2x = (contactPoint - rbp2->pos + shift2);
> + Vector3r relativeVelocity =
> (rbp2->vel+rbp2->angVel.cross(c2x)) - (rbp1->vel+rbp1->angVel.cross(c1x));
> + relativeVelocity+=shiftVel;
> + return relativeVelocity;}
> }
>
> Vector3r ScGeom::getIncidentVel(const State* rbp1, const State* rbp2, Real
> dt, bool avoidGranularRatcheting){
> @@ -90,10 +83,8 @@
> if (isNew) {
> initRotations(rbp1,rbp2);
> } else {
> -#ifdef YADE_SCGEOM_DEBUG
> - rbp1.ori.normalize(); rbp2.ori.normalize();
> initialOrientation2.normalize(); initialOrientation1.normalize();
> -#endif
> Quaternionr delta((rbp1.ori *
> (initialOrientation1.conjugate())) * (initialOrientation2 *
> (rbp2.ori.conjugate())));
> + delta.normalize();
> if (creep) delta = delta * twistCreep;
> AngleAxisr aa(delta); // axis of rotation - this is the
> Moment direction UNIT vector; // angle represents the power of resistant
> ELASTIC moment
> //Eigen::AngleAxisr(q) returns nan's when q close to
> identity, next tline fixes the pb.
> @@ -123,10 +114,10 @@
> {
> initialOrientation1 = state1.ori;
> initialOrientation2 = state2.ori;
> -
> initialContactOrientation.setFromTwoVectors(Vector3r(1.0,0.0,0.0),normal);
> - currentContactOrientation = initialContactOrientation;
> - orientationToContact1 = initialOrientation1.conjugate() *
> initialContactOrientation;
> - orientationToContact2 = initialOrientation2.conjugate() *
> initialContactOrientation;
> +//
> initialContactOrientation.setFromTwoVectors(Vector3r(1.0,0.0,0.0),normal);
> +// currentContactOrientation = initialContactOrientation;
> +// orientationToContact1 = initialOrientation1.conjugate() *
> initialContactOrientation;
> +// orientationToContact2 = initialOrientation2.conjugate() *
> initialContactOrientation;
> twist=0;
> bending=Vector3r::Zero();
> twistCreep=Quaternionr(1.0,0.0,0.0,0.0);
>
> === modified file 'pkg/dem/ScGeom.hpp'
> --- pkg/dem/ScGeom.hpp 2010-11-07 11:46:20 +0000
> +++ pkg/dem/ScGeom.hpp 2010-12-10 20:20:16 +0000
> @@ -22,7 +22,7 @@
> Vector3r twist_axis;//rotation vector around normal
> Vector3r orthonormal_axis;//rotation vector in contact plane
> public:
> - // inherited from GenericSpheresContact: Vector3r& normal;
> + // inherited from GenericSpheresContact: Vector3r& normal;
> Real &radius1, &radius2;
> virtual ~ScGeom();
> inline ScGeom& operator= (const ScGeom& source){
> @@ -31,7 +31,7 @@
> radius1=source.radius1; radius2=source.radius2;
> penetrationDepth=source.penetrationDepth;
> shearInc=source.shearInc;
> return *this;}
> -
> +
> //!precompute values of shear increment and interaction
> rotation data. Update contact normal to the currentNormal value.
> Precondition : the value of normal is not updated outside (and before) this
> function.
> void precompute(const State& rbp1, const State& rbp2, const
> Scene* scene, const shared_ptr<Interaction>& c, const Vector3r&
> currentNormal, bool isNew, const Vector3r& shift2, bool
> avoidGranularRatcheting=true);
>
> @@ -43,10 +43,10 @@
> Vector3r getIncidentVel(const State* rbp1, const State*
> rbp2, Real dt, const Vector3r& shiftVel, const Vector3r& shift2, bool
> avoidGranularRatcheting=true);
> // Implement another version of getIncidentVel which does
> not handle periodicity.
> Vector3r getIncidentVel(const State* rbp1, const State*
> rbp2, Real dt, bool avoidGranularRatcheting=true);
> -
> +
> // convenience version to be called from python
> - Vector3r getIncidentVel_py(shared_ptr<Interaction> i, bool
> avoidGranularRatcheting);
> -
> + Vector3r getIncidentVel_py(shared_ptr<Interaction> i, bool
> avoidGranularRatcheting);
> +
>
> YADE_CLASS_BASE_DOC_ATTRS_INIT_CTOR_PY(ScGeom,GenericSpheresContact,"Class
> representing :yref:`geometry<IGeom>` of two :yref:`bodies<Body>` in contact.
> The contact has 3 DOFs (normal and 2×shear) and uses incremental algorithm
> for updating shear. (For shear formulated in total displacements and
> rotations, see :yref:`Dem3DofGeom` and related classes).\n\nWe use symbols
> $\\vec{x}$, $\\vec{v}$, $\\vec{\\omega}$ respectively for position, linear
> and angular velocities (all in global coordinates) and $r$ for particles
> radii; subscripted with 1 or 2 to distinguish 2 spheres in contact. Then we
> compute unit contact normal\n\n..
> math::\n\n\t\\vec{n}=\\frac{\\vec{x}_2-\\vec{x}_1}{||\\vec{x}_2-\\vec{x}_1||}\n\nRelative
> velocity of spheres is then\n\n..
> math::\n\n\t\\vec{v}_{12}=(\\vec{v}_2+\\vec{\\omega}_2\\times(-r_2\\vec{n}))-(\\vec{v}_1+\\vec{\\omega}_1\\times(r_1\\vec{n}))\n\nand
> its shear component\n\n..
> math::\n\n\t\\Delta\\vec{v}_{12}^s=\\vec{v}_{12}-(\\vec{n}\\cdot\\vec{v}_{12})\\vec{n}.\n\nTangential
> displacement increment over last step then reads\n\n..
> math::\n\n\t\\vec{x}_{12}^s=\\Delta t \\vec{v}_{12}^s.",
>
> ((Real,penetrationDepth,NaN,(Attr::noSave|Attr::readonly),"Penetration
> distance of spheres (positive if overlapping)"))
>
> ((Vector3r,shearInc,Vector3r::Zero(),(Attr::noSave|Attr::readonly),"Shear
> displacement increment in the last step"))
> @@ -67,15 +67,15 @@
>
> void precomputeRotations(const State& rbp1, const State&
> rbp2, bool isNew, bool creep=false);
> void initRotations(const State& rbp1, const State& rbp2);
> -
> +
>
> YADE_CLASS_BASE_DOC_ATTRS_INIT_CTOR_PY(ScGeom6D,ScGeom,"Class representing
> :yref:`geometry<IGeom>` of two :yref:`bodies<Body>` in contact. The contact
> has 6 DOFs (normal, 2×shear, twist, 2xbending) and uses :yref:`ScGeom`
> incremental algorithm for updating shear.",
> -
> ((Quaternionr,orientationToContact1,Quaternionr(1.0,0.0,0.0,0.0),,""))
> -
> ((Quaternionr,orientationToContact2,Quaternionr(1.0,0.0,0.0,0.0),,""))
> +//
> ((Quaternionr,orientationToContact1,Quaternionr(1.0,0.0,0.0,0.0),,""))
> +//
> ((Quaternionr,orientationToContact2,Quaternionr(1.0,0.0,0.0,0.0),,""))
>
> ((Quaternionr,initialOrientation1,Quaternionr(1.0,0.0,0.0,0.0),,""))
>
> ((Quaternionr,initialOrientation2,Quaternionr(1.0,0.0,0.0,0.0),,""))
> ((Quaternionr,twistCreep,Quaternionr(1.0,0.0,0.0,0.0),,""))
> -
> ((Quaternionr,currentContactOrientation,Quaternionr(1.0,0.0,0.0,0.0),,""))
> -
> ((Quaternionr,initialContactOrientation,Quaternionr(1.0,0.0,0.0,0.0),,""))
> +//
> ((Quaternionr,currentContactOrientation,Quaternionr(1.0,0.0,0.0,0.0),,""))
> +//
> ((Quaternionr,initialContactOrientation,Quaternionr(1.0,0.0,0.0,0.0),,""))
> ((Real,twist,0,(Attr::noSave | Attr::readonly),"Elastic
> twist angle of the contact."))
> ((Vector3r,bending,Vector3r::Zero(),(Attr::noSave |
> Attr::readonly),"Bending at contact as a vector defining axis of rotation
> and angle (angle=norm)."))
> ,
>
> === modified file 'py/tests/pbc.py'
> --- py/tests/pbc.py 2010-10-26 13:41:30 +0000
> +++ py/tests/pbc.py 2010-12-10 20:20:16 +0000
> @@ -20,7 +20,7 @@
> O.reset(); O.periodic=True;
> O.cell.refSize=Vector3(2.5,2.5,3)
> self.cellDist=Vector3i(0,0,10) # how many cells away we go
> - self.relDist=Vector3(0,.999,0) # rel position of the 2nd
> ball within the cell
> + self.relDist=Vector3(0,.999999999999999999,0) # rel
> position of the 2nd ball within the cell
> self.initVel=Vector3(0,0,5)
> O.bodies.append(utils.sphere((1,1,1),.5))
>
> self.initPos=Vector3([O.bodies[0].state.pos[i]+self.relDist[i]+self.cellDist[i]*O.cell.refSize[i]
> for i in (0,1,2)])
>
>
> _______________________________________________
> Mailing list: https://launchpad.net/~yade-dev<https://launchpad.net/%7Eyade-dev>
> Post to : yade-dev@xxxxxxxxxxxxxxxxxxx
> Unsubscribe : https://launchpad.net/~yade-dev<https://launchpad.net/%7Eyade-dev>
> More help : https://help.launchpad.net/ListHelp
>
>
References
