------------------------------------------------------------ revno: 2049 committer: Chiara Modenese <c...@engs-018373> branch nick: trunk timestamp: Tue 2010-02-23 14:22:35 +0000 message: Added constitutive law with Hertz formulation for contact stiffnesses added: pkg/dem/Engine/GlobalEngine/HertzMindlin.cpp pkg/dem/Engine/GlobalEngine/HertzMindlin.hpp scripts/mindlin_test.py
-- 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/HertzMindlin.cpp' --- pkg/dem/Engine/GlobalEngine/HertzMindlin.cpp 1970-01-01 00:00:00 +0000 +++ pkg/dem/Engine/GlobalEngine/HertzMindlin.cpp 2010-02-23 14:22:35 +0000 @@ -0,0 +1,120 @@ +// 2010 © Chiara Modenese <[email protected]> + + +#include"HertzMindlin.hpp" +#include<yade/pkg-dem/ScGeom.hpp> +#include<yade/core/Omega.hpp> +#include<yade/core/Scene.hpp> + + +YADE_PLUGIN((MindlinPhys)(Ip2_FrictMat_FrictMat_MindlinPhys)(Law2_ScGeom_MindlinPhys_Mindlin)); + + +/******************** MindlinPhys *****************************/ +CREATE_LOGGER(MindlinPhys); + +MindlinPhys::~MindlinPhys(){}; // destructor + + +/******************** Ip2_FrictMat_FrictMat_MindlinPhys *******/ +CREATE_LOGGER(Ip2_FrictMat_FrictMat_MindlinPhys); + +void Ip2_FrictMat_FrictMat_MindlinPhys::go(const shared_ptr<Material>& b1,const shared_ptr<Material>& b2, const shared_ptr<Interaction>& interaction){ + if(interaction->interactionPhysics) return; // no updates of an already existing contact necessary + shared_ptr<MindlinPhys> mindlinPhys(new MindlinPhys()); + interaction->interactionPhysics = mindlinPhys; + FrictMat* mat1 = YADE_CAST<FrictMat*>(b1.get()); + FrictMat* mat2 = YADE_CAST<FrictMat*>(b2.get()); + + + /* from interaction physics */ + Real Ea = mat1->young; + Real Eb = mat2->young; + Real Va = mat1->poisson; + Real Vb = mat2->poisson; + Real fa = mat1->frictionAngle; + Real fb = mat2->frictionAngle; + + + /* from interaction geometry */ + ScGeom* scg = YADE_CAST<ScGeom*>(interaction->interactionGeometry.get()); + Real Da = scg->radius1; + Real Db = scg->radius2; + Vector3r normal=scg->normal; + + + /* calculate stiffness coefficients */ + Real Ga = Ea/(2*(1+Va)); + Real Gb = Eb/(2*(1+Vb)); + Real G = (Ga+Gb)/2; // average of shear modulus + Real V = (Va+Vb)/2; // average of poisson's ratio + Real R = 2*Da*Db/(Da+Db); // harmonic average of the radii + Real Kno = 2*G*Mathr::Sqrt(2*R)/(3*(1-V)); // coefficient for normal stiffness + Real Kso = 2*std::pow((std::pow(G,2)*3*(1-V)*R),1/3)/(2-V); // coefficient for shear stiffness + Real frictionAngle = std::min(fa,fb); + + + /* pass values calculated from above to MindlinPhys */ + mindlinPhys->tangensOfFrictionAngle = std::tan(frictionAngle); + mindlinPhys->prevNormal = scg->normal; + mindlinPhys->kno = Kno; // this is just a coeff + mindlinPhys->kso = Kso; // this is just a coeff +} + + +/******************** Law2_ScGeom_MindlinPhys_Mindlin *********/ +CREATE_LOGGER(Law2_ScGeom_MindlinPhys_Mindlin); + +void Law2_ScGeom_MindlinPhys_Mindlin::go(shared_ptr<InteractionGeometry>& ig, shared_ptr<InteractionPhysics>& ip, Interaction* contact, Scene* ncb){ + Real dt = Omega::instance().getTimeStep(); // get time step + + int id1 = contact->getId1(); // get id body 1 + int id2 = contact->getId2(); // get id body 2 + + State* de1 = Body::byId(id1,ncb)->state.get(); + State* de2 = Body::byId(id2,ncb)->state.get(); + + ScGeom* scg = static_cast<ScGeom*>(ig.get()); + MindlinPhys* phys = static_cast<MindlinPhys*>(ip.get()); + + + /*** NORMAL FORCE ***/ + Real uN = scg->penetrationDepth; // get overlapping + if (uN<0) {ncb->interactions->requestErase(contact->getId1(),contact->getId2()); return;} + /*** Hertz-Mindlin's formulation (PFC) ***/ + phys->kn = phys->kno*Mathr::Sqrt(uN); // normal stiffness + Real Fn = phys->kn*uN; // normal Force (scalar) + phys->normalForce = Fn*scg->normal; // normal Force (vector) + + + /*** SHEAR FORCE ***/ + phys->ks = phys->kso*std::pow(Fn,1/3); // normal stiffness + Vector3r& trialFs = phys->shearForce; + scg->updateShearForce(trialFs, phys->ks, phys->prevNormal, de1, de2, dt, preventGranularRatcheting); + + + /*** MOHR-COULOMB LAW ***/ + Real maxFs = phys->normalForce.SquaredLength(); + if (trialFs.SquaredLength() > maxFs) + {Real ratio = Mathr::Sqrt(maxFs)/trialFs.Length(); trialFs *= ratio;} + phys->shearForce = trialFs; // store shearForce now that is at the actual value + + + /*** APPLY FORCES ***/ + applyForceAtContactPoint(-phys->normalForce-trialFs , scg->contactPoint , id1, de1->se3.position, id2, de2->se3.position, ncb); + phys->prevNormal = scg->normal; +} + + + + + + + + + + + + + + === added file 'pkg/dem/Engine/GlobalEngine/HertzMindlin.hpp' --- pkg/dem/Engine/GlobalEngine/HertzMindlin.hpp 1970-01-01 00:00:00 +0000 +++ pkg/dem/Engine/GlobalEngine/HertzMindlin.hpp 2010-02-23 14:22:35 +0000 @@ -0,0 +1,66 @@ +// 2010 © Chiara Modenese <[email protected]> +// +/* +=== HIGH LEVEL OVERVIEW OF Mindlin === + +Mindlin is a set of classes to include the Hertz-Mindlin formulation for the contact stiffnesses. + +*/ + +#pragma once + +#include<yade/pkg-common/ElastMat.hpp> +#include<yade/pkg-dem/ScGeom.hpp> +#include<yade/pkg-common/InteractionPhysicsFunctor.hpp> +#include<yade/pkg-dem/FrictPhys.hpp> +#include<yade/pkg-common/PeriodicEngines.hpp> +#include<yade/pkg-common/NormShearPhys.hpp> +#include<yade/pkg-common/LawFunctor.hpp> + + +/******************** MindlinPhys *********************************/ +class MindlinPhys: public FrictPhys{ + public: + virtual ~MindlinPhys(); + YADE_CLASS_BASE_DOC_ATTRS_CTOR(MindlinPhys,FrictPhys,"Representation of an interaction of the Mindlin type.", + ((Real,kno,0.0,"Constant value in the formulation of the normal stiffness")) + ((Real,kso,0.0,"Constant value in the formulation of the tangential stiffness")), + createIndex()); + REGISTER_CLASS_INDEX(MindlinPhys,FrictPhys); + DECLARE_LOGGER; +}; +REGISTER_SERIALIZABLE(MindlinPhys); + + +/******************** Ip2_FrictMat_FrictMat_MindlinPhys *******/ +class Ip2_FrictMat_FrictMat_MindlinPhys: public InteractionPhysicsFunctor{ + public : + virtual void go(const shared_ptr<Material>& b1, const shared_ptr<Material>& b2, const shared_ptr<Interaction>& interaction); + FUNCTOR2D(FrictMat,FrictMat); + YADE_CLASS_BASE_DOC(Ip2_FrictMat_FrictMat_MindlinPhys,InteractionPhysicsFunctor,"Calculate some physical parameters needed to obtain the normal and shear stiffnesses according to the Hertz-Mindlin's formulation (as implemented in PFC)."); + DECLARE_LOGGER; +}; +REGISTER_SERIALIZABLE(Ip2_FrictMat_FrictMat_MindlinPhys); + + +/******************** Law2_ScGeom_MindlinPhys_Mindlin *********/ +class Law2_ScGeom_MindlinPhys_Mindlin: public LawFunctor{ + public: + virtual void go(shared_ptr<InteractionGeometry>& _geom, shared_ptr<InteractionPhysics>& _phys, Interaction* I, Scene* rootBody); + FUNCTOR2D(ScGeom,MindlinPhys); + YADE_CLASS_BASE_DOC_ATTRS(Law2_ScGeom_MindlinPhys_Mindlin,LawFunctor,"Constitutive law for the Mindlin's formulation.", + ((bool,preventGranularRatcheting,false,"bool to avoid granular ratcheting")) + + ); + DECLARE_LOGGER; +}; +REGISTER_SERIALIZABLE(Law2_ScGeom_MindlinPhys_Mindlin); + + + + + + + + + === added file 'scripts/mindlin_test.py' --- scripts/mindlin_test.py 1970-01-01 00:00:00 +0000 +++ scripts/mindlin_test.py 2010-02-23 14:22:35 +0000 @@ -0,0 +1,73 @@ +#!/usr/local/bin/yade-trunk -x +# -*- coding: utf-8 -*- +# -*- encoding=utf-8 -*- +## +## SCRIPT TO TEST A NEW CONSTITUTIVE LAW (MINDLIN - nonlinear elastic model) + +O.initializers=[BoundDispatcher([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()])] + +## list of engines +O.engines=[ + ForceResetter(), + BoundDispatcher([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]), + InsertionSortCollider(), + InteractionDispatchers( + [Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()], + [Ip2_FrictMat_FrictMat_MindlinPhys()], + [Law2_ScGeom_MindlinPhys_Mindlin()] + ), + GravityEngine(gravity=(-100,0,0)), + NewtonIntegrator(damping=0.0), + ### + ### NOTE this extra engine: + ### + ### You want snapshot to be taken every 1 sec (realTimeLim) or every 50 iterations (iterLim), + ### whichever comes soones. virtTimeLim attribute is unset, hence virtual time period is not taken into account. + PeriodicPythonRunner(iterPeriod=1,command='myAddPlotData()') +] + +## define and append material +mat=FrictMat(young=600.0e6,poisson=0.6,density=2.60e3,frictionAngle=26,label='Friction') +O.materials.append(mat) + +## create two spheres (one will be fixed) and append them +from yade import utils +s0=utils.sphere([0,0,0],1,color=[0,1,0],dynamic=False,wire=True,material='Friction') +s1=utils.sphere([2,0,0],1,color=[0,2,0],dynamic=True,wire=True,material='Friction') +O.bodies.append(s0) +O.bodies.append(s1) + +## time step +O.dt=.2*utils.PWaveTimeStep() +O.saveTmp('Mindlin') + +from yade import qt +qt.View() +qt.Controller() + +############################################ +##### now the part pertaining to plots ##### +############################################ + +from math import * +from yade import plot +## make one plot: step as function of fn +plot.plots={'un':('fn')} + +## this function is called by plotDataCollector +## it should add data with the labels that we will plot +## if a datum is not specified (but exists), it will be NaN and will not be plotted + +def myAddPlotData(): + i=O.interactions[0,1] + ## store some numbers under some labels + plot.addData(fn=i.phys.normalForce[0],step=O.iter,un=2*s0.shape['radius']-s1.state.pos[0]-s0.state.pos[0],kn=i.phys.kn) + +O.run(50,True); +print "Now calling plot.plot() to show the figure." + +## We will have: +## 1) data in graphs (if you call plot.plot()) +## 2) data in file (if you call plot.saveGnuplot('/tmp/a') +## 3) data in memory as plot.data['step'], plot.data['fn'], plot.data['un'], etc. under the labels they were saved +
_______________________________________________ Mailing list: https://launchpad.net/~yade-dev Post to : [email protected] Unsubscribe : https://launchpad.net/~yade-dev More help : https://help.launchpad.net/ListHelp
