yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #10019
[Branch ~yade-pkg/yade/git-trunk] Rev 3707: Modify GlobalStiffnessTimeStepper for visco-elastic contact law.
------------------------------------------------------------
revno: 3707
author: Raphaël Maurin <raphael.maurin@xxxxxxxxx>
committer: Anton Gladky <gladky.anton@xxxxxxxxx>
timestamp: Thu 2013-09-26 09:51:11 +0200
message:
Modify GlobalStiffnessTimeStepper for visco-elastic contact law.
modified:
pkg/dem/GlobalStiffnessTimeStepper.cpp
pkg/dem/GlobalStiffnessTimeStepper.hpp
--
lp:yade
https://code.launchpad.net/~yade-pkg/yade/git-trunk
Your team Yade developers is subscribed to branch lp:yade.
To unsubscribe from this branch go to https://code.launchpad.net/~yade-pkg/yade/git-trunk/+edit-subscription
=== modified file 'pkg/dem/GlobalStiffnessTimeStepper.cpp'
--- pkg/dem/GlobalStiffnessTimeStepper.cpp 2013-05-24 16:48:19 +0000
+++ pkg/dem/GlobalStiffnessTimeStepper.cpp 2013-09-26 07:51:11 +0000
@@ -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 @@
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 @@
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 @@
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 @@
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 @@
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 @@
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);
+ }
+
}
}
=== modified file 'pkg/dem/GlobalStiffnessTimeStepper.hpp'
--- pkg/dem/GlobalStiffnessTimeStepper.hpp 2013-04-09 14:14:48 +0000
+++ pkg/dem/GlobalStiffnessTimeStepper.hpp 2013-09-26 07:51:11 +0000
@@ -24,6 +24,8 @@
private :
vector<Vector3r> stiffnesses;
vector<Vector3r> Rstiffnesses;
+ vector<Vector3r> viscosities;
+ vector<Vector3r> Rviscosities;
void computeStiffnesses(Scene*);
Real newDt;
@@ -37,13 +39,14 @@
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;
};
Follow ups
