yade-dev team mailing list archive

[Klaus Thoeni <klaus.thoeni@xxxxxxxxx>]

Hi Chiara,

thanks for your comment. I am not after viscous damping either. But I think 
the current implementation does exactly what I want, it's the same as the PFC 
option b_damp(). So I can deactivate the damping coefficient for individual 
bodies.

I can commit the code if anyone is interested , after the diff has been 
approved ;-)

Klaus

On Mon, 23 Jan 2012 08:44:02 PM Chiara wrote:
> 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
> 
> _______________________________________________
> 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