yade-dev team mailing list archive

[noreply@xxxxxxxxxxxxx]

------------------------------------------------------------
revno: 2069
committer: Luc Scholtes <sch50p@fluent-ph>
branch nick: trunk
timestamp: Wed 2010-03-10 09:43:53 +1000
message:
  NEW TRY to add the updated version of cohesiveFrictionalContactLaw called cohesiveFrictionalPM (see revision 2068)
added:
  pkg/dem/Engine/GlobalEngine/CohesiveFrictionalPM.cpp
  pkg/dem/Engine/GlobalEngine/CohesiveFrictionalPM.hpp


--
lp:yade
https://code.launchpad.net/~yade-dev/yade/trunk

Your team Yade developers is subscribed to branch lp:yade.
To unsubscribe from this branch go to https://code.launchpad.net/~yade-dev/yade/trunk/+edit-subscription.

=== added file 'pkg/dem/Engine/GlobalEngine/CohesiveFrictionalPM.cpp'
--- pkg/dem/Engine/GlobalEngine/CohesiveFrictionalPM.cpp	1970-01-01 00:00:00 +0000
+++ pkg/dem/Engine/GlobalEngine/CohesiveFrictionalPM.cpp	2010-03-09 23:43:53 +0000
@@ -0,0 +1,189 @@
+/* LucScholtes2010  */
+
+#include"CohesiveFrictionalPM.hpp"
+#include<yade/core/Scene.hpp>
+#include<yade/pkg-dem/ScGeom.hpp>
+#include<yade/core/Omega.hpp>
+
+YADE_PLUGIN((Law2_ScGeom_CFpmPhys_CohesiveFrictionalPM)(Ip2_CFpmMat_CFpmMat_CFpmPhys)(CFpmPhys)(CFpmMat));
+
+/********************** Law2_ScGeom_CFpmPhys_CohesiveFrictionalPM ****************************/
+CREATE_LOGGER(Law2_ScGeom_CFpmPhys_CohesiveFrictionalPM);
+
+void Law2_ScGeom_CFpmPhys_CohesiveFrictionalPM::go(shared_ptr<InteractionGeometry>& ig, shared_ptr<InteractionPhysics>& ip, Interaction* contact, Scene* rootBody){
+	
+	ScGeom* geom = static_cast<ScGeom*>(ig.get()); 
+	CFpmPhys* phys = static_cast<CFpmPhys*>(ip.get());
+	int id1 = contact->getId1();
+        int id2 = contact->getId2();
+	shared_ptr<BodyContainer>& bodies = rootBody->bodies;
+	Body* b1 = ( *bodies ) [id1].get();
+	Body* b2 = ( *bodies ) [id2].get();
+	Real dt = Omega::instance().getTimeStep();
+	
+	Real displN = geom->penetrationDepth; // NOTE: the sign for penetrationdepth is different from ScGeom and Dem3DofGeom: geom->penetrationDepth>0 when spheres interpenetrate
+	Real D = 0; // initial interparticular distance
+	Real Dtensile=phys->FnMax/phys->kn;
+	Real Dsoftening = phys->strengthSoftening*Dtensile; 
+	
+	/*to set the equilibrium distance between all cohesive elements when they first meet*/
+	if ( contact->isFresh(rootBody) && phys->isCohesive) { phys->initD = displN; phys->normalForce = Vector3r::ZERO; phys->shearForce = Vector3r::ZERO;}
+	
+	D = displN - phys->initD; // displacement is computed depending on the equilibrium distance
+	
+	/* Determination of interaction */
+	if (D < 0){ //spheres do not touch 
+	  if (!phys->isCohesive){ rootBody->interactions->requestErase(contact->getId1(),contact->getId2()); return; } // destroy the interaction before calculation
+	  if ((phys->isCohesive) && (abs(D) > (Dtensile + Dsoftening))) { // spheres are bonded and the interacting distance is greater than the one allowed ny the defined cohesion
+	    phys->isCohesive=false; 
+	    rootBody->interactions->requestErase(contact->getId1(),contact->getId2()); return;
+	  }
+	}	  
+	
+	/*NormalForce*/
+	Real Fn=0, Dsoft=0;
+		
+	if ((D < 0) && (abs(D) > Dtensile)) { //to take into account strength softening
+	  Dsoft = D+Dtensile; // Dsoft<0 for a negative value of Fn (attractive force)
+	  Fn = phys->FnMax+(phys->kn/phys->strengthSoftening)*Dsoft; // computes FnMax - FnSoftening
+	}
+	else {
+	  Fn = phys->kn*D;
+	}
+		
+	phys->normalForce = Fn*geom->normal;  // NOTE normal is position2-position1 - It is directed from particle1 to particle2
+	        
+	/*ShearForce*/
+	Vector3r& shearForce = phys->shearForce; 
+		
+	// using scGeom function updateShear(Vector3r& shearForce, Real ks, const Vector3r& prevNormal, const State* rbp1, const State* rbp2, Real dt, bool avoidGranularRatcheting)	
+	State* st1 = Body::byId(id1,rootBody)->state.get();
+	State* st2 = Body::byId(id2,rootBody)->state.get();
+	geom->updateShearForce(shearForce, phys->ks, phys->prevNormal, st1, st2, dt, preventGranularRatcheting);
+	
+	// needed for the next timestep
+	phys->prevNormal = geom->normal;
+		
+	/* Morh-Coulomb criterion */
+	Vector3r Fs = shearForce; // NOTE the following could be done directly with shearForce?
+	Real maxFs = abs(Fn)*phys->tanFrictionAngle + phys->FsMax;
+		
+	if (Fs.SquaredLength() > maxFs*maxFs){ 
+	  phys->isCohesive=false; //if not already false to delete cohesive interactions
+	  Fs*=maxFs/Fs.Length(); // to fix the shear force to its yielding value with the right sign
+	}
+	phys->shearForce = Fs; 
+	
+	/* Apply forces */
+	Vector3r f = phys->normalForce + phys->shearForce;
+	
+	// NOTE applyForceAtContactPoint computes torque induced by normal and shear force too
+ 	applyForceAtContactPoint(f, geom->contactPoint, contact->getId2(), b2->state->pos, contact->getId1(), b1->state->pos, rootBody); // NOTE id1 and id2 must suit the orientation of the force
+
+	/* Moment Rotation Law */
+	// NOTE this part could probably be computed in ScGeom to avoid copy/paste multiplication !!!
+	Quaternionr delta( b1->state->ori * phys->initialOrientation1.Conjugate() *phys->initialOrientation2 * b2->state->ori.Conjugate()); //relative orientation
+	Vector3r axisM;	// axis of rotation - this is the Moment direction UNIT vector.
+	Real angleM;	// angle represents the power of resistant ELASTIC moment
+	delta.ToAxisAngle(axisM,angleM);
+	if(angleM > Mathr::PI) angleM -= Mathr::TWO_PI; // angle is between 0 and 2*pi, but should be between -pi and pi 
+	  
+	phys->cumulativeRotation = angleM;
+	  
+	//Find angle*axis. That's all.  But first find angle about contact normal. Result is scalar. Axis is contact normal.
+	Real angle_twist(angleM * axisM.Dot(geom->normal) ); //rotation about normal
+	Vector3r axis_twist(angle_twist * geom->normal);
+	Vector3r moment_twist(axis_twist * phys->kr);
+	  
+	Vector3r axis_bending(angleM*axisM - axis_twist); //total rotation minus rotation about normal
+	Vector3r moment_bending(axis_bending * phys->kr);
+	Vector3r moment = moment_twist + moment_bending;
+
+	Real MomentMax = phys->maxBend*std::fabs(phys->normalForce.Length());
+	Real scalarMoment = moment.Length();
+
+	/*Plastic moment */
+	if(scalarMoment > MomentMax) 
+	{
+	  Real ratio = 0;
+	  ratio *= MomentMax/scalarMoment; // to fix the moment to its yielding value
+	  moment *= ratio;		
+	  moment_twist *=  ratio;
+	  moment_bending *= ratio;
+	}
+
+	phys->moment_twist = moment_twist;
+	phys->moment_bending = moment_bending;
+	  
+	rootBody->forces.addTorque(id1,-moment);
+	rootBody->forces.addTorque(id2, moment);
+		
+}
+
+CREATE_LOGGER(Ip2_CFpmMat_CFpmMat_CFpmPhys);
+
+void Ip2_CFpmMat_CFpmMat_CFpmPhys::go(const shared_ptr<Material>& b1, const shared_ptr<Material>& b2, const shared_ptr<Interaction>& interaction){
+	
+	/* avoid any updates if the interaction already exists */
+	if(interaction->interactionPhysics) return; 
+
+	ScGeom* geom=dynamic_cast<ScGeom*>(interaction->interactionGeometry.get());
+	assert(geom);
+
+	const shared_ptr<CFpmMat>& yade1 = YADE_PTR_CAST<CFpmMat>(b1);
+	const shared_ptr<CFpmMat>& yade2 = YADE_PTR_CAST<CFpmMat>(b2);
+			
+	shared_ptr<CFpmPhys> contactPhysics(new CFpmPhys()); 
+	
+	/* From interaction physics */
+	Real E1 	= yade1->young;
+	Real E2 	= yade2->young;
+	Real V1 	= yade1->poisson;
+	Real V2 	= yade2->poisson;
+	Real f1 	= yade1->frictionAngle;
+	Real f2 	= yade2->frictionAngle;
+
+	/* From interaction geometry */
+	Real R1= geom->radius1;
+	Real R2= geom->radius2;
+	Real rMean = 0.5*(R1+R2); 
+	Real crossSection = Mathr::PI*pow(min(R1,R2),2); 
+	
+	/* calculate stiffness */
+	Real kNormal=0, kShear=0, kRotate=0;
+	
+	if(useAlphaBeta == true){
+	  kNormal	= 2.*E1*R1*E2*R2/(E1*R1+E2*R2); // harmonic average of two stiffnesses
+	  kShear 	= Alpha*kNormal;
+	  kRotate	= Beta*kShear*rMean*rMean; 
+	}
+	else {
+	  kNormal	= 2.*E1*R1*E2*R2/(E1*R1+E2*R2); // harmonic average of two stiffnesses
+	  kShear	= 2.*E1*R1*V1*E2*R2*V2/(E1*R1*V1+E2*R2*V2); // harmonic average of two stiffnesses with ks=V*kn for each sphere
+	  kRotate	= 0.;
+	}
+
+	/* Pass values calculated from above to CFpmPhys */
+	contactPhysics->kn = kNormal;
+	contactPhysics->ks = kShear;
+	contactPhysics->kr = kRotate;
+	contactPhysics->frictionAngle		= std::min(f1,f2); 
+	contactPhysics->tanFrictionAngle	= std::tan(contactPhysics->frictionAngle); 
+	contactPhysics->maxBend 		= eta*rMean;
+	contactPhysics->prevNormal 		= geom->normal;
+	contactPhysics->initialOrientation1	= Body::byId(interaction->getId1())->state->ori;
+	contactPhysics->initialOrientation2	= Body::byId(interaction->getId2())->state->ori;
+	
+	///to set if the contact is cohesive or not
+	if ( (scene->currentIteration < cohesiveTresholdIteration) && ((tensileStrength>0) || (cohesion>0)) && (yade1->type == yade2->type)){ contactPhysics->isCohesive=true; }
+	
+	if ( contactPhysics->isCohesive ) {
+	  contactPhysics->FnMax = tensileStrength*crossSection;
+	  contactPhysics->strengthSoftening = strengthSoftening;
+	  contactPhysics->FsMax = cohesion*crossSection;
+	}
+	
+	interaction->interactionPhysics = contactPhysics;
+}
+
+CFpmPhys::~CFpmPhys(){}
\ No newline at end of file

=== added file 'pkg/dem/Engine/GlobalEngine/CohesiveFrictionalPM.hpp'
--- pkg/dem/Engine/GlobalEngine/CohesiveFrictionalPM.hpp	1970-01-01 00:00:00 +0000
+++ pkg/dem/Engine/GlobalEngine/CohesiveFrictionalPM.hpp	2010-03-09 23:43:53 +0000
@@ -0,0 +1,101 @@
+/* lucScholtes2010 */
+
+#pragma once
+
+#include<yade/pkg-common/ElastMat.hpp>
+#include<yade/pkg-common/InteractionPhysicsFunctor.hpp>
+#include<yade/pkg-common/NormShearPhys.hpp>
+#include<yade/pkg-common/LawFunctor.hpp>
+#include<yade/pkg-dem/ScGeom.hpp>
+
+/*
+=== OVERVIEW OF CohesiveFrictionalPM ===
+
+CohesiveFrictional Particle Model (CohesiveFrictionalPM, CFpm) is a set of classes for modelling mechanical behavior of cohesive frictional material. It is a quite general contact model starting from the basic frictional model with the possible addition of a moment between particles, a tensile strength and a cohesion.
+
+Features of the interaction law:
+
+1.  If useAlphaBeta=False, Kn, Ks are calculated from the harmonic average of the "macro" material parameters young (E) and poisson (V) as defined in SimpleElasticRelationships.cpp and Kr=0.
+
+2.  If useAlphaBeta=True, users have to input Alpha=Ks/Kn, Beta=Kr/(Ks*meanRadius^2), eta=MtPlastic/(meanRadius*Fn) as defined in Plassiard et al. (Granular Matter, 2009) A spherical discrete element model.
+
+3.  Users can input a tensile strength ( FnMax = tensileStrength*) and a cohesion ( FsMax = cohesion*(Ks*meanRadius))
+
+  Remark: - This contact law works well in the case of sphere/sphere and sphere/box interaction as it uses ScGeom to compute the interaction geometry (suitable for the triaxialTest) 
+	  - It has not been tested for sphere/facet or sphere/wall interactions and could be updated to be used by DemXDofGeom
+*/
+
+/** This class holds information associated with each body */
+class CFpmMat: public FrictMat {
+	YADE_CLASS_BASE_DOC_ATTRS_CTOR(CFpmMat,FrictMat,"cohesive frictional material, for use with other CFpm classes",
+	  ((int,type,0,"Type of the particle. If particles of two different types interact, it will be with friction only (no cohesion).[-]")),
+	  createIndex();
+	);
+	REGISTER_CLASS_INDEX(CFpmMat,FrictMat);
+};
+REGISTER_SERIALIZABLE(CFpmMat);
+
+/** This class holds information associated with each interaction */
+class CFpmPhys: public NormShearPhys {
+	public:
+		virtual ~CFpmPhys();
+	
+		YADE_CLASS_BASE_DOC_ATTRS_CTOR_PY(CFpmPhys,NormShearPhys,"Representation of a single interaction of the CFpm type, storage for relevant parameters",
+		  ((Real,initD,0,"equilibrium distance for particles. Computed as the initial interparticular distance when bonded particle interact. initD=0 for non cohesive interactions."))
+		  ((bool,isCohesive,false,"If false, particles interact in a frictional way. If true, particles are bonded regarding the given cohesion and tensileStrength."))
+		  ((Real, frictionAngle,0,"defines Coulomb friction. [deg]"))
+		  ((Real,tanFrictionAngle,0,"Tangent of frictionAngle. [-]"))
+		  ((Real,FnMax,0,"Defines the maximum admissible normal force in traction FnMax=tensileStrength*crossSection, with crossSection=pi*Rmin^2. [Pa]"))
+		  ((Real,FsMax,0,"Defines the maximum admissible tangential force in shear FsMax=cohesion*FnMax, with crossSection=pi*Rmin^2. [Pa]"))
+		  ((Real,strengthSoftening,0,"Defines the softening when Dtensile is reached to avoid explosion. Typically, when D > Dtensile, Fn=FnMax - (kn/strengthSoftening)*(Dtensile-D). [-]"))
+		  ((Real,cumulativeRotation,0,"Cumulated rotation... [-]"))
+		  ((Real,kr,0,"Defines the stiffness to compute the resistive moment in rotation. [-]"))
+		  ((Real,maxBend,0,"Defines the maximum admissible resistive moment in rotation Mtmax=maxBend*Fn, maxBend=eta*meanRadius. [m]"))
+		  ((Vector3r,prevNormal,Vector3r::ZERO,"Normal to the contact at previous time step."))
+		  ((Vector3r,moment_twist,Vector3r::ZERO," [N.m]"))
+		  ((Vector3r,moment_bending,Vector3r::ZERO," [N.m]"))
+		  ((Quaternionr,initialOrientation1,Quaternionr(1.0,0.0,0.0,0.0),"Use for moment computation."))
+		  ((Quaternionr,initialOrientation2,Quaternionr(1.0,0.0,0.0,0.0),"Use for moment computation."))
+		  ,
+		  createIndex();
+		  ,
+		);
+	DECLARE_LOGGER;
+	REGISTER_CLASS_INDEX(CFpmPhys,NormShearPhys);
+};
+REGISTER_SERIALIZABLE(CFpmPhys);
+
+/** 2d functor creating InteractionPhysics (Ip2) taking CFpmMat and CFpmMat of 2 bodies, returning type CFpmPhys */
+class Ip2_CFpmMat_CFpmMat_CFpmPhys: public InteractionPhysicsFunctor{
+	public:
+		virtual void go(const shared_ptr<Material>& pp1, const shared_ptr<Material>& pp2, const shared_ptr<Interaction>& interaction);
+		
+		FUNCTOR2D(CFpmMat,CFpmMat);
+		DECLARE_LOGGER;
+		
+		YADE_CLASS_BASE_DOC_ATTRS(Ip2_CFpmMat_CFpmMat_CFpmPhys,InteractionPhysicsFunctor,"Converts 2 CFpmmat instances to CFpmPhys with corresponding parameters.",
+		  ((int,cohesiveTresholdIteration,1,"Should new contacts be cohesive? They will before this iter, they won't afterward."))
+		  ((bool,useAlphaBeta,false,"If true, stiffnesses are computed based on Alpha and Beta."))
+		  ((Real,Alpha,0,"Defines the ratio ks/kn."))
+		  ((Real,Beta,0,"Defines the ratio kr/(ks*meanRadius^2) to compute the resistive moment in rotation. [-]"))
+		  ((Real,eta,0,"Defines the maximum admissible resistive moment in rotation MtMax=eta*meanRadius*Fn. [-]"))
+		  ((Real,tensileStrength,0,"Defines the maximum admissible normal force in traction FnMax=tensileStrength*crossSection. [Pa]"))
+		  ((Real,cohesion,0,"Defines the maximum admissible tangential force in shear FsMax=cohesion*crossSection. [Pa]"))
+		  ((Real,strengthSoftening,0,"Defines the softening when Dtensile is reached to avoid explosion of the contact. Typically, when D > Dtensile, Fn=FnMax - (kn/strengthSoftening)*(Dtensile-D). [-]"))
+		);
+};
+REGISTER_SERIALIZABLE(Ip2_CFpmMat_CFpmMat_CFpmPhys);
+
+
+/** 2d functor creating the interaction law (Law2) based on SphereContactGeometry (ScGeom) and CFpmPhys of 2 bodies, returning type CohesiveFrictionalPM */
+class Law2_ScGeom_CFpmPhys_CohesiveFrictionalPM: public LawFunctor{
+	public:
+		virtual void go(shared_ptr<InteractionGeometry>& _geom, shared_ptr<InteractionPhysics>& _phys, Interaction* I, Scene* rootBody);
+		FUNCTOR2D(ScGeom,CFpmPhys);
+
+	YADE_CLASS_BASE_DOC_ATTRS(Law2_ScGeom_CFpmPhys_CohesiveFrictionalPM,LawFunctor,"Constitutive law for the CFpm model.",
+		  ((bool,preventGranularRatcheting,false,"If true rotations are computed such as granular ratcheting is prevented. See F. ALONSO-MARROQUIN, R. GARCIA-ROJO, H.J. HERRMANN, Micro-mechanical investigation of granular ratcheting, in Cyclic Behaviour of Soils and Liquefaction Phenomena, ed. T. Triantafyllidis (Balklema, London, 2004), p. 3-10 - and a lot more papers from the same authors)."))
+	);
+	DECLARE_LOGGER;	
+};
+REGISTER_SERIALIZABLE(Law2_ScGeom_CFpmPhys_CohesiveFrictionalPM);