yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #09996
RE : RE : Modified GlobalStiffnessTimeStepper for visco-elastic contact law
Here it is!
Thanks Anton
Raphael
-------- Message d'origine--------
De: Anton Gladky [mailto:gladky.anton@xxxxxxxxx]
Date: lun. 23/09/2013 16:47
À: Raphaël Maurin
Cc: yade-dev
Objet : Re: [Yade-dev] RE : Modified GlobalStiffnessTimeStepper for visco-elastic contact law
Hi Raphaël,
could you, please, attach a diff (or patch) instead of copying it into
the mail-body.
Thanks,
Anton
2013/9/20 Raphaël Maurin <raphael.maurin@xxxxxxxxx>:
> Hi all,
> Here is the final version of GlobalStiffnessTimeStepper that I modified to
> take into account the viscous contribution when using
> Law2_ScGeom_ViscElPhys_Basic.
> I didn't modify the elastic part, so that it will work exactly the same for
> people using it as before (considering only elastic contribution). I just
> added an option : viscEl which, when True, allows to calculate the time step
> for ViscElPhys as described before.
> I tested it and everything works well.
>
> I put the patch below for GlobalStiffnessTimeStepper.cpp and
> GlobalStiffnessTimeStepper.hpp, can someone add the change on github ?
> Thanks
>
> Raphael
>
>
>
> GlobalStiffnessTimeStepper.cpp
>
>
>
> diff --git a/pkg/dem/GlobalStiffnessTimeStepper.cpp
> b/pkg/dem/GlobalStiffnessTimeStepper.cpp
> index 8a41ceb..0205655 100644
> --- a/pkg/dem/GlobalStiffnessTimeStepper.cpp
> +++ b/pkg/dem/GlobalStiffnessTimeStepper.cpp
> @@ -14,6 +14,7 @@
> #include<yade/core/Scene.hpp>
> #include<yade/core/Clump.hpp>
> #include<yade/pkg/dem/Shop.hpp>
> +#include<yade/pkg/dem/ViscoelasticPM.hpp>
>
> CREATE_LOGGER(GlobalStiffnessTimeStepper);
> YADE_PLUGIN((GlobalStiffnessTimeStepper));
> @@ -25,13 +26,16 @@ void
> GlobalStiffnessTimeStepper::findTimeStepFromBody(const shared_ptr<Body>& bo
> State* sdec=body->state.get();
> Vector3r& stiffness= stiffnesses[body->getId()];
> Vector3r& Rstiffness=Rstiffnesses[body->getId()];
> -
> if(body->isClump()) {// if clump, we sum stifnesses of all members
> const shared_ptr<Clump>&
> clump=YADE_PTR_CAST<Clump>(body->shape);
> FOREACH(Clump::MemberMap::value_type& B, clump->members){
> const shared_ptr<Body>& b =
> Body::byId(B.first,scene);
> stiffness+=stiffnesses[b->getId()];
> Rstiffness+=Rstiffnesses[b->getId()];
> + if (viscEl == true){
> +
> viscosities[body->getId()]+=viscosities[b->getId()];
> +
> Rviscosities[body->getId()]+=Rviscosities[b->getId()];
> + }
> }
> }
>
> @@ -39,13 +43,15 @@ void
> GlobalStiffnessTimeStepper::findTimeStepFromBody(const shared_ptr<Body>& bo
> if (densityScaling) sdec->densityScaling =
> min(1.0001*sdec->densityScaling,
> timestepSafetyCoefficient*pow(defaultDt/targetDt,2.0));
> return; // not possible to compute!
> }
> -
> +
> + //Determine the elastic minimum eigenperiod (and if required
> determine also the viscous one separately and take the minimum of the two)
> +
> + //Elastic
> Real dtx, dty, dtz;
> Real dt = max( max (stiffness.x(), stiffness.y()), stiffness.z() );
> if (dt!=0) {
> dt = sdec->mass/dt; computedSomething = true;}//dt =
> squared eigenperiod of translational motion
> else dt = Mathr::MAX_REAL;
> -
> if (Rstiffness.x()!=0) {
> dtx = sdec->inertia.x()/Rstiffness.x(); computedSomething =
> true;}//dtx = squared eigenperiod of rotational motion around x
> else dtx = Mathr::MAX_REAL;
> @@ -56,15 +62,43 @@ void
> GlobalStiffnessTimeStepper::findTimeStepFromBody(const shared_ptr<Body>& bo
> dtz = sdec->inertia.z()/Rstiffness.z(); computedSomething =
> true;}
> else dtz = Mathr::MAX_REAL;
>
> - Real Rdt = std::min( std::min (dtx, dty), dtz );//Rdt = smallest
> squared eigenperiod for rotational motions
> + Real Rdt = std::min( std::min (dtx, dty), dtz );//Rdt = smallest
> squared eigenperiod for elastic rotational motions
> dt =
> 1.41044*timestepSafetyCoefficient*std::sqrt(std::min(dt,Rdt));//1.41044 =
> sqrt(2)
> +
> + //Viscous
> + if (viscEl == true){
> + Vector3r& viscosity = viscosities[body->getId()];
> + Vector3r& Rviscosity = Rviscosities[body->getId()];
> + Real dtx_visc, dty_visc, dtz_visc;
> + Real dt_visc = max(max(viscosity.x(), viscosity.y()),
> viscosity.z() );
> + if (dt_visc!=0) {
> + dt_visc = sdec->mass/dt_visc; computedSomething =
> true;}//dt = eigenperiod of the viscous translational motion
> + else {dt_visc = Mathr::MAX_REAL;}
> +
> + if (Rviscosity.x()!=0) {
> + dtx_visc = sdec->inertia.x()/Rviscosity.x();
> computedSomething = true;}//dtx = eigenperiod of viscous rotational motion
> around x
> + else dtx_visc = Mathr::MAX_REAL;
> + if (Rviscosity.y()!=0) {
> + dty_visc = sdec->inertia.y()/Rviscosity.y();
> computedSomething = true;}
> + else dty_visc = Mathr::MAX_REAL;
> + if (Rviscosity.z()!=0) {
> + dtz_visc = sdec->inertia.z()/Rviscosity.z();
> computedSomething = true;}
> + else dtz_visc = Mathr::MAX_REAL;
> +
> + Real Rdt_visc = std::min( std::min (dtx_visc, dty_visc),
> dtz_visc );//Rdt = smallest squared eigenperiod for viscous rotational
> motions
> + dt_visc =
> 2*timestepSafetyCoefficient*std::min(dt_visc,Rdt_visc);
> +
> + //Take the minimum between the elastic and viscous minimum
> eigenperiod.
> + dt = std::min(dt,dt_visc);
> + }
> +
> //if there is a target dt, then we apply density scaling on the
> body, the inertia used in Newton will be mass*scaling, the weight is
> unmodified
> if (densityScaling) {
> sdec->densityScaling =
> min(sdec->densityScaling,timestepSafetyCoefficient*pow(dt /targetDt,2.0));
> newDt=targetDt;
> }
> //else we update dt normaly
> - else {newDt = std::min(dt,newDt);}
> + else {newDt = std::min(dt,newDt);}
> }
>
> bool GlobalStiffnessTimeStepper::isActivated()
> @@ -97,7 +131,7 @@ void GlobalStiffnessTimeStepper::computeTimeStep(Scene*
> ncb)
> computedOnce = true;}
> else if (!computedOnce) scene->dt=defaultDt;
> LOG_INFO("computed timestep " << newDt <<
> - (scene->dt==newDt ? string(", appplied") :
> + (scene->dt==newDt ? string(", applied") :
> string(", BUT timestep is
> ")+lexical_cast<string>(scene->dt))<<".");
> }
>
> @@ -107,15 +141,24 @@ void
> GlobalStiffnessTimeStepper::computeStiffnesses(Scene* rb){
> if(size<rb->bodies->size()){
> size=rb->bodies->size();
> stiffnesses.resize(size); Rstiffnesses.resize(size);
> + if (viscEl == true){
> + viscosities.resize(size); Rviscosities.resize(size);
> + }
> }
> /* reset stored values */
> memset(& stiffnesses[0],0,sizeof(Vector3r)*size);
> memset(&Rstiffnesses[0],0,sizeof(Vector3r)*size);
> + if (viscEl == true){
> + memset(& viscosities[0],0,sizeof(Vector3r)*size);
> + memset(&Rviscosities[0],0,sizeof(Vector3r)*size);
> + }
> +
> FOREACH(const shared_ptr<Interaction>& contact, *rb->interactions){
> if(!contact->isReal()) continue;
>
> GenericSpheresContact*
> geom=YADE_CAST<GenericSpheresContact*>(contact->geom.get()); assert(geom);
> NormShearPhys*
> phys=YADE_CAST<NormShearPhys*>(contact->phys.get()); assert(phys);
> +
> // all we need for getting stiffness
> Vector3r& normal=geom->normal; Real& kn=phys->kn; Real&
> ks=phys->ks; Real& radius1=geom->refR1; Real& radius2=geom->refR2;
> Real fn = (static_cast<NormShearPhys *>
> (contact->phys.get()))->normalForce.squaredNorm();
> @@ -134,10 +177,35 @@ void
> GlobalStiffnessTimeStepper::computeStiffnesses(Scene* rb){
>
> std::pow(normal.x(),2)+std::pow(normal.z(),2),
>
> std::pow(normal.x(),2)+std::pow(normal.y(),2));
> diag_Rstiffness *= ks;
> +
> +
> //NOTE : contact laws with moments would be handled
> correctly by summing directly bending+twisting stiffness to diag_Rstiffness.
> The fact that there is no problem currently with e.g. cohesiveFrict law is
> probably because final computed dt is constrained by translational motion,
> not rotations.
> stiffnesses [contact->getId1()]+=diag_stiffness;
>
> Rstiffnesses[contact->getId1()]+=diag_Rstiffness*pow(radius1,2);
> stiffnesses [contact->getId2()]+=diag_stiffness;
> -
> Rstiffnesses[contact->getId2()]+=diag_Rstiffness*pow(radius2,2);
> +
> Rstiffnesses[contact->getId2()]+=diag_Rstiffness*pow(radius2,2);
> +
> + //Same for the Viscous part, if required
> + if (viscEl == true){
> + ViscElPhys* viscPhys =
> YADE_CAST<ViscElPhys*>(contact->phys.get()); assert(viscPhys);
> + Real& cn=viscPhys->cn; Real& cs=viscPhys->cs;
> + //Diagonal terms of the translational viscous matrix
> + Vector3r diag_viscosity =
> Vector3r(std::pow(normal.x(),2),std::pow(normal.y(),2),std::pow(normal.z(),2));
> + diag_viscosity *= cn-cs;
> + diag_viscosity = diag_viscosity +
> Vector3r(1,1,1)*cs;
> + //diagonal terms of the rotational viscous matrix
> + Vector3r diag_Rviscosity =
> +
> Vector3r(std::pow(normal.y(),2)+std::pow(normal.z(),2),
> +
> std::pow(normal.x(),2)+std::pow(normal.z(),2),
> +
> std::pow(normal.x(),2)+std::pow(normal.y(),2));
> + diag_Rviscosity *= cs;
> +
> + // Add the contact stiffness matrix to the two
> particles one
> + viscosities [contact->getId1()]+=diag_viscosity;
> +
> Rviscosities[contact->getId1()]+=diag_Rviscosity*pow(radius1,2);
> + viscosities [contact->getId2()]+=diag_viscosity;
> +
> Rviscosities[contact->getId2()]+=diag_Rviscosity*pow(radius2,2);
> + }
> +
> }
> }
>
>
>
>
>
>
> GlobalStiffnessTimeStepper.hpp
>
>
>
>
>
>
> diff --git a/pkg/dem/GlobalStiffnessTimeStepper.hpp
> b/pkg/dem/GlobalStiffnessTimeStepper.hpp
> index 1d60f88..bd657fa 100644
> --- a/pkg/dem/GlobalStiffnessTimeStepper.hpp
> +++ b/pkg/dem/GlobalStiffnessTimeStepper.hpp
> @@ -24,6 +24,8 @@ class GlobalStiffnessTimeStepper : public TimeStepper
> private :
> vector<Vector3r> stiffnesses;
> vector<Vector3r> Rstiffnesses;
> + vector<Vector3r> viscosities;
> + vector<Vector3r> Rviscosities;
> void computeStiffnesses(Scene*);
>
> Real newDt;
> @@ -37,13 +39,14 @@ class GlobalStiffnessTimeStepper : public TimeStepper
> virtual void computeTimeStep(Scene*);
> virtual bool isActivated();
> YADE_CLASS_BASE_DOC_ATTRS_CTOR(
> - GlobalStiffnessTimeStepper,TimeStepper,"An engine
> assigning the time-step as a fraction of the minimum eigen-period in the
> problem. The derivation is detailed in the chapter on `DEM formulation
> <https://www.yade-dem.org/doc/formulation.html#dem-simulations>`_",
> + GlobalStiffnessTimeStepper,TimeStepper,"An engine
> assigning the time-step as a fraction of the minimum eigen-period in the
> problem. The derivation is detailed in the chapter on `DEM formulation
> <https://www.yade-dem.org/doc/formulation.html#dem-simulations>`_. The
> viscEl option enables to evaluate the timestep in a similar way for the
> visco-elastic contact law :yref:`Law2_ScGeom_ViscElPhys_Basic`, more detail
> in :yref:`GlobalStiffnessTimestepper::viscEl`. ",
> ((Real,defaultDt,-1,,"used as the initial value of
> the timestep (especially useful in the first steps when no contact exist).
> If negative, it will be defined by :yref:`utils.PWaveTimeStep` *
> :yref:`GlobalStiffnessTimeStepper::timestepSafetyCoefficient`"))
> ((Real,maxDt,Mathr::MAX_REAL,,"if positive, used as
> max value of the timestep whatever the computed value"))
> ((Real,previousDt,1,,"last computed dt |yupdate|"))
> ((Real,timestepSafetyCoefficient,0.8,,"safety factor
> between the minimum eigen-period and the final assigned dt (less than 1))"))
> ((bool,densityScaling,false,,"|yupdate| don't modify
> this value if you don't plan to modify the scaling factor manually for some
> bodies. In most cases, it is enough to set
> :yref:`NewtonIntegrator::densityScaling` and let this one be adjusted
> automatically."))
> - ((Real,targetDt,1,,"if
> :yref:`NewtonIntegrator::densityScaling` is active, this value will be used
> as the simulation timestep and the scaling will use this value of dt as the
> target value. The value of targetDt is arbitrary and should have no effect
> in the result in general. However if some bodies have imposed velocities,
> for instance, they will move more or less per each step depending on this
> value.")),
> + ((Real,targetDt,1,,"if
> :yref:`NewtonIntegrator::densityScaling` is active, this value will be used
> as the simulation timestep and the scaling will use this value of dt as the
> target value. The value of targetDt is arbitrary and should have no effect
> in the result in general. However if some bodies have imposed velocities,
> for instance, they will move more or less per each step depending on this
> value."))
> + ((bool,viscEl,false,,"To use with
> :yref:`ViscElPhys`. if True, evaluate separetly the minimum eigen-period in
> the problem considering only the elastic contribution on one hand (spring
> only), and only the viscous contribution on the other hand (dashpot only).
> Take then the minimum of the two and use the safety coefficient
> :yref:`GlobalStiffnessTimestepper::timestepSafetyCoefficient` to take into
> account the possible coupling between the two contribution.")),
> computedOnce=false;)
> DECLARE_LOGGER;
> };
>
>
>
>
>
>
>
>
>
>
> -------- Message d'origine--------
> De: Raphaël Maurin
> Date: mar. 03/09/2013 12:21
> À: yade-dev@xxxxxxxxxxxxxxxxxxx
> Objet : Modified GlobalStiffnessTimeStepper for visco-elastic contact law
>
> Hi all,
>
> I developed a "copy" of GlobalStiffnessTimeStepper which evaluate the time
> step of a system of particle with Law2_ScGeom_ViscElPhys_Basic contact law
> (spring + viscous damping).
>
> The principle of the time step evaluation is :
> We consider that the elastic and the viscous degrees of freedom are
> decoupled. From that, we evaluate the critical time step of the system of
> particles considering the elastic part only (dt_el) on one hand, and
> considering the viscous part only (dt_visc) on the other hand. Then, we take
> the minimum of the two timestep and multiply it by a safety factor to
> account for the possible coupling between the two effect (dt =
> alpha*min(dt_el,dt_visc).
>
> To do so, I took GlobalStiffnessTimeStepper (GSTS) which is already
> evaluating the elastic critical time step, and implemented similarly a
> calculation for the viscous critical time step of the system from the the
> viscosity matrix of each particles. The changes are not so important between
> the version I implemented and the last one (see the code below). I tested it
> in cases where the normal viscous part is dominant and the new version gives
> a lower time step which allows to keep the caculus stable, while in the old
> version everything explode, so it seems to work.
>
> I have some questions about that :
> - Does the evaluation of the time step seems reasonable to you?
> - What can I do to test it quantitavely ? (if I try a collision between two
> beads, GSTS does not work..)
> - To add it to the code, should I add it as a new "file"/"function" or
> should I implement it as an option of GSTS ?
>
> And about GSTS in general :
> - why is GSTS not working when we consider a gravity deposition or a
> collision ?? (even when defaultDt is small enough)
>
> Thanks,
>
> Raphael
>
>
>
> GlobalStiffnessTimeStepper.cpp
>
> /*************************************************************************
> * Copyright (C) 2006 by Bruno Chareyre *
> * bruno.chareyre@xxxxxxxxxxx *
> * *
> * This program is free software; it is licensed under the terms of the *
> * GNU General Public License v2 or later. See file LICENSE for details. *
> *************************************************************************/
>
> #include"GlobalStiffnessTimeStepper.hpp"
> #include<yade/pkg/dem/FrictPhys.hpp>
> #include<yade/pkg/dem/ScGeom.hpp>
> #include<yade/pkg/dem/DemXDofGeom.hpp>
> #include<yade/core/Interaction.hpp>
> #include<yade/core/Scene.hpp>
> #include<yade/core/Clump.hpp>
> #include<yade/pkg/dem/Shop.hpp>
> #include<yade/pkg/dem/ViscoelasticPM.hpp>
>
> CREATE_LOGGER(GlobalStiffnessTimeStepper);
> YADE_PLUGIN((GlobalStiffnessTimeStepper));
>
> GlobalStiffnessTimeStepper::~GlobalStiffnessTimeStepper() {}
>
> void GlobalStiffnessTimeStepper::findTimeStepFromBody(const
> shared_ptr<Body>& body, Scene * ncb)
> {
> State* sdec=body->state.get();
> Vector3r& stiffness= stiffnesses[body->getId()];
> Vector3r& Rstiffness=Rstiffnesses[body->getId()];
> Vector3r& viscosity= viscosities[body->getId()];
> Vector3r& Rviscosity= Rviscosities[body->getId()];
>
> if(body->isClump()) {// if clump, we sum stifnesses of all members
> const shared_ptr<Clump>&
> clump=YADE_PTR_CAST<Clump>(body->shape);
> FOREACH(Clump::MemberMap::value_type& B, clump->members){
> const shared_ptr<Body>& b =
> Body::byId(B.first,scene);
> stiffness+=stiffnesses[b->getId()];
> Rstiffness+=Rstiffnesses[b->getId()];
> viscosity+=viscosities[b->getId()];
> Rviscosity+=Rviscosities[b->getId()];
> }
> }
>
> if(!sdec || stiffness==Vector3r::Zero()){
> if (densityScaling) sdec->densityScaling =
> min(1.0001*sdec->densityScaling,
> timestepSafetyCoefficient*pow(defaultDt/targetDt,2.0));
> return; // not possible to compute!
> }
>
> //Determine the elastic critical time step and the viscous one
> separately (as if they were decoupled) and take the minimum time step
> with a safety coefficient.
>
> //Elastic
> Real dtx, dty, dtz;
> Real dt = max( max (stiffness.x(), stiffness.y()), stiffness.z() );
> if (dt!=0) {
> dt = sdec->mass/dt; computedSomething = true;}//dt =
> squared eigenperiod of translational motion
> else dt = Mathr::MAX_REAL;
> if (Rstiffness.x()!=0) {
> dtx = sdec->inertia.x()/Rstiffness.x(); computedSomething =
> true;}//dtx = squared eigenperiod of rotational motion around x
> else dtx = Mathr::MAX_REAL;
> if (Rstiffness.y()!=0) {
> dty = sdec->inertia.y()/Rstiffness.y(); computedSomething =
> true;}
> else dty = Mathr::MAX_REAL;
> if (Rstiffness.z()!=0) {
> dtz = sdec->inertia.z()/Rstiffness.z(); computedSomething =
> true;}
> else dtz = Mathr::MAX_REAL;
>
> Real Rdt = std::min( std::min (dtx, dty), dtz );//Rdt = smallest
> squared eigenperiod for elastic rotational motions
> dt =
> 1.41044*timestepSafetyCoefficient*std::sqrt(std::min(dt,Rdt));//1.41044 =
> sqrt(2)
>
> //Viscous
> Real dtx_visc, dty_visc, dtz_visc;
> Real dt_visc = max(max(viscosity.x(), viscosity.y()), viscosity.z()
> );
> if (dt_visc!=0) {
> dt_visc = sdec->mass/dt_visc; computedSomething =
> true;}//dt = eigenperiod of the viscous translational motion
> else {dt_visc = Mathr::MAX_REAL;}
>
> if (Rviscosity.x()!=0) {
> dtx_visc = sdec->inertia.x()/Rviscosity.x();
> computedSomething = true;}//dtx = eigenperiod of viscous rotational motion
> around x
> else dtx_visc = Mathr::MAX_REAL;
> if (Rviscosity.y()!=0) {
> dty_visc = sdec->inertia.y()/Rviscosity.y();
> computedSomething = true;}
> else dty_visc = Mathr::MAX_REAL;
> if (Rviscosity.z()!=0) {
> dtz_visc = sdec->inertia.z()/Rviscosity.z();
> computedSomething = true;}
> else dtz_visc = Mathr::MAX_REAL;
>
> Real Rdt_visc = std::min( std::min (dtx_visc, dty_visc), dtz_visc
> );//Rdt = smallest squared eigenperiod for viscous rotational motions
> dt_visc = 2*timestepSafetyCoefficient*std::min(dt_visc,Rdt_visc);
>
> //Take the minimum between the elastic and viscous critical time
> step.
> dt = std::min(dt,dt_visc);
>
> //if there is a target dt, then we apply density scaling on the
> body, the inertia used in Newton will be mass*scaling, the weight is
> unmodified
> if (densityScaling) {
> sdec->densityScaling =
> min(sdec->densityScaling,timestepSafetyCoefficient*pow(dt /targetDt,2.0));
> newDt=targetDt;
> }
> //else we update dt normaly
> else {newDt = std::min(dt,newDt);}
> }
>
> bool GlobalStiffnessTimeStepper::isActivated()
> {
> return (active && ((!computedOnce) || (scene->iter %
> timeStepUpdateInterval == 0) || (scene->iter < (long int) 2) ));
> }
>
> void GlobalStiffnessTimeStepper::computeTimeStep(Scene* ncb)
> {
> // for some reason, this line is necessary to have correct
> functioning (no idea _why_)
> // see scripts/test/compare-identical.py, run with or without
> active=active.
> active=active;
> if (defaultDt<0) defaultDt=
> timestepSafetyCoefficient*Shop::PWaveTimeStep(Omega::instance().getScene());
> computeStiffnesses(ncb);
>
> shared_ptr<BodyContainer>& bodies = ncb->bodies;
> newDt = Mathr::MAX_REAL;
> computedSomething=false;
> BodyContainer::iterator bi = bodies->begin();
> BodyContainer::iterator biEnd = bodies->end();
> for( ; bi!=biEnd ; ++bi ){
> shared_ptr<Body> b = *bi;
> if (b->isDynamic() && !b->isClumpMember())
> findTimeStepFromBody(b, ncb);
>
> }
> if(densityScaling) (newDt=targetDt);
> if(computedSomething || densityScaling){
> previousDt = min ( min(newDt , maxDt), 1.05*previousDt );//
> at maximum, dt will be multiplied by 1.05 in one iterration, this is to
> prevent brutal switches from 0.000... to 1 in some computations
> scene->dt=previousDt;
> computedOnce = true;}
> else if (!computedOnce) scene->dt=defaultDt;
> LOG_INFO("computed timestep " << newDt <<
> (scene->dt==newDt ? string(", applied") :
> string(", BUT timestep is
> ")+lexical_cast<string>(scene->dt))<<".");
> }
>
> void GlobalStiffnessTimeStepper::computeStiffnesses(Scene* rb){
> /* check size */
> size_t size=stiffnesses.size();
> if(size<rb->bodies->size()){
> size=rb->bodies->size();
> stiffnesses.resize(size); Rstiffnesses.resize(size);
> viscosities.resize(size); Rviscosities.resize(size);
> }
> /* reset stored values */
> memset(& stiffnesses[0],0,sizeof(Vector3r)*size);
> memset(&Rstiffnesses[0],0,sizeof(Vector3r)*size);
> memset(& viscosities[0],0,sizeof(Vector3r)*size);
> memset(& Rviscosities[0],0,sizeof(Vector3r)*size);
>
> FOREACH(const shared_ptr<Interaction>& contact, *rb->interactions){
> if(!contact->isReal()) continue;
>
> GenericSpheresContact*
> geom=YADE_CAST<GenericSpheresContact*>(contact->geom.get()); assert(geom);
> ViscElPhys* phys =
> YADE_CAST<ViscElPhys*>(contact->phys.get()); assert(phys);
> // all we need for getting stiffness
> Vector3r& normal=geom->normal; Real& kn=phys->kn; Real&
> ks=phys->ks; Real& radius1=geom->refR1; Real& radius2=geom->refR2;
> Real& cn=phys->cn; Real& cs=phys->cs;
> Real fn = (static_cast<NormShearPhys *>
> (contact->phys.get()))->normalForce.squaredNorm();
> if (fn==0) continue;//Is it a problem with some laws? I
> still don't see why.
>
> //Diagonal terms of the translational stiffness matrix
> Vector3r diag_stiffness =
> Vector3r(std::pow(normal.x(),2),std::pow(normal.y(),2),std::pow(normal.z(),2));
> diag_stiffness *= kn-ks;
> diag_stiffness = diag_stiffness + Vector3r(1,1,1)*ks;
>
> //diagonal terms of the rotational stiffness matrix
> // Vector3r branch1 =
> currentContactGeometry->normal*currentContactGeometry->radius1;
> // Vector3r branch2 =
> currentContactGeometry->normal*currentContactGeometry->radius2;
> Vector3r diag_Rstiffness =
>
> Vector3r(std::pow(normal.y(),2)+std::pow(normal.z(),2),
>
> std::pow(normal.x(),2)+std::pow(normal.z(),2),
>
> std::pow(normal.x(),2)+std::pow(normal.y(),2));
> diag_Rstiffness *= ks;
>
> //Diagonal terms of the translational viscous matrix
> Vector3r diag_viscosity =
> Vector3r(std::pow(normal.x(),2),std::pow(normal.y(),2),std::pow(normal.z(),2));
> diag_viscosity *= cn-cs;
> diag_viscosity = diag_viscosity + Vector3r(1,1,1)*cs;
> //diagonal terms of the rotational viscous matrix
> Vector3r diag_Rviscosity =
>
> Vector3r(std::pow(normal.y(),2)+std::pow(normal.z(),2),
>
> std::pow(normal.x(),2)+std::pow(normal.z(),2),
>
> std::pow(normal.x(),2)+std::pow(normal.y(),2));
> diag_Rviscosity *= cs;
>
> //NOTE : contact laws with moments would be handled
> correctly by summing directly bending+twisting stiffness to diag_Rstiffness.
> The fact that there is no problem currently with e.g. cohesiveFrict law is
> probably because final computed dt is constrained by translational motion,
> not rotations.
> stiffnesses [contact->getId1()]+=diag_stiffness;
>
> Rstiffnesses[contact->getId1()]+=diag_Rstiffness*pow(radius1,2);
> viscosities [contact->getId1()]+=diag_viscosity;
>
> Rviscosities[contact->getId1()]+=diag_Rviscosity*pow(radius1,2);
> stiffnesses [contact->getId2()]+=diag_stiffness;
>
> Rstiffnesses[contact->getId2()]+=diag_Rstiffness*pow(radius2,2);
> viscosities [contact->getId2()]+=diag_viscosity;
>
> Rviscosities[contact->getId2()]+=diag_Rviscosity*pow(radius2,2);
> }
> }
>
>
>
>
>
>
>
>
>
>
>
> GlobalStiffnessTimeStepper.hpp
>
>
>
> /*************************************************************************
> * Copyright (C) 2006 by Bruno Chareyre *
> * bruno.chareyre@xxxxxxxxxxx *
> * *
> * This program is free software; it is licensed under the terms of the *
> * GNU General Public License v2 or later. See file LICENSE for details. *
> *************************************************************************/
>
> #pragma once
>
> #include<yade/core/TimeStepper.hpp>
>
>
> /*! \brief Compute the critical timestep of the leap-frog scheme based on
> global stiffness of bodies
> See usage details in TriaxialTest documentation (TriaxialTest is
> also a good example of how to use this class)
> */
>
> class Interaction;
> class BodyContainer;
> class Scene;
>
> class GlobalStiffnessTimeStepper : public TimeStepper
> {
> private :
> vector<Vector3r> stiffnesses;
> vector<Vector3r> Rstiffnesses;
> vector<Vector3r> viscosities;
> vector<Vector3r> Rviscosities;
> void computeStiffnesses(Scene*);
>
> Real newDt;
> bool computedSomething,
> computedOnce;
> void findTimeStepFromBody(const shared_ptr<Body>& body,
> Scene * ncb);
>
> public :
> virtual ~GlobalStiffnessTimeStepper();
>
> virtual void computeTimeStep(Scene*);
> virtual bool isActivated();
> YADE_CLASS_BASE_DOC_ATTRS_CTOR(
> GlobalStiffnessTimeStepper,TimeStepper,"An engine
> assigning the time-step as a fraction of the minimum eigen-period in the
> problem. The derivation is detailed in the chapter on `DEM formulation
> <https://www.yade-dem.org/doc/formulation.html#dem-simulations>`_",
> ((Real,defaultDt,-1,,"used as the initial value of
> the timestep (especially useful in the first steps when no contact exist).
> If negative, it will be defined by :yref:`utils.PWaveTimeStep` *
> :yref:`GlobalStiffnessTimeStepper::timestepSafetyCoefficient`"))
> ((Real,maxDt,Mathr::MAX_REAL,,"if positive, used as
> max value of the timestep whatever the computed value"))
> ((Real,previousDt,1,,"last computed dt |yupdate|"))
> ((Real,timestepSafetyCoefficient,0.8,,"safety factor
> between the minimum eigen-period and the final assigned dt (less than 1))"))
> ((bool,densityScaling,false,,"|yupdate| don't modify
> this value if you don't plan to modify the scaling factor manually for some
> bodies. In most cases, it is enough to set
> :yref:`NewtonIntegrator::densityScaling` and let this one be adjusted
> automatically."))
> ((Real,targetDt,1,,"if
> :yref:`NewtonIntegrator::densityScaling` is active, this value will be used
> as the simulation timestep and the scaling will use this value of dt as the
> target value. The value of targetDt is arbitrary and should have no effect
> in the result in general. However if some bodies have imposed velocities,
> for instance, they will move more or less per each step depending on this
> value.")),
> computedOnce=false;)
> DECLARE_LOGGER;
> };
>
> REGISTER_SERIALIZABLE(GlobalStiffnessTimeStepper);
>
>
>
>
>
> _______________________________________________
> Mailing list: https://launchpad.net/~yade-dev
> Post to : yade-dev@xxxxxxxxxxxxxxxxxxx
> Unsubscribe : https://launchpad.net/~yade-dev
> More help : https://help.launchpad.net/ListHelp
>
diff --git a/pkg/dem/GlobalStiffnessTimeStepper.cpp b/pkg/dem/GlobalStiffnessTimeStepper.cpp
index 8a41ceb..0205655 100644
--- a/pkg/dem/GlobalStiffnessTimeStepper.cpp
+++ b/pkg/dem/GlobalStiffnessTimeStepper.cpp
@@ -14,6 +14,7 @@
#include<yade/core/Scene.hpp>
#include<yade/core/Clump.hpp>
#include<yade/pkg/dem/Shop.hpp>
+#include<yade/pkg/dem/ViscoelasticPM.hpp>
CREATE_LOGGER(GlobalStiffnessTimeStepper);
YADE_PLUGIN((GlobalStiffnessTimeStepper));
@@ -25,13 +26,16 @@ void GlobalStiffnessTimeStepper::findTimeStepFromBody(const shared_ptr<Body>& bo
State* sdec=body->state.get();
Vector3r& stiffness= stiffnesses[body->getId()];
Vector3r& Rstiffness=Rstiffnesses[body->getId()];
-
if(body->isClump()) {// if clump, we sum stifnesses of all members
const shared_ptr<Clump>& clump=YADE_PTR_CAST<Clump>(body->shape);
FOREACH(Clump::MemberMap::value_type& B, clump->members){
const shared_ptr<Body>& b = Body::byId(B.first,scene);
stiffness+=stiffnesses[b->getId()];
Rstiffness+=Rstiffnesses[b->getId()];
+ if (viscEl == true){
+ viscosities[body->getId()]+=viscosities[b->getId()];
+ Rviscosities[body->getId()]+=Rviscosities[b->getId()];
+ }
}
}
@@ -39,13 +43,15 @@ void GlobalStiffnessTimeStepper::findTimeStepFromBody(const shared_ptr<Body>& bo
if (densityScaling) sdec->densityScaling = min(1.0001*sdec->densityScaling, timestepSafetyCoefficient*pow(defaultDt/targetDt,2.0));
return; // not possible to compute!
}
-
+
+ //Determine the elastic minimum eigenperiod (and if required determine also the viscous one separately and take the minimum of the two)
+
+ //Elastic
Real dtx, dty, dtz;
Real dt = max( max (stiffness.x(), stiffness.y()), stiffness.z() );
if (dt!=0) {
dt = sdec->mass/dt; computedSomething = true;}//dt = squared eigenperiod of translational motion
else dt = Mathr::MAX_REAL;
-
if (Rstiffness.x()!=0) {
dtx = sdec->inertia.x()/Rstiffness.x(); computedSomething = true;}//dtx = squared eigenperiod of rotational motion around x
else dtx = Mathr::MAX_REAL;
@@ -56,15 +62,43 @@ void GlobalStiffnessTimeStepper::findTimeStepFromBody(const shared_ptr<Body>& bo
dtz = sdec->inertia.z()/Rstiffness.z(); computedSomething = true;}
else dtz = Mathr::MAX_REAL;
- Real Rdt = std::min( std::min (dtx, dty), dtz );//Rdt = smallest squared eigenperiod for rotational motions
+ Real Rdt = std::min( std::min (dtx, dty), dtz );//Rdt = smallest squared eigenperiod for elastic rotational motions
dt = 1.41044*timestepSafetyCoefficient*std::sqrt(std::min(dt,Rdt));//1.41044 = sqrt(2)
+
+ //Viscous
+ if (viscEl == true){
+ Vector3r& viscosity = viscosities[body->getId()];
+ Vector3r& Rviscosity = Rviscosities[body->getId()];
+ Real dtx_visc, dty_visc, dtz_visc;
+ Real dt_visc = max(max(viscosity.x(), viscosity.y()), viscosity.z() );
+ if (dt_visc!=0) {
+ dt_visc = sdec->mass/dt_visc; computedSomething = true;}//dt = eigenperiod of the viscous translational motion
+ else {dt_visc = Mathr::MAX_REAL;}
+
+ if (Rviscosity.x()!=0) {
+ dtx_visc = sdec->inertia.x()/Rviscosity.x(); computedSomething = true;}//dtx = eigenperiod of viscous rotational motion around x
+ else dtx_visc = Mathr::MAX_REAL;
+ if (Rviscosity.y()!=0) {
+ dty_visc = sdec->inertia.y()/Rviscosity.y(); computedSomething = true;}
+ else dty_visc = Mathr::MAX_REAL;
+ if (Rviscosity.z()!=0) {
+ dtz_visc = sdec->inertia.z()/Rviscosity.z(); computedSomething = true;}
+ else dtz_visc = Mathr::MAX_REAL;
+
+ Real Rdt_visc = std::min( std::min (dtx_visc, dty_visc), dtz_visc );//Rdt = smallest squared eigenperiod for viscous rotational motions
+ dt_visc = 2*timestepSafetyCoefficient*std::min(dt_visc,Rdt_visc);
+
+ //Take the minimum between the elastic and viscous minimum eigenperiod.
+ dt = std::min(dt,dt_visc);
+ }
+
//if there is a target dt, then we apply density scaling on the body, the inertia used in Newton will be mass*scaling, the weight is unmodified
if (densityScaling) {
sdec->densityScaling = min(sdec->densityScaling,timestepSafetyCoefficient*pow(dt /targetDt,2.0));
newDt=targetDt;
}
//else we update dt normaly
- else {newDt = std::min(dt,newDt);}
+ else {newDt = std::min(dt,newDt);}
}
bool GlobalStiffnessTimeStepper::isActivated()
@@ -97,7 +131,7 @@ void GlobalStiffnessTimeStepper::computeTimeStep(Scene* ncb)
computedOnce = true;}
else if (!computedOnce) scene->dt=defaultDt;
LOG_INFO("computed timestep " << newDt <<
- (scene->dt==newDt ? string(", appplied") :
+ (scene->dt==newDt ? string(", applied") :
string(", BUT timestep is ")+lexical_cast<string>(scene->dt))<<".");
}
@@ -107,15 +141,24 @@ void GlobalStiffnessTimeStepper::computeStiffnesses(Scene* rb){
if(size<rb->bodies->size()){
size=rb->bodies->size();
stiffnesses.resize(size); Rstiffnesses.resize(size);
+ if (viscEl == true){
+ viscosities.resize(size); Rviscosities.resize(size);
+ }
}
/* reset stored values */
memset(& stiffnesses[0],0,sizeof(Vector3r)*size);
memset(&Rstiffnesses[0],0,sizeof(Vector3r)*size);
+ if (viscEl == true){
+ memset(& viscosities[0],0,sizeof(Vector3r)*size);
+ memset(&Rviscosities[0],0,sizeof(Vector3r)*size);
+ }
+
FOREACH(const shared_ptr<Interaction>& contact, *rb->interactions){
if(!contact->isReal()) continue;
GenericSpheresContact* geom=YADE_CAST<GenericSpheresContact*>(contact->geom.get()); assert(geom);
NormShearPhys* phys=YADE_CAST<NormShearPhys*>(contact->phys.get()); assert(phys);
+
// all we need for getting stiffness
Vector3r& normal=geom->normal; Real& kn=phys->kn; Real& ks=phys->ks; Real& radius1=geom->refR1; Real& radius2=geom->refR2;
Real fn = (static_cast<NormShearPhys *> (contact->phys.get()))->normalForce.squaredNorm();
@@ -134,10 +177,35 @@ void GlobalStiffnessTimeStepper::computeStiffnesses(Scene* rb){
std::pow(normal.x(),2)+std::pow(normal.z(),2),
std::pow(normal.x(),2)+std::pow(normal.y(),2));
diag_Rstiffness *= ks;
+
+
//NOTE : contact laws with moments would be handled correctly by summing directly bending+twisting stiffness to diag_Rstiffness. The fact that there is no problem currently with e.g. cohesiveFrict law is probably because final computed dt is constrained by translational motion, not rotations.
stiffnesses [contact->getId1()]+=diag_stiffness;
Rstiffnesses[contact->getId1()]+=diag_Rstiffness*pow(radius1,2);
stiffnesses [contact->getId2()]+=diag_stiffness;
- Rstiffnesses[contact->getId2()]+=diag_Rstiffness*pow(radius2,2);
+ Rstiffnesses[contact->getId2()]+=diag_Rstiffness*pow(radius2,2);
+
+ //Same for the Viscous part, if required
+ if (viscEl == true){
+ ViscElPhys* viscPhys = YADE_CAST<ViscElPhys*>(contact->phys.get()); assert(viscPhys);
+ Real& cn=viscPhys->cn; Real& cs=viscPhys->cs;
+ //Diagonal terms of the translational viscous matrix
+ Vector3r diag_viscosity = Vector3r(std::pow(normal.x(),2),std::pow(normal.y(),2),std::pow(normal.z(),2));
+ diag_viscosity *= cn-cs;
+ diag_viscosity = diag_viscosity + Vector3r(1,1,1)*cs;
+ //diagonal terms of the rotational viscous matrix
+ Vector3r diag_Rviscosity =
+ Vector3r(std::pow(normal.y(),2)+std::pow(normal.z(),2),
+ std::pow(normal.x(),2)+std::pow(normal.z(),2),
+ std::pow(normal.x(),2)+std::pow(normal.y(),2));
+ diag_Rviscosity *= cs;
+
+ // Add the contact stiffness matrix to the two particles one
+ viscosities [contact->getId1()]+=diag_viscosity;
+ Rviscosities[contact->getId1()]+=diag_Rviscosity*pow(radius1,2);
+ viscosities [contact->getId2()]+=diag_viscosity;
+ Rviscosities[contact->getId2()]+=diag_Rviscosity*pow(radius2,2);
+ }
+
}
}
diff --git a/pkg/dem/GlobalStiffnessTimeStepper.hpp b/pkg/dem/GlobalStiffnessTimeStepper.hpp
index 1d60f88..bd657fa 100644
--- a/pkg/dem/GlobalStiffnessTimeStepper.hpp
+++ b/pkg/dem/GlobalStiffnessTimeStepper.hpp
@@ -24,6 +24,8 @@ class GlobalStiffnessTimeStepper : public TimeStepper
private :
vector<Vector3r> stiffnesses;
vector<Vector3r> Rstiffnesses;
+ vector<Vector3r> viscosities;
+ vector<Vector3r> Rviscosities;
void computeStiffnesses(Scene*);
Real newDt;
@@ -37,13 +39,14 @@ class GlobalStiffnessTimeStepper : public TimeStepper
virtual void computeTimeStep(Scene*);
virtual bool isActivated();
YADE_CLASS_BASE_DOC_ATTRS_CTOR(
- GlobalStiffnessTimeStepper,TimeStepper,"An engine assigning the time-step as a fraction of the minimum eigen-period in the problem. The derivation is detailed in the chapter on `DEM formulation <https://www.yade-dem.org/doc/formulation.html#dem-simulations>`_",
+ GlobalStiffnessTimeStepper,TimeStepper,"An engine assigning the time-step as a fraction of the minimum eigen-period in the problem. The derivation is detailed in the chapter on `DEM formulation <https://www.yade-dem.org/doc/formulation.html#dem-simulations>`_. The viscEl option enables to evaluate the timestep in a similar way for the visco-elastic contact law :yref:`Law2_ScGeom_ViscElPhys_Basic`, more detail in :yref:`GlobalStiffnessTimestepper::viscEl`. ",
((Real,defaultDt,-1,,"used as the initial value of the timestep (especially useful in the first steps when no contact exist). If negative, it will be defined by :yref:`utils.PWaveTimeStep` * :yref:`GlobalStiffnessTimeStepper::timestepSafetyCoefficient`"))
((Real,maxDt,Mathr::MAX_REAL,,"if positive, used as max value of the timestep whatever the computed value"))
((Real,previousDt,1,,"last computed dt |yupdate|"))
((Real,timestepSafetyCoefficient,0.8,,"safety factor between the minimum eigen-period and the final assigned dt (less than 1))"))
((bool,densityScaling,false,,"|yupdate| don't modify this value if you don't plan to modify the scaling factor manually for some bodies. In most cases, it is enough to set :yref:`NewtonIntegrator::densityScaling` and let this one be adjusted automatically."))
- ((Real,targetDt,1,,"if :yref:`NewtonIntegrator::densityScaling` is active, this value will be used as the simulation timestep and the scaling will use this value of dt as the target value. The value of targetDt is arbitrary and should have no effect in the result in general. However if some bodies have imposed velocities, for instance, they will move more or less per each step depending on this value.")),
+ ((Real,targetDt,1,,"if :yref:`NewtonIntegrator::densityScaling` is active, this value will be used as the simulation timestep and the scaling will use this value of dt as the target value. The value of targetDt is arbitrary and should have no effect in the result in general. However if some bodies have imposed velocities, for instance, they will move more or less per each step depending on this value."))
+ ((bool,viscEl,false,,"To use with :yref:`ViscElPhys`. if True, evaluate separetly the minimum eigen-period in the problem considering only the elastic contribution on one hand (spring only), and only the viscous contribution on the other hand (dashpot only). Take then the minimum of the two and use the safety coefficient :yref:`GlobalStiffnessTimestepper::timestepSafetyCoefficient` to take into account the possible coupling between the two contribution.")),
computedOnce=false;)
DECLARE_LOGGER;
};
References
