yade-users team mailing list archive

[Rémi Chassagne <question695457@xxxxxxxxxxxxxxxxxxxxx>]

New question #695457 on Yade:
https://answers.launchpad.net/yade/+question/695457

Dear all,

I am trying to compute the dissipated power of a granular system in a YADE simulation. In my simulations, I use a linear spring dashpot system (Law2_ScGeom_ViscElPhys_Basic()), that is to say:
Fn = kn delta_n + cn d(delta_n)/dt    (normal direction of contact)
Ft = -min(ks delta_t, muFn)                   (tangential direction of contact)

In the normal direction, only the damper dissipates energy so the normal dissipated power at contact is
Pn = - cn d(delta_n)/dt * d(delta_n)/dt
Considering a collision of two particles with only a normal component (I ran a yade simulation), I have been able to verify that the integrate of Pn during the contact time indeed corresponds to the theoretical energy dissipated for the corresponding restitution coefficient.

My question concerns more the tangential dissipation. Following the same reasoning, I would say that:
if Ft = -muFn:
    Pt = -muFn*shearVel
where shearVel is the tangential component of the relative velocity between the two particles. However, in some of the contact laws (those for which the energy dissipation is already coded, for example in ElasticContactLaw.cpp), the tangential dissipation is computed differently. It is not the shearVel velocity which is considered but a velocity based on the excess of elastic force compared to the sliding force:
if Ft = -muFn:
    Pt = - muFn * ((Ft[t-1]-ks delta_t) - Ft)/(ks*dt)
where Ft[t-1] is the shear force at previous time step and delta_t corresponds to a the shear increment (I assume delta_t = shearVel*dt). You can find the corresponding portion of code of ElasticContactLaw.cpp just below. 

I don't really understand why the sliding velocity has to be computed on the difference of forces and what this physically means. If someone understands that or have a reference to share which gives an explanation I would be grateful.

Thanks in advance,
Rémi



Vector3r&       shearForce = geom->rotate(phys->shearForce);
const Vector3r& shearDisp  = geom->shearIncrement();
shearForce -= phys->ks * shearDisp;
Real maxFs = phys->normalForce.squaredNorm() * math::pow(phys->tangensOfFrictionAngle, 2);

if (!scene->trackEnergy && !traceEnergy) { //Update force but don't compute energy terms (see below))
	// PFC3d SlipModel, is using friction angle. CoulombCriterion
	if (shearForce.squaredNorm() > maxFs) {
		Real ratio = sqrt(maxFs) / shearForce.norm();
		shearForce *= ratio;
	}
} else {
	//almost the same with additional Vector3r instatinated for energy tracing,
	//duplicated block to make sure there is no cost for the instanciation of the vector when traceEnergy==false
	if (shearForce.squaredNorm() > maxFs) {
		Real     ratio      = sqrt(maxFs) / shearForce.norm();
		Vector3r trialForce = shearForce; //store prev force for definition of plastic slip
		//define the plastic work input and increment the total plastic energy dissipated
		shearForce *= ratio;
		Real dissip = ((1 / phys->ks) * (trialForce - shearForce)) /*plastic disp*/.dot(shearForce) /*active force*/;
        }
}

-- 
You received this question notification because your team yade-users is
an answer contact for Yade.