yade-dev team mailing list archive

[noreply@xxxxxxxxxxxxx]

------------------------------------------------------------
revno: 2600
committer: Klaus Thoeni <klaus.thoeni@xxxxxxxxx>
branch nick: yade
timestamp: Thu 2010-12-09 18:46:51 +1300
message:
  New implementation of a particle model for wires and rockfall meshes (initial state)
added:
  pkg/dem/WirePM.cpp
  pkg/dem/WirePM.hpp
modified:
  doc/references.bib


--
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

=== modified file 'doc/references.bib'
--- doc/references.bib	2010-11-03 15:50:02 +0000
+++ doc/references.bib	2010-12-09 05:46:51 +0000
@@ -14,6 +14,26 @@
 
 ***************************
 
+@Article{ Bertrand2008,
+	title = "Discrete element method (DEM) numerical modeling of double-twisted hexagonal mesh ",
+	author = "D. Bertrand and F. Nicot and P. Gotteland and S. Lambert",
+	journal = "Canadian Geotechnical Journal",
+	pages = "1104--1117",
+	volume = "45",
+	number = "8",
+	year = "2008"
+}
+
+@article{ Bertrand2005,
+	title = "Modelling a geo-composite cell using discrete analysis",
+	author = "D. Bertrand and F. Nicot and P. Gotteland and S. Lambert",
+	journal = "Computers and Geotechnics",
+	pages = "564--577",
+	volume = "32",
+	number = "8",
+	year = "2005"
+}
+
 @Article{ Chareyre2002a,
 	title = "Theoretical versus experimental modeling of the anchorage capacity of geotextiles in trenches.",
 	author = "B. Chareyre and L. Briancon and P. Villard",

=== added file 'pkg/dem/WirePM.cpp'
--- pkg/dem/WirePM.cpp	1970-01-01 00:00:00 +0000
+++ pkg/dem/WirePM.cpp	2010-12-09 05:46:51 +0000
@@ -0,0 +1,214 @@
+/* Klaus Thoeni 2010  */
+
+#include"WirePM.hpp"
+#include<yade/core/Scene.hpp>
+#include<yade/pkg/dem/ScGeom.hpp>
+#include<yade/core/Omega.hpp>
+
+YADE_PLUGIN((WireMat)(WireState)(WirePhys)(Ip2_WireMat_WireMat_WirePhys)(Law2_ScGeom_WirePhys_WirePM));
+
+
+/********************** WireMat ****************************/
+CREATE_LOGGER(WireMat);
+
+void WireMat::postLoad(WireMat&){
+	if(StrainStressValues.size() < 2)
+		throw invalid_argument("WireMat.StrainStressValues: at least two points must be given.");
+	if(StrainStressValues[0](0) == 0. && StrainStressValues[0](1) == 0.)
+		throw invalid_argument("WireMat.StrainStressValues: Definition must start with values greather then zero (strain>0,stress>0)");
+	// compute cross-sectin area
+	As = pow(diameter*0.5,2)*Mathr::PI;
+
+	// debug printout
+// 	Vector2r vbegin = *StrainStressValues.begin();
+// // 	Vector2r vbegin = StrainStressValues.front();
+// 	std::cerr << vbegin(0) << endl;
+// 	std::cerr << vbegin(1) << endl;
+// 	std::cerr << endl;
+// 
+// 	Real v0 = (*(StrainStressValues.begin()))(0);
+// 	Real v1 = (*StrainStressValues.begin())(1);
+// 	std::cerr << v0 << endl;
+// 	std::cerr << v1 << endl;
+
+	//BUG ????? postLoad is called twice,
+}
+
+
+/********************** Law2_ScGeom_WirePhys_WirePM ****************************/
+CREATE_LOGGER(Law2_ScGeom_WirePhys_WirePM);
+
+void Law2_ScGeom_WirePhys_WirePM::go(shared_ptr<IGeom>& ig, shared_ptr<IPhys>& ip, Interaction* contact){
+	ScGeom* geom = static_cast<ScGeom*>(ig.get()); 
+	WirePhys* phys = static_cast<WirePhys*>(ip.get());
+	const int &id1 = contact->getId1();
+	const int &id2 = contact->getId2();
+	Body* b1 = Body::byId(id1,scene).get();
+	Body* b2 = Body::byId(id2,scene).get();
+	
+	Real displN = geom->penetrationDepth; // NOTE: ScGeom->penetrationDepth>0 when spheres interpenetrate, and therefore, for wire always negative
+
+	/* get reference to values since values are updated for unloading */
+	vector<Vector2r> &DFValues = phys->DisplForceValues;
+	vector<Real> &kValues = phys->StiffnessValues;
+
+	Real D = displN - phys->initD; // interparticular distance is computed depending on the equilibrium distance
+
+	/* check whether the particles are linked or not */
+	if ( !phys->isLinked ) { // destroy the interaction before calculation
+		scene->interactions->requestErase(contact->getId1(),contact->getId2());
+		return;
+	}
+	if ( (phys->isLinked) && (D < DFValues.back()(0)) ) { // spheres are linked but failure because of reaching maximal admissible displacement 
+		phys->isLinked=false; 
+		// update body state with the number of broken links
+		WireState* st1=dynamic_cast<WireState*>(b1->state.get());
+		WireState* st2=dynamic_cast<WireState*>(b2->state.get());
+		st1->numBrokenLinks+=1;
+		st2->numBrokenLinks+=1;
+		scene->interactions->requestErase(contact->getId1(),contact->getId2());
+		return;
+	}
+	
+	/* compute normal force Fn */
+	Real Fn = 0.;
+	if ( D > DFValues[0](0) )
+		Fn = kValues[0] * (D-phys->Dplast);
+	else {
+		bool isDone = false;
+		unsigned int i = 0;
+		while (!isDone && i < DFValues.size()) { 
+			i++;
+			if ( D > DFValues[i](0) ) {
+				Fn = DFValues[i-1](1) + (D-DFValues[i-1](0))*kValues[i];
+				phys->Dplast = D - Fn/kValues[0];
+				isDone = true;
+				// update values for unloading
+				DFValues[0](0) = D;
+				DFValues[0](1) = Fn;
+			}
+		}
+	}
+	
+	TRVAR3( displN, D, Fn );
+
+	/* compression forces cannot be applied to wires */
+	if (Fn > 0.) Fn = 0.;
+	
+	phys->normalForce = Fn*geom->normal; // NOTE: normal is position2-position1 - It is directed from particle1 to particle2
+	        
+	State* st1 = Body::byId(id1,scene)->state.get();
+	State* st2 = Body::byId(id2,scene)->state.get();
+	
+	/* apply forces */
+	Vector3r f = phys->normalForce;
+	// these lines to adapt to periodic boundary conditions
+	if ( !scene->isPeriodic )  
+		applyForceAtContactPoint(f , geom->contactPoint , id2, st2->se3.position, id1, st1->se3.position);
+	else { // in scg we do not wrap particles positions, hence "applyForceAtContactPoint" cannot be used when scene is periodic
+		scene->forces.addForce(id1,-f);
+		scene->forces.addForce(id2,f);
+	}
+
+}
+
+CREATE_LOGGER(Ip2_WireMat_WireMat_WirePhys);
+
+void Ip2_WireMat_WireMat_WirePhys::go(const shared_ptr<Material>& b1, const shared_ptr<Material>& b2, const shared_ptr<Interaction>& interaction){
+	
+	/* avoid any updates if interactions which already exist */
+	if(interaction->phys) return; 
+	
+	ScGeom* geom=dynamic_cast<ScGeom*>(interaction->geom.get());
+	assert(geom);
+
+	/* set equilibrium distance, e.g. initial distance between particle (stress free state) */
+	shared_ptr<WirePhys> contactPhysics(new WirePhys()); 
+	contactPhysics->initD = geom->penetrationDepth;
+	contactPhysics->normalForce = Vector3r::Zero();
+
+	/* get values from material */
+	const shared_ptr<WireMat>& mat1 = YADE_PTR_CAST<WireMat>(b1);
+	const shared_ptr<WireMat>& mat2 = YADE_PTR_CAST<WireMat>(b2);
+
+	Real crossSection;
+	vector<Vector2r> SSValues;
+	
+	/* ckeck properties of interaction */
+	if ( mat1->id == mat2->id ) { // interaction of two bodies of the same material
+		crossSection = mat1->As;
+		SSValues = mat1->StrainStressValues;
+		if ( (mat1->isDoubleTwist) && (abs(interaction->getId1()-interaction->getId2())==1) ) // bodies which id differs by 1 are double twisted
+			contactPhysics->isDoubleTwist = true;
+		else
+			contactPhysics->isDoubleTwist = false;
+	}
+	else { // interaction of two bodies of two different materials, take weaker material and no double-twist
+		contactPhysics->isDoubleTwist = false;
+		if ( mat1->diameter <= mat2->diameter){
+			crossSection = mat1->As;
+			SSValues = mat1->StrainStressValues;
+		}
+		else {
+			crossSection = mat2->As;
+			SSValues = mat2->StrainStressValues;
+		}
+	}
+	
+	Real R1 = geom->radius1;
+	Real R2 = geom->radius2;
+	
+	Real l0 = R1 + R2 - contactPhysics->initD; // initial lenght of the wire (can be single or double twisted)
+
+	/* compute displscement-force values (tension negative since ScGem is used!) */
+	vector<Vector2r> DFValues;
+	for ( vector<Vector2r>::iterator it = SSValues.begin(); it != SSValues.end(); it++ ) {
+		Vector2r values = Vector2r::Zero();
+		values(0) = -(*it)(0)*l0;
+		values(1) = -(*it)(1)*crossSection;
+		DFValues.push_back(values);
+	}
+
+	/* compute elastic stiffness of single wire */
+	vector<Real> kValues;
+	Real k = DFValues[0](1) / DFValues[0](0);
+
+	/* update values if the interaction is double twiseted */
+	if ( contactPhysics->isDoubleTwist ) {
+		Real alpha = atan( l0 / (3.*Mathr::PI*mat1->diameter) );
+		Real kh = k * ( l0*mat1->diameter/crossSection ) / ( 48.*cos(alpha) * ( 41./9.*(1.+mat1->poisson) + 17./4.*pow(tan(alpha),2) ) );
+		k = 2. * ( mat1->lambdak*kh + (1-mat1->lambdak)*k );
+		Real F = k * DFValues[0](0);
+		Real mappingF = F/DFValues[0](1);
+		DFValues[0](1) = F;
+		for (unsigned int i = 1; i<DFValues.size(); i++) {
+			DFValues[i](0) *= mat1->lambdaEps;
+			DFValues[i](1) *= mappingF;
+		}
+	}
+	
+	/* store displscement-force values in physics */
+	contactPhysics->DisplForceValues = DFValues;
+
+	/* compute stiffness-values of wire */
+	contactPhysics->kn = k;
+	kValues.push_back(k);
+	for( unsigned int i = 1 ; i < DFValues.size()-1; i++ ) {
+		Real deltau = -DFValues[i](0) + DFValues[i-1](0);
+		Real deltaF = -DFValues[i](1) + DFValues[i-1](1);
+		k = deltaF/deltau;
+		kValues.push_back(k);
+	}
+	contactPhysics->StiffnessValues = kValues;
+
+	/* set particles as linked */
+	if ( (scene->iter < linkThresholdIteration))
+		contactPhysics->isLinked=true;
+	else
+		contactPhysics->isLinked=false;
+
+	interaction->phys = contactPhysics;
+
+}
+
+WirePhys::~WirePhys(){}

=== added file 'pkg/dem/WirePM.hpp'
--- pkg/dem/WirePM.hpp	1970-01-01 00:00:00 +0000
+++ pkg/dem/WirePM.hpp	2010-12-09 05:46:51 +0000
@@ -0,0 +1,106 @@
+/* Klaus Thoeni 2010 */
+
+/**
+=== OVERVIEW OF WirePM ===
+
+A particle model to simulate single wires and rockfall meshes (see Bertrad et al. 2005, Bertrad et al. 2008).
+
+Features of the interaction law:
+
+1. The law is designed for particles which do not touch. A link will be created by defining an interaction radius and running the threshold iteration.
+
+2. The contact law is for tension only. Compressive forces never exist between the particles. However, elastic unloading is considered.
+
+3. The force displacement curve which defines the interaction forces is piecewise linear and defined by the stress-strain curve of the wire material. Any piecewise linear curve can be used. 
+
+Remarks:
+- The contact law is new and still needs some testing :-) 
+*/
+
+#pragma once
+
+#include<yade/pkg/common/ElastMat.hpp>
+#include<yade/pkg/common/Dispatching.hpp>
+#include<yade/pkg/common/NormShearPhys.hpp>
+#include<yade/pkg/dem/ScGeom.hpp>
+
+/** This class holds information associated with each body state*/
+class WireState: public State {
+	YADE_CLASS_BASE_DOC_ATTRS_CTOR(WireState,State,"Wire state information of each body.\n\nNone of that is used for computation (at least not now), only for post-processing.",
+		((int,numBrokenLinks,0,,"Number of broken links (e.g. number of wires connected to the body which are broken). [-]"))
+		,
+		createIndex();
+	);
+	REGISTER_CLASS_INDEX(WireState,State);
+};
+REGISTER_SERIALIZABLE(WireState);
+
+/** This class holds information associated with each body */
+class WireMat: public ElastMat {
+	public:
+		virtual shared_ptr<State> newAssocState() const { return shared_ptr<State>(new WireState); }
+		virtual bool stateTypeOk(State* s) const { return (bool)dynamic_cast<WireState*>(s); }
+		void postLoad(WireMat&);
+	DECLARE_LOGGER;
+	YADE_CLASS_BASE_DOC_ATTRS_CTOR(WireMat,ElastMat,"Material for use with the Wire classes",
+		((Real,diameter,0.0027,," (Diameter of the single wire in [m] (the diameter is used to compute the cross-section area of the wire)."))
+		((vector<Vector2r>,StrainStressValues,,Attr::triggerPostLoad,"Piecewise linear definition of the stress-strain curve by set of points (strain[-]>0,stress[Pa]>0) for one single wire. Tension only is considered and the point (0,0) is not needed!"))
+		((bool,isDoubleTwist,false,,"Type of the mesh. If true two particles of the same material which body ids differ by one will be considered as double-twisted interaction."))
+		((Real,lambdaEps,0.4,,"Parameter between 0 and 1 to reduce the failure strain of the double-twisted wire (as used by [Bertrand2008]_). [-]"))
+		((Real,lambdak,0.21,,"Parameter between 0 and 1 to compute the elastic stiffness of the double-twisted wire (as used by [Bertrand2008]_): $\\k^D=2(\\lambda_k k_h + (1-\\lambda_k)k^S). [-]"))
+		((Real,As,0,Attr::readonly,"Cross-section area of a single wire used for the computation of the limit normal contact forces. [m²]"))
+		,
+		createIndex();
+	);
+	REGISTER_CLASS_INDEX(WireMat,ElastMat);
+};
+REGISTER_SERIALIZABLE(WireMat);
+
+/** This class holds information associated with each interaction */
+class WirePhys: public NormPhys {
+	public:
+		virtual ~WirePhys();
+	
+		YADE_CLASS_BASE_DOC_ATTRS_CTOR_PY(WirePhys,NormPhys,"Representation of a single interaction of the WirePM type, storage for relevant parameters",
+			((Real,initD,0,,"Equilibrium distance for particles. Computed as the initial inter-particular distance when particle are linked."))
+			((bool,isLinked,false,,"If true particles are linked and will interact. Interactions are linked automatically by the definition of the corresponding interaction radius. The value is false if the wire breaks (no more interaction)."))
+			((bool,isDoubleTwist,false,,"If true the properties of the interaction will be defined as a double-twisted wire."))
+			((vector<Vector2r>,DisplForceValues,,(Attr::readonly|Attr::noSave),"Defines the values for force-displacement curve."))
+			((vector<Real>,StiffnessValues,,(Attr::readonly|Attr::noSave),"Defines the values for the different stiffness (first value corresponds to elastic stiffness kn)."))
+			((Real,Dplast,0,Attr::readonly,"Plastic part of the inter-particular distance of the previous step. \n\n.. note::\n\t Only elastic displacements are reversible (the elastic stiffness is used for unloading) and compressive forces are inadmissible. The compressive stiffness is assumed to be equal to zero (see [Bertrand2005]_).\n\n.."))
+			,
+			createIndex();
+			,
+		);
+	DECLARE_LOGGER;
+	REGISTER_CLASS_INDEX(WirePhys,NormShearPhys);
+};
+REGISTER_SERIALIZABLE(WirePhys);
+
+/** 2d functor creating IPhys (Ip2) taking WireMat and WireMat of 2 bodies, returning type WirePhys */
+class Ip2_WireMat_WireMat_WirePhys: public IPhysFunctor{
+	public:
+		virtual void go(const shared_ptr<Material>& pp1, const shared_ptr<Material>& pp2, const shared_ptr<Interaction>& interaction);
+		
+		FUNCTOR2D(WireMat,WireMat);
+		DECLARE_LOGGER;
+		
+		YADE_CLASS_BASE_DOC_ATTRS(Ip2_WireMat_WireMat_WirePhys,IPhysFunctor,"Converts 2 :yref:`path<WireMat.path>` instances to :yref:`path<WirePhys.path>` with corresponding parameters.",
+			((int,linkThresholdIteration,1,,"Iteration to create the link."))
+		);
+};
+REGISTER_SERIALIZABLE(Ip2_WireMat_WireMat_WirePhys);
+
+/** 2d functor creating the interaction law (Law2) based on SphereContactGeometry (ScGeom) and WirePhys of 2 bodies, returning type WirePM */
+class Law2_ScGeom_WirePhys_WirePM: public LawFunctor{
+	public:
+		virtual void go(shared_ptr<IGeom>& _geom, shared_ptr<IPhys>& _phys, Interaction* I);
+		
+		FUNCTOR2D(ScGeom,WirePhys);
+
+		YADE_CLASS_BASE_DOC_ATTRS(Law2_ScGeom_WirePhys_WirePM,LawFunctor,"Constitutive law for the wire model.",
+			((int,linkThresholdIteration,1,,"Iteration to create the link."))
+		);
+	DECLARE_LOGGER;
+};
+REGISTER_SERIALIZABLE(Law2_ScGeom_WirePhys_WirePM);