yade-dev team mailing list archive

[luc scholtes <lscholtes63@xxxxxxxxx>]

Hi all,

Before any commit, I would like to have your opinion on what I did for the
moment concerning my project (jointed rock mass simulation). This is really
basic stuff for the moment due to my very poor knowledge in programming, but
I hope you guys could give me some advices. Please have a look if you can.

I did not attach all the files needed to run a simulation due to a previous
sending that have never arrived, but do not hesitate to ask for them.  I
attached an animated gif as a very simple illustration too.

Here is the idea:

+ For the moment, the JointedRockPM is based on the RockPM. The formulation
is exactly the same except that discontinuities (joints) can be taken into
account in the computation. The near future will be a specific law for
spheres on joint to allow a smoother behavior (a modification of the
interacting geometry would be the first step -> probably some questions for
you guys), and probably a more complete law for the cohesive behaviour with
shear and moment cohesion (-> probably some more questions here too).

+ to take into account joints, I developed a procedure (maybe not very
efficient) based on python scripting (IdentificationSpheresOnJoint.
py and ImportJointedRockMass.py that I can send you if you want). This can
appear a bit complicated for the given illustration, but combining this with
some imports from a CAD would give some interesting stuffs (import of
complete Discrete Fracture Network for example...)

  - first, given a rock (or whatever you want ndlr) sample geometry, and a
"joints set" that cut through this sample (defined by facets), spheres
belonging to each side of the joints surface are identified with the joint
normal and the side to which they belong as attributes. An ascii files is
then saved through a modified function of the sphereToFile() function
defined in utils.py (oncemore, I can give you the concerned file if you
want).

  - second, the jointed rock mass is loaded from this ascii files through
the asciiJointedRockPM() and the simulation can start. Here is still some
weird coding for the import (ImportjointedRockMass.py), but it works. I
choose this saving/loading procedure to have a reset of the simulation
without the presence of the joints (which would give some unappropriate
behavior). Maybe this can be done on the same simulation, but I don't
know...


Oncemore, this is just a base. Every comments are welcome (...)


Regards

/**
=== HIGH LEVEL OVERVIEW OF JointedRockPM ===

Jointed Rock Particle Model (JointedRockPM) is a set of classes for modelling
mechanical behavior of jointed rock masses (directly derived from RockPM).
*/

#pragma once

#include<yade/pkg-common/InteractionPhysicsFunctor.hpp>
#include<yade/pkg-dem/ScGeom.hpp>
#include<yade/pkg-common/PeriodicEngines.hpp>
#include<yade/pkg-common/NormalShearInteractions.hpp>
#include<yade/pkg-common/LawFunctor.hpp>
#include<yade/pkg-common/ElasticMat.hpp>


class Law2_Dem3DofGeom_JointedRockPMPhys_JRpm: public LawFunctor{
	public:
		virtual void go(shared_ptr<InteractionGeometry>& _geom, shared_ptr<InteractionPhysics>& _phys, Interaction* I, Scene* rootBody);
		FUNCTOR2D(Dem3DofGeom,JRpmPhys);
		REGISTER_CLASS_AND_BASE(Law2_Dem3DofGeom_JointedRockPMPhys_JRpm,LawFunctor);
		REGISTER_ATTRIBUTES(LawFunctor,/*nothing here*/);
		DECLARE_LOGGER;	
};
REGISTER_SERIALIZABLE(Law2_Dem3DofGeom_JointedRockPMPhys_JRpm);

/** This class holds information associated with each body */
class JRpmMat: public GranularMat {
	public:
		int exampleNumber; ///<Number of "stone"
		bool initCohesive, isDamaged;
		Real stressCompressMax, Brittleness, G_over_E; ///<Parameters for damaging
		
		/// +++ jointed interactions
		bool onJoint, // is the sphere located on a joint surface?
		     side; // which side of the joint surface?
		Vector3r jointNormal; // normal to the joint plane
		
		JRpmMat(): 
			exampleNumber(0.), 
			initCohesive(false), 
			isDamaged(false), 
			stressCompressMax(0), 
			Brittleness(0), 
			G_over_E(1),
			onJoint(0), side(0) {createIndex(); jointNormal=Vector3r::ZERO;};
		REGISTER_ATTRIBUTES(GranularMat, 
			(exampleNumber)
			(initCohesive)
			(isDamaged)
			(stressCompressMax)
			(Brittleness)
			(G_over_E)
			(onJoint) (side) (jointNormal));
		REGISTER_CLASS_AND_BASE(JRpmMat,GranularMat);
		REGISTER_CLASS_INDEX(JRpmMat,GranularMat);
};
REGISTER_SERIALIZABLE(JRpmMat);


class Ip2_JRpmMat_JRpmMat_JRpmPhys: public InteractionPhysicsFunctor{
	private:
	public:
		Real initDistance;

		Ip2_JRpmMat_JRpmMat_JRpmPhys(){
			initDistance = 0;
		}

