yade-dev team mailing list archive

[noreply@xxxxxxxxxxxxx]

------------------------------------------------------------
revno: 2049
committer: Chiara Modenese <chia@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 <c.modenese@xxxxxxxxx> 
+
+
+#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 <c.modenese@xxxxxxxxx> 
+// 
+/*
+=== 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
+