Hi Bruno, hi Anton,
I think density scaling is not what I am looking for since I am interested in
real dynamic simulations. I need different damping parameters in
NewtonIntegrator for the block and the mesh. I have to consider free flight
under gravity so damping=0 (and I am interested in the real flight time). For
the mesh I have to consider additional energy absorption which is not
considered in my model via damping.
It is similar to b_damp() in PFC which is used to remove damping for certain
particle. I implemented it the same way now. The state has a member which is
true by default which means that damping is used. It can be set to false and
damping=0 will be used for this particle. So basically damping in
NewtonIntegrator can be activated and deactivated for individual particles.
@Anton: It has to be checked for all particles.
Here is the diff of my latest implementation, tell me what you think:
=== modified file core/State.hpp
--- core/State.hpp 2011-02-14 08:05:09 +0000
+++ core/State.hpp 2012-01-22 23:56:31 +0000
@@ -59,7 +59,8 @@
((Vector3r,inertia,Vector3r::Zero(),,"Inertia of associated
body, in
local coordinate system."))
((Vector3r,refPos,Vector3r::Zero(),,"Reference position"))
((Quaternionr,refOri,Quaternionr::Identity(),,"Reference
orientation"))
- ((unsigned,blockedDOFs,,,"[Will be overridden]")),
+ ((unsigned,blockedDOFs,,,"[Will be overridden]"))
+ ((bool,isDamped,true,,"Damping in :yref:`Newtonintegrator` can
be
deactivated for individual particles by setting this variable to FALSE. E.g.
damping is inappropriate for particles in free flight under gravity but it
might still be applicable to other particles in the same simulation.")),
/* additional initializers */
((pos,se3.position))
((ori,se3.orientation)),
=== modified file pkg/dem/NewtonIntegrator.cpp
--- pkg/dem/NewtonIntegrator.cpp 2011-11-30 17:39:33 +0000
+++ pkg/dem/NewtonIntegrator.cpp 2012-01-22 23:56:31 +0000
@@ -11,17 +11,18 @@
#include<yade/pkg/dem/Clump.hpp>
#include<yade/pkg/common/VelocityBins.hpp>
#include<yade/lib/base/Math.hpp>
YADE_PLUGIN((NewtonIntegrator));
CREATE_LOGGER(NewtonIntegrator);
// 1st order numerical damping
-void NewtonIntegrator::cundallDamp1st(Vector3r& force, const Vector3r& vel){
- for(int i=0; i<3; i++) force[i]*=1-damping*Mathr::Sign(force[i]*vel[i]);
+void NewtonIntegrator::cundallDamp1st(Vector3r& force, const Vector3r& vel,
const Real& dampcoeff){
+ for(int i=0; i<3; i++)
force[i]*=1-dampcoeff*Mathr::Sign(force[i]*vel[i]);
}
// 2nd order numerical damping
-void NewtonIntegrator::cundallDamp2nd(const Real& dt, const Vector3r& force,
const Vector3r& vel, Vector3r& accel){
- for(int i=0; i<3; i++) accel[i]*= 1 - damping*Mathr::Sign (
force[i]*(vel[i] + 0.5*dt*accel[i]) );
+void NewtonIntegrator::cundallDamp2nd(const Real& dt, const Vector3r& force,
const Vector3r& vel, Vector3r& accel, const Real& dampcoeff){
+ for(int i=0; i<3; i++) accel[i]*= 1 - dampcoeff*Mathr::Sign (
force[i]*(vel[i] + 0.5*dt*accel[i]) );
}
Vector3r NewtonIntegrator::computeAccel(const Vector3r& force, const Real&
mass, int blockedDOFs){
@@ -39,11 +40,13 @@
void NewtonIntegrator::updateEnergy(const shared_ptr<Body>& b, const State*
state, const Vector3r& fluctVel, const Vector3r& f, const Vector3r& m){
assert(b->isStandalone() || b->isClump());
- // always positive dissipation, by-component: |F_i|*|v_i|*damping*dt (|
T_i|*|ω_i|*damping*dt for rotations)
- if(damping!=0.){
- scene->energy-
>add(fluctVel.cwise().abs().dot(f.cwise().abs())*damping*scene-
>dt,"nonviscDamp",nonviscDampIx,/*non-incremental*/false);
+ // check if damping for this body is activated
+ Real dampcoeff=(state->isDamped ? damping : 0);
+ // always positive dissipation, by-component: |F_i|*|v_i|*dampcoeff*dt
(|
T_i|*|ω_i|*dampcoeff*dt for rotations)
+ if(dampcoeff!=0.){
+ scene->energy-
>add(fluctVel.cwise().abs().dot(f.cwise().abs())*dampcoeff*scene-
>dt,"nonviscDamp",nonviscDampIx,/*non-incremental*/false);
// when the aspherical integrator is used, torque is damped
instead of
ang acceleration; this code is only approximate
- scene->energy->add(state-
>angVel.cwise().abs().dot(m.cwise().abs())*damping*scene-
>dt,"nonviscDamp",nonviscDampIx,false);
+ scene->energy->add(state-
>angVel.cwise().abs().dot(m.cwise().abs())*dampcoeff*scene-
>dt,"nonviscDamp",nonviscDampIx,false);
}
// kinetic energy
Real Etrans=.5*state->mass*fluctVel.squaredNorm();
@@ -106,9 +109,13 @@
YADE_PARALLEL_FOREACH_BODY_BEGIN(const shared_ptr<Body>& b,
scene->bodies)
{
// clump members are handled inside clumps
if(unlikely(b->isClumpMember())) continue;
-
+
State* state=b->state.get(); const Body::id_t&
id=b->getId();
Vector3r f=gravity*state->mass, m=Vector3r::Zero();
+
+ // check if damping for this body is activated
+ Real dampcoeff=(state->isDamped ? damping : 0);
+
// clumps forces
if(b->isClump()) {
b->shape-
>cast<Clump>().addForceTorqueFromMembers(state,scene,f,m);
@@ -146,7 +153,7 @@
if (state->blockedDOFs!=State::DOF_ALL) {
// linear acceleration
Vector3r
linAccel=computeAccel(f,state->mass,state->blockedDOFs);
- cundallDamp2nd(dt,f,fluctVel,linAccel);
+
cundallDamp2nd(dt,f,fluctVel,linAccel,dampcoeff);
//This is the convective term, appearing in the
time derivation
of Cundall/Thornton expression (dx/dt=velGrad*pos ->
d²x/dt²=dvelGrad/dt*pos+velGrad*vel), negligible in many cases but not for
high speed large deformations (gaz or turbulent flow).
linAccel+=prevVelGrad*state->vel;
//finally update velocity
@@ -154,11 +161,11 @@
// angular acceleration
if(!useAspherical){ // uses angular velocity
Vector3r
angAccel=computeAngAccel(m,state->inertia,state-
>blockedDOFs);
-
cundallDamp2nd(dt,m,state->angVel,angAccel);
+
cundallDamp2nd(dt,m,state->angVel,angAccel,dampcoeff);
state->angVel+=dt*angAccel;
} else { // uses torque
for(int i=0; i<3; i++)
if(state->blockedDOFs &
State::axisDOF(i,true)) m[i]=0; // block DOFs here
- cundallDamp1st(m,state->angVel);
+
cundallDamp1st(m,state->angVel,dampcoeff);
}
}
=== modified file pkg/dem/NewtonIntegrator.hpp
--- pkg/dem/NewtonIntegrator.hpp 2011-09-20 10:58:18 +0000
+++ pkg/dem/NewtonIntegrator.hpp 2012-01-12 05:26:05 +0000
@@ -22,8 +22,8 @@
class VelocityBins;
class NewtonIntegrator : public GlobalEngine{
- inline void cundallDamp1st(Vector3r& force, const Vector3r& vel);
- inline void cundallDamp2nd(const Real& dt, const Vector3r& force, const
Vector3r& vel, Vector3r& accel);
+ inline void cundallDamp1st(Vector3r& force, const Vector3r& vel, const
Real& dampcoeff);
+ inline void cundallDamp2nd(const Real& dt, const Vector3r& force, const
Vector3r& vel, Vector3r& accel, const Real& dampcoeff);
bool haveBins;
inline void leapfrogTranslate(State*, const Body::id_t& id, const Real&
dt); // leap-frog translate
inline void leapfrogSphericalRotate(State*, const Body::id_t& id, const
Real& dt); // leap-frog rotate of spherical body
On Fri, 20 Jan 2012 02:06:51 AM Bruno Chareyre wrote:
> Klaus,
> If what you want is to increase the timestep (because your steel-net
> induces very small dt, right?), this change will not work. You need to
> modify density in the good place in Newton, not damping. It can be done
> by adding dscale parameter in the material class (hence no need to check
> if it exists), which would be 1 by default, hence no need to check if it
> is NaN.
> isnan is not in C++ standard, by the way. I know it is used in many
> places but we should avoid it if possible. I'm also not sure of the cost
> of this function, it may be more expensive than comparing two doubles.
>
> Bruno
>
>
> _______________________________________________
> Mailing list: https://launchpad.net/~yade-dev
> Post to : yade-dev@lists.launchpad.net
> Unsubscribe : https://launchpad.net/~yade-dev
> More help : https://help.launchpad.net/ListHelp
_______________________________________________
Mailing list: https://launchpad.net/~yade-dev
Post to : yade-dev@lists.launchpad.net
Unsubscribe : https://launchpad.net/~yade-dev
More help : https://help.launchpad.net/ListHelp