		virtual void go(const shared_ptr<Material>& pp1, const shared_ptr<Material>& pp2, const shared_ptr<Interaction>& interaction);
		REGISTER_ATTRIBUTES(InteractionPhysicsFunctor,
			(initDistance)
		);

		FUNCTOR2D(JRpmMat,JRpmMat);
		REGISTER_CLASS_AND_BASE(Ip2_JRpmMat_JRpmMat_JRpmPhys,InteractionPhysicsFunctor);
		DECLARE_LOGGER;
};
REGISTER_SERIALIZABLE(Ip2_JRpmMat_JRpmMat_JRpmPhys);


class JRpmPhys: public NormalShearInteraction {
	private:
	public:
		Real crossSection,E,G,tanFrictionAngle,lengthMaxCompression,lengthMaxTension;
		Real omega, Fn, sigmaN, epsN;
		Vector3r epsT, sigmaT, Fs;
		 
		bool isCohesive;

		/// +++ jointed interactions
		Vector3r modifiedNormal; // to assign the joint normal to the interaction
		bool isOnJoint; // is the interaction intersect by a joint?

		JRpmPhys(): NormalShearInteraction(),
			crossSection(0),
			E(0),
			G(0), 
			tanFrictionAngle(0), 
			lengthMaxCompression(0), 
			lengthMaxTension(0),
			isCohesive(false),
			isOnJoint(false) { createIndex(); epsT=Vector3r::ZERO; omega=0; Fn=0; Fs=Vector3r::ZERO; modifiedNormal=Vector3r::ZERO;}
		virtual ~JRpmPhys();

		REGISTER_ATTRIBUTES(NormalShearInteraction,
			(E)
			(G)
			(tanFrictionAngle)
			(isCohesive)
			(isOnJoint)
			(modifiedNormal)
			(lengthMaxCompression)
			(lengthMaxTension)
		);
	REGISTER_CLASS_AND_BASE(JRpmPhys,NormalShearInteraction);
	REGISTER_CLASS_INDEX(JRpmPhys,NormalShearInteraction);
};
REGISTER_SERIALIZABLE(JRpmPhys);

#include"JointedRockPM.hpp"
#include<yade/core/Scene.hpp>
#include<yade/pkg-dem/DemXDofGeom.hpp>
#include<yade/pkg-dem/Shop.hpp>


YADE_PLUGIN((Law2_Dem3DofGeom_JointedRockPMPhys_JRpm)(JRpmMat)(Ip2_JRpmMat_JRpmMat_JRpmPhys)(JRpmPhys));


/********************** Law2_Dem3DofGeom_JointedRockPMPhys_JRpm ****************************/
CREATE_LOGGER(Law2_Dem3DofGeom_JointedRockPMPhys_JRpm);

void Law2_Dem3DofGeom_JointedRockPMPhys_JRpm::go(shared_ptr<InteractionGeometry>& ig, shared_ptr<InteractionPhysics>& ip, Interaction* contact, Scene* rootBody){
	Dem3DofGeom* geom=static_cast<Dem3DofGeom*>(ig.get());
	JRpmPhys* phys=static_cast<JRpmPhys*>(ip.get());
	
	Real displN=geom->displacementN();
	const Real& crossSection=phys->crossSection;
	
	Real& epsN=phys->epsN;
	Vector3r& epsT=phys->epsT;
	
	Real& sigmaN=phys->sigmaN;
	Vector3r& sigmaT=phys->sigmaT;
	
	const Real& E = phys->E;
	const Real& G = phys->G;
	
	Real& Fn=phys->Fn;
	Vector3r& Fs=phys->Fs;
	

	if(displN<=0){
		/**Normal Interaction*/
		epsN=geom->strainN();
		epsT=geom->strainT();
		
		
		sigmaN=E*epsN;
		Fn=sigmaN*crossSection; /// Fn=E*phys->epsN*(pi*Rmin^2): why phys->epsN and not phys->displN????
		
		phys->normalForce=Fn*geom->normal;
		
		/// to do-> modify the interaction geometry in order to have a smooth interactions along the plane of the joint
// 		normal and shear directions should be oriented along the joint plane in order to have interpenetration only if the contact point move from this plane
// 		-> probably needs an initialization of the new Un and Us in relation with the initial contact point Fn+=Kn*(Un-Un(0))
// 		-> if Us(t=0) is stored, a dilation can probably be take into account (Fn+=(Us-Us(0))*tan(psi))
// 		if (phys->isOnJoint) {
// 			if (geom->normal.Dot(phys->modifiedNormal) > 0) { 
// 			    phys->normalForce=Fn*phys->modifiedNormal;
// 			}
// 			else { 
// 			    phys->normalForce= -Fn*phys->modifiedNormal;
// 			}
//			// maybe an other way to modify the interaction geomatry...
// 		}
// 		else phys->normalForce=Fn*geom->normal;
		

		sigmaT=G*epsT;

		/**Check, whether the shear stress more, than normal force multiplicated to tanFrictionAngle*/
		
		Real maxFsSq = phys->normalForce.SquaredLength()*phys->tanFrictionAngle*phys->tanFrictionAngle;
		
		if(sigmaT.SquaredLength()>maxFsSq) {
			sigmaT*=sqrt(maxFsSq/(sigmaT.SquaredLength()));
		}
		
		Fs=sigmaT*crossSection;
		phys->shearForce = Fs;

		/**Normal Interaction*/
		applyForceAtContactPoint(phys->normalForce + phys->shearForce, geom->contactPoint, contact->getId1(), geom->se31.position, contact->getId2(), geom->se32.position, rootBody);
		
 		/**Destruction*/
		if ((phys->isCohesive)&&(displN<(-phys->lengthMaxCompression))) {
			phys->isCohesive = false;
		}
		return;
	} else {
		/**If spheres do not touch, check, whether they are cohesive*/
		if (phys->isCohesive) {
			/**If the distance between spheres more than critical and they are cohesive,
			 * we delete the interaction Destruction **/
			if (displN>(phys->lengthMaxTension)) {
				rootBody->interactions->requestErase(contact->getId1(),contact->getId2());
				return; 
			} else {
			/**If the distance between spheres less than critical and they are cohesive,
			 * we aply additional forces to keep particles together **/
				phys->normalForce=phys->kn*displN*geom->normal; // phys->displN is used as phys->epsN is not defined for non contacting particles (displN<0 here -> Negative Fn)
				applyForceAtContactPoint(phys->normalForce, geom->contactPoint, contact->getId1(), geom->se31.position, contact->getId2(), geom->se32.position, rootBody);
				return;
			}
		} else {
			/** Delete interactions */ 
			rootBody->interactions->requestErase(contact->getId1(),contact->getId2());
			return;
		}
	}
}


