← Back to team overview

yade-dev team mailing list archive

  1. yade-dev team
  2. Mailing list archive
  3. Message #08239

Re: : Add damping coefficient to material

  Thread Previous • Date Previous • • Thread Next 
On 23/01/12 00:12, Klaus Thoeni wrote:

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.
If you are not after any density scaling, why damping in NewtonIntegrator should be acceptable in a realistic dynamic situation? Would not be easier/more meaningful for you to implement damping inside your contact model? Sorry, just curious. 
Chiara


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@xxxxxxxxxxxxxxxxxxx
Unsubscribe : https://launchpad.net/~yade-dev
More help   : https://help.launchpad.net/ListHelp

_______________________________________________
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

Follow ups

References

  Thread Previous • Date Previous • • Thread Next