------------------------------------------------------------ revno: 3840 committer: Jan Stransky <jan.stran...@fsv.cvut.cz> 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.stran...@fsv.cvut.cz> +#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.stran...@fsv.cvut.cz> + +#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);
_______________________________________________ Mailing list: https://launchpad.net/~yade-dev Post to : yade-dev@lists.launchpad.net Unsubscribe : https://launchpad.net/~yade-dev More help : https://help.launchpad.net/ListHelp