CREATE_LOGGER(Ip2_JRpmMat_JRpmMat_JRpmPhys);

void Ip2_JRpmMat_JRpmMat_JRpmPhys::go(const shared_ptr<Material>& pp1, const shared_ptr<Material>& pp2, const shared_ptr<Interaction>& interaction){
	if(interaction->interactionPhysics) return; 

	Dem3DofGeom* contGeom=YADE_CAST<Dem3DofGeom*>(interaction->interactionGeometry.get());
	Omega& OO=Omega::instance();
	assert(contGeom);
	//LOG_WARN(Omega::instance().getCurrentIteration());
	const shared_ptr<JRpmMat>& rpm1=YADE_PTR_CAST<JRpmMat>(pp1);
	const shared_ptr<JRpmMat>& rpm2=YADE_PTR_CAST<JRpmMat>(pp2);
	
	long cohesiveThresholdIter=10;
	
	bool initCohesive = rpm1->initCohesive*rpm2->initCohesive;
	
	Real E12=2*rpm1->young*rpm2->young/(rpm1->young+rpm2->young);
	Real minRad=(contGeom->refR1<=0?contGeom->refR2:(contGeom->refR2<=0?contGeom->refR1:min(contGeom->refR1,contGeom->refR2)));
	Real S12=Mathr::PI*pow(minRad,2);
	Real G_over_E = (rpm1->G_over_E + rpm2->G_over_E)/2;
	shared_ptr<JRpmPhys> contPhys(new JRpmPhys());
	contPhys->E=E12;
	contPhys->G=E12*G_over_E;
	contPhys->tanFrictionAngle=tan(.5*(rpm1->frictionAngle+rpm2->frictionAngle));
	contPhys->crossSection=S12;
	contPhys->kn=contPhys->E*contPhys->crossSection;
	contPhys->ks=contPhys->G*contPhys->crossSection;
	
	contPhys->lengthMaxCompression=S12*(0.5*(rpm1->stressCompressMax+rpm2->stressCompressMax))/(contPhys->kn);
	contPhys->lengthMaxTension=S12*(0.5*(rpm1->stressCompressMax*rpm1->Brittleness+rpm2->stressCompressMax*rpm2->Brittleness))/(contPhys->kn);
	
	initDistance = contGeom->displacementN();

	if ((rpm1->exampleNumber==rpm2->exampleNumber)&&(initDistance<(contPhys->lengthMaxTension))&&(initCohesive)&&(OO.getCurrentIteration()<=cohesiveThresholdIter)) {
		contPhys->isCohesive=true;
	}
	
	/// +++ Jointed interactions
	if (( rpm1->onJoint && rpm2->onJoint ) && ( rpm1->jointNormal == rpm2->jointNormal ) && ( rpm1->side != rpm2->side )) {
		  contPhys->isCohesive = false;
		  contPhys->isOnJoint = true;
// 		  contPhys->modifiedNormal = rpm1->jointNormal;
		}

	interaction->interactionPhysics=contPhys;
}

JRpmPhys::~JRpmPhys(){};

Attachment: 2joints_Comp.gif
Description: GIF image