yade-dev team mailing list archive

[noreply@xxxxxxxxxxxxx]

------------------------------------------------------------
revno: 3840
committer: Jan Stransky <jan.stransky@xxxxxxxxxxx>
timestamp: Wed 2016-04-13 18:45:33 +0200
message:
  Added a material model for mortar layer
  
  added one reference
  added example script
added:
  examples/mortar/
  examples/mortar/matModel.py
  pkg/dem/MortarMat.cpp
  pkg/dem/MortarMat.hpp
modified:
  doc/references.bib


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

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

=== modified file 'doc/references.bib'
--- doc/references.bib	2016-03-24 22:08:26 +0000
+++ doc/references.bib	2016-04-13 16:45:33 +0000
@@ -892,3 +892,10 @@
   url = {http://docsdrive.com/pdfs/sciencepublications/jmssp/2005/8-11.pdf}
 }
 
+@misc{Lourenco1994,
+	author = {P B Lourenço},
+	title = {Analysis of masonry structures with interface elements. theory and applications},
+	year = {1994},
+	note = {Report No. 03-21-22-0-01, Delft University of Technology, Faculty of Civil Engineering},
+	url = {http://www.csarmento.uminho.pt/docs/ncr/de_civil/1994_Lourenco.pdf}
+}

=== added directory 'examples/mortar'
=== added file 'examples/mortar/matModel.py'
--- examples/mortar/matModel.py	1970-01-01 00:00:00 +0000
+++ examples/mortar/matModel.py	2016-04-13 16:45:33 +0000
@@ -0,0 +1,47 @@
+from yade import plot
+
+dt = 1e-5
+
+mortar = MortarMat()
+polyMat = PolyhedraMat()
+for mat in (mortar,polyMat):
+	O.materials.append(mat)
+
+def sim(angle):
+	O.reset()
+	O.dt = dt
+	bs = b1,b2 = [polyhedron(((-1,-1,-1),(+1,-1,-1),(-1,+1,-1),(+1,+1,-1),(-1,-1,+1),(+1,-1,+1),(-1,+1,+1),(+1,+1,+1)),material=mortar) for i in (0,1)]
+	b2.state.pos = (0,0,2)
+	for b in bs:
+		b.state.blockedDOFs = 'xyzXYZ'
+	O.bodies.append(bs)
+	#
+	factor=1.1
+	O.engines=[
+		ForceResetter(),
+		InsertionSortCollider([Bo1_Polyhedra_Aabb(aabbEnlargeFactor=factor,label='bo1')]),
+		InteractionLoop(
+			[Ig2_Polyhedra_Polyhedra_PolyhedraGeomOrScGeom(label='ig2')], 
+			[Ip2_PolyhedraMat_PolyhedraMat_PolyhedraPhys(),Ip2_MortarMat_MortarMat_MortarPhys()],
+			[Law2_PolyhedraGeom_PolyhedraPhys_Volumetric(),Law2_ScGeom_MortarPhys_Lourenco()]
+		),
+		NewtonIntegrator(),
+	]
+	ig2.ig2scGeom.interactionDetectionFactor = factor
+	O.step()
+	ig2.ig2scGeom.interactionDetectionFactor = bo1.aabbEnlargeFactor = 1
+	b2.state.vel = (sin(angle),0,cos(angle))
+	while len([i for i in O.interactions]):
+		sn,st = i.phys.sigmaN, i.phys.sigmaT.norm()
+		O.step()
+		if O.iter > 1e6:
+			raise RuntimeError, "TODO"
+	plot.addData(sn=sn,st=st)
+
+n = 50
+for i in xrange(n):
+	sim(i*pi/(n-1))
+
+plot.plots = {'sn':'st'}
+plot.matplotlib.pyplot.axes().set_aspect(1)
+p = plot.plot()

=== added file 'pkg/dem/MortarMat.cpp'
--- pkg/dem/MortarMat.cpp	1970-01-01 00:00:00 +0000
+++ pkg/dem/MortarMat.cpp	2016-04-13 16:45:33 +0000
@@ -0,0 +1,127 @@
+// 2016 © Jan Stránský <jan.stransky@xxxxxxxxxxx> 
+#include"MortarMat.hpp"
+
+
+YADE_PLUGIN((MortarMat)(Ip2_MortarMat_MortarMat_MortarPhys)(MortarPhys)(Law2_ScGeom_MortarPhys_Lourenco))
+
+
+CREATE_LOGGER(Ip2_MortarMat_MortarMat_MortarPhys);
+void Ip2_MortarMat_MortarMat_MortarPhys::go(const shared_ptr<Material>& material1, const shared_ptr<Material>& material2, const shared_ptr<Interaction>& interaction){
+	if (interaction->phys) return;
+	if (scene->iter >= cohesiveThresholdIter) return;
+	shared_ptr<MortarPhys> phys(new MortarPhys());
+	interaction->phys = phys;
+	MortarMat* mat1 = YADE_CAST<MortarMat*>(material1.get());
+	MortarMat* mat2 = YADE_CAST<MortarMat*>(material2.get());
+
+	if (mat1->id>=0 && mat1->id == mat2->id) {
+		#define _CPATTR(a) phys->a=mat1->a
+			_CPATTR(tensileStrength);
+			_CPATTR(compressiveStrength);
+			_CPATTR(cohesion);
+			_CPATTR(ellAspect);
+		#undef _CPATTR
+		phys->tangensOfFrictionAngle = std::tan(mat1->frictionAngle);
+	} else {
+		// averaging over both materials
+		#define _MINATTR(a) phys->a=std::min(mat1->a,mat2->a)
+		#define _AVGATTR(a) phys->a=.5*(mat1->a+mat2->a)
+			_MINATTR(tensileStrength);
+			_MINATTR(compressiveStrength);
+			_MINATTR(cohesion);
+			_AVGATTR(ellAspect);
+		#undef _AVGATTR
+		#undef _MINATTR
+		phys->tangensOfFrictionAngle = std::tan(.5*(mat1->frictionAngle+mat2->frictionAngle));
+		// E, G, kn, ks, crosssection, refPD, refLength to be computed in Law2
+	}
+}
+
+
+
+
+CREATE_LOGGER(MortarPhys);
+MortarPhys::~MortarPhys(){};
+
+
+
+
+
+
+
+/********************** Law2_ScGeom_MortarPhys_Lourenco ****************************/
+CREATE_LOGGER(Law2_ScGeom_MortarPhys_Lourenco);
+bool Law2_ScGeom_MortarPhys_Lourenco::go(shared_ptr<IGeom>& iGeom, shared_ptr<IPhys>& iPhys, Interaction* interaction){
+	ScGeom* geom=static_cast<ScGeom*>(iGeom.get());
+	MortarPhys* phys=static_cast<MortarPhys*>(iPhys.get());
+
+	Body::id_t id1 = interaction->getId1();
+ 	Body::id_t id2 = interaction->getId2();
+	const shared_ptr<Body> b1 = Body::byId(id1,scene);
+	const shared_ptr<Body> b2 = Body::byId(id2,scene);
+
+	/* just the first time */
+	if (interaction->isFresh(scene)) {
+		const Vector3r& pos1 = b1->state->pos;
+		const Vector3r& pos2 = b2->state->pos;
+		const Real& r1 = geom->refR1;
+		const Real& r2 = geom->refR2;
+		Vector3r shift2 = scene->isPeriodic? Vector3r(scene->cell->hSize*interaction->cellDist.cast<Real>()) : Vector3r::Zero();
+		phys->refLength = (pos2 - pos1 + shift2).norm();
+		Real minRad = r1 <= 0 ? r2 : r2 <= 0 ? r1 : std::min(r1,r2);
+		phys->crossSection = std::pow(minRad,2);
+		phys->refPD = geom->refR1 + geom->refR2 - phys->refLength;
+		const shared_ptr<MortarMat> mat1 = YADE_PTR_CAST<MortarMat>(b1->material);
+		const shared_ptr<MortarMat> mat2 = YADE_PTR_CAST<MortarMat>(b2->material);
+		const Real& E1(mat1->young);
+		const Real& E2(mat2->young);
+		const Real& n1(mat1->poisson);
+		const Real& n2(mat2->poisson);
+		phys->kn = 2*E1*r1*E2*r2/(E1*r1+E2*r2);
+		phys->ks = 2*E1*r1*n1*E2*r2*n2/(E1*r1*n1+E2*r2*n2);
+		phys->E = phys->kn * phys->refLength / phys->crossSection;
+		phys->G = phys->ks * phys->refLength / phys->crossSection;
+	}
+	
+	/* shorthands */
+	Real& epsN(phys->epsN);
+	Vector3r& epsT(phys->epsT);
+	const Real& E(phys->E); \
+	const Real& G(phys->G);
+	const Real& crossSection(phys->crossSection);
+	Real& sigmaN(phys->sigmaN);
+	Vector3r& sigmaT(phys->sigmaT);
+
+	epsN = - (-phys->refPD + geom->penetrationDepth) / phys->refLength;
+	geom->rotate(epsT);
+	epsT += geom->shearIncrement() / phys->refLength;
+
+	/* constitutive law */
+	sigmaN = E*epsN;
+	sigmaT = G*epsT;
+
+	Real st = sigmaT.norm();
+	bool cond1 = sigmaN - phys->tensileStrength > 0;
+	bool cond2 = st + sigmaN*phys->tangensOfFrictionAngle - phys->cohesion > 0;
+	bool cond3 = std::pow(sigmaN,2) + std::pow(phys->ellAspect*st,2) - std::pow(phys->compressiveStrength,2) > 0;
+	if (cond1 || cond2 || cond3) {
+		return false;
+	}
+   
+	phys->normalForce = -sigmaN*crossSection*geom->normal;
+	phys->shearForce = -sigmaT*crossSection;
+
+	State* s1 = b1->state.get();
+	State* s2 = b2->state.get();	
+
+	Vector3r f = -phys->normalForce - phys->shearForce;
+	if (!scene->isPeriodic) {
+		applyForceAtContactPoint(f, geom->contactPoint , id1, s1->se3.position, id2, s2->se3.position);
+	} else {
+		scene->forces.addForce(id1,f);
+		scene->forces.addForce(id2,-f);
+		scene->forces.addTorque(id1,(geom->radius1+.5*(phys->refPD-geom->penetrationDepth))*geom->normal.cross(f));
+		scene->forces.addTorque(id2,(geom->radius2+.5*(phys->refPD-geom->penetrationDepth))*geom->normal.cross(f));
+	}
+	return true;
+}

=== added file 'pkg/dem/MortarMat.hpp'
--- pkg/dem/MortarMat.hpp	1970-01-01 00:00:00 +0000
+++ pkg/dem/MortarMat.hpp	2016-04-13 16:45:33 +0000
@@ -0,0 +1,95 @@
+// 2016 © Jan Stránský <jan.stransky@xxxxxxxxxxx> 
+
+#pragma once
+
+#include<pkg/common/ElastMat.hpp>
+#include<pkg/common/Dispatching.hpp>
+#include<pkg/common/Sphere.hpp>
+#include<pkg/common/PeriodicEngines.hpp>
+#include<pkg/common/NormShearPhys.hpp>
+#include<pkg/dem/DemXDofGeom.hpp>
+#include<pkg/dem/ScGeom.hpp>
+#include<pkg/dem/FrictPhys.hpp>
+
+namespace py=boost::python;
+
+
+
+class MortarMat: public FrictMat {
+	public:
+		YADE_CLASS_BASE_DOC_ATTRS_CTOR(MortarMat,FrictMat,"Material for mortar interface, used in Ip2_MortarMat_MortarMat_MortarPhys and Law2_ScGeom_MortarPhys_Lourenco. Default values according to ",
+			((Real,young,1e9,,"Normal elastic modulus [Pa]"))
+			((Real,poisson,1,,"Shear to normal modulus ratio"))
+			((Real,frictionAngle,.25,,"Friction angle"))
+			//
+			((Real,tensileStrength,1e6,,"tensileStrength [Pa]"))
+			((Real,compressiveStrength,10e6,,"compressiveStrength [Pa]"))
+			((Real,cohesion,1e6,,"cohesion [Pa]"))
+			((Real,ellAspect,3,,"aspect ratio of elliptical 'cap'. Value >1 means the ellipse is longer along normal stress axis."))
+			,
+			createIndex();
+		);
+		REGISTER_CLASS_INDEX(MortarMat,FrictMat);
+};
+REGISTER_SERIALIZABLE(MortarMat);
+
+
+
+
+
+
+class MortarPhys: public FrictPhys {
+	public:
+		Real epsN, sigmaN;
+		Vector3r epsT, sigmaT;
+		virtual ~MortarPhys();
+		YADE_CLASS_BASE_DOC_ATTRS_CTOR_PY(MortarPhys,FrictPhys,"IPhys class containing parameters of MortarMat. Used by Law2_ScGeom_MortarPhys_Lourenco.",
+			((Real,tensileStrength,NaN,,"tensileStrength [Pa]"))
+			((Real,compressiveStrength,NaN,,"compressiveStrength [Pa]"))
+			((Real,cohesion,NaN,,"cohesion [Pa]"))
+			((Real,ellAspect,NaN,,"aspect ratio of elliptical 'cap'. Value >1 means the ellipse is longer along normal stress axis."))
+			((Real,crossSection,NaN,,"Crosssection of interaction"))
+			((Real,refLength,NaN,,"Initial length of interaction [m]"))
+			((Real,refPD,NaN,,"Reference penetration depth [m]"))
+			((Real,E,NaN,,"interaction Young's modulus [Pa]"))
+			((Real,G,NaN,,"interaction shear modulus [Pa]"))
+			, // ctors
+			epsT = Vector3r::Zero();
+			createIndex();
+			,
+			.def_readonly("epsN",&MortarPhys::epsN,"Current normal strain |yupdate|")
+			.def_readonly("epsT",&MortarPhys::epsT,"Current shear strain |yupdate|")
+			.def_readonly("sigmaN",&MortarPhys::sigmaN,"Current normal stress |yupdate|")
+			.def_readonly("sigmaT",&MortarPhys::sigmaT,"Current shear stress |yupdate|")
+		);
+		DECLARE_LOGGER;
+		REGISTER_CLASS_INDEX(MortarPhys,NormShearPhys);
+};
+REGISTER_SERIALIZABLE(MortarPhys);
+
+
+
+
+
+class Ip2_MortarMat_MortarMat_MortarPhys: public IPhysFunctor{
+	public:
+		virtual void go(const shared_ptr<Material>& material1, const shared_ptr<Material>& material2, const shared_ptr<Interaction>& interaction);
+		FUNCTOR2D(MortarMat,MortarMat);
+		DECLARE_LOGGER;
+		YADE_CLASS_BASE_DOC_ATTRS(Ip2_MortarMat_MortarMat_MortarPhys,IPhysFunctor,"Ip2 creating MortarPhys from two MortarMat instances.",
+			((long,cohesiveThresholdIter,2,,"Should new contacts be cohesive? They will before this iter#, they will not be afterwards. If <=0, they will never be."))
+		);
+};
+REGISTER_SERIALIZABLE(Ip2_MortarMat_MortarMat_MortarPhys);
+
+
+
+class Law2_ScGeom_MortarPhys_Lourenco: public LawFunctor{
+	public:
+		bool go(shared_ptr<IGeom>& iGeom, shared_ptr<IPhys>& iPhys, Interaction* interaction);
+		FUNCTOR2D(GenericSpheresContact,MortarPhys);
+		YADE_CLASS_BASE_DOC(Law2_ScGeom_MortarPhys_Lourenco,LawFunctor,"Material law for mortar layer according to [Lourenco1994]_. The contact behaves elastic until brittle failure when reaching strength envelope. The envelope has three parts.\n\nTensile with condition $\\sigma_N-f_t$.\n\nShear part with Mohr-Coulomb condition $|\\sigma_T|+\\sigma_N\\tan\\varphi-c$.\n\nCompressive part with condition $\\sigma_N^2+A^2\\sigma_T^2-f_c^2$\n\nThe main idea is to begin simulation with this model and when the contact is broken, to use standard non-cohesive Law2_PolyhedraGeom_PolyhedraPhys_Volumetric."
+		);
+		DECLARE_LOGGER;
+};
+REGISTER_SERIALIZABLE(Law2_ScGeom_MortarPhys_Lourenco);