yade-dev team mailing list archive

Thread
Date

[Branch ~yade-dev/yade/trunk] Rev 2717: - Add a function to HMLaw which computes the ratio of sliding contacts.

To: Yade developers <yade-dev@xxxxxxxxxxxxxxxxxxx>
From: noreply@xxxxxxxxxxxxx
Date: Tue, 08 Feb 2011 16:37:28 -0000
Reply-to: noreply@xxxxxxxxxxxxx
Sender: bounces@xxxxxxxxxxxxx

------------------------------------------------------------
revno: 2717
committer: Chiara Modenese <c.modenese@xxxxxxxxx>
branch nick: yade
timestamp: Tue 2011-02-08 16:37:12 +0000
message:
  - Add a function to HMLaw which computes the ratio of sliding contacts. 
  - Adjust the sign convention in the computation of the fabric tensor. Now the tensor can also be split into tensile and compressive part (simply setting the threshold value equal to 0 - added documentation about that).
  - Rename (for the fabric tensor function only) compressionPositive to revertSign (that is the meaning indeed - work in progress to use the scene flag and make automatic distinction whichever law is used).
modified:
  pkg/dem/HertzMindlin.cpp
  pkg/dem/HertzMindlin.hpp
  pkg/dem/Shop.cpp
  pkg/dem/Shop.hpp
  py/_utils.cpp


--
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 'pkg/dem/HertzMindlin.cpp'
--- pkg/dem/HertzMindlin.cpp	2011-02-04 18:00:01 +0000
+++ pkg/dem/HertzMindlin.cpp	2011-02-08 16:37:12 +0000
@@ -108,6 +108,20 @@
 	return contactsAdhesive;
 }
 
+/* Function which returns the ratio between the number of sliding contacts to the total number at a given time */
+Real Law2_ScGeom_MindlinPhys_Mindlin::ratioSlidingContacts()
+{
+	Real ratio(0); int count(0);
+	FOREACH(const shared_ptr<Interaction>& I, *scene->interactions){
+		if(!I->isReal()) continue;
+		MindlinPhys* phys = dynamic_cast<MindlinPhys*>(I->phys.get());
+		if (phys->isSliding) {ratio+=1;}
+		count++;
+	}  
+	ratio/=count;
+	return ratio;
+}
+
 /* Function to get the NORMAL elastic potential energy of the system */
 Real Law2_ScGeom_MindlinPhys_Mindlin::normElastEnergy()
 {
@@ -358,6 +372,7 @@
 	if (!includeAdhesion) {
 		Real maxFs = Fn*phys->tangensOfFrictionAngle;
 		if (shearElastic.squaredNorm() > maxFs*maxFs){
+			phys->isSliding=true;
 			noShearDamp = true; // no damping is added in the shear direction, hence no need to account for shear damping dissipation
 			Real ratio = maxFs/shearElastic.norm();
 			shearElastic *= ratio; phys->shearForce = shearElastic; /*store only elastic shear displacement*/ us_elastic*= ratio;
@@ -368,9 +383,10 @@
 			phys->shearForce = shearElastic - phys->shearViscous;}
 		else if (!useDamping) {phys->shearForce = shearElastic;} // update the shear force at the elastic value if no damping is present and if we passed MC
 	}
-	else { // Mohr-Coulomb formulation adpated due to the presence of adhesion (see Thornton, 1991). FIXME: is this correct?
+	else { // Mohr-Coulomb formulation adpated due to the presence of adhesion (see Thornton, 1991).
 		Real maxFs = phys->tangensOfFrictionAngle*(phys->adhesionForce+Fn); // adhesionForce already included in normalForce (above)
 		if (shearElastic.squaredNorm() > maxFs*maxFs){
+			phys->isSliding=true;
 			noShearDamp = true; // no damping is added in the shear direction, hence no need to account for shear damping dissipation
 			Real ratio = maxFs/shearElastic.norm(); shearElastic *= ratio; phys->shearForce = shearElastic; /*store only elastic shear displacement*/ us_elastic *= ratio;
 			if (calcEnergy) {frictionDissipation += (us_total-prevUs_tot).dot(shearElastic);} // calculate energy dissipation due to sliding behavior

=== modified file 'pkg/dem/HertzMindlin.hpp'
--- pkg/dem/HertzMindlin.hpp	2011-02-04 18:00:01 +0000
+++ pkg/dem/HertzMindlin.hpp	2011-02-08 16:37:12 +0000
@@ -51,6 +51,7 @@
 			//((Real,gamma,0.0,"Surface energy parameter [J/m^2] per each unit contact surface, to derive DMT formulation from HM"))
 			((Real,adhesionForce,0.0,,"Force of adhesion as predicted by DMT"))
 			((bool,isAdhesive,false,,"bool to identify if the contact is adhesive, that is to say if the contact force is attractive"))
+			((bool,isSliding,false,,"check if the contact is sliding (useful to calculate the ratio of sliding contacts)"))
 
 			// Contact damping coefficients as for linear elastic contact law
 			((Real,betan,0.0,,"Fraction of the viscous damping coefficient (normal direction) equal to $\\frac{c_{n}}{C_{n,crit}}$."))
@@ -129,6 +130,7 @@
 		Real getnormDampDissip();
 		Real getshearDampDissip();
 		Real contactsAdhesive();
+		Real ratioSlidingContacts();
 
 		FUNCTOR2D(ScGeom,MindlinPhys);
 		YADE_CLASS_BASE_DOC_ATTRS_DEPREC_INIT_CTOR_PY(Law2_ScGeom_MindlinPhys_Mindlin,LawFunctor,"Constitutive law for the Hertz-Mindlin formulation. It includes non linear elasticity in the normal direction as predicted by Hertz for two non-conforming elastic contact bodies. In the shear direction, instead, it reseambles the simplified case without slip discussed in Mindlin's paper, where a linear relationship between shear force and tangential displacement is provided. Finally, the Mohr-Coulomb criterion is employed to established the maximum friction force which can be developed at the contact. Moreover, it is also possible to include the effect of linear viscous damping through the definition of the parameters $\\beta_{n}$ and $\\beta_{s}$.",
@@ -151,6 +153,7 @@
 			, /* ctor */
 			, /* py */
 			.def("contactsAdhesive",&Law2_ScGeom_MindlinPhys_Mindlin::contactsAdhesive,"Compute total number of adhesive contacts.")
+			.def("ratioSlidingContacts",&Law2_ScGeom_MindlinPhys_Mindlin::ratioSlidingContacts,"Return the ratio between the number of contacts sliding to the total number at a given time.")
 			.def("normElastEnergy",&Law2_ScGeom_MindlinPhys_Mindlin::normElastEnergy,"Compute normal elastic potential energy. It handles the DMT formulation if :yref:`Law2_ScGeom_MindlinPhys_Mindlin::includeAdhesion` is set to true.")	
 	);
 	DECLARE_LOGGER;

=== modified file 'pkg/dem/Shop.cpp'
--- pkg/dem/Shop.cpp	2011-02-08 10:17:07 +0000
+++ pkg/dem/Shop.cpp	2011-02-08 16:37:12 +0000
@@ -533,7 +533,7 @@
 
 /* Return the fabric tensor as according to [Satake1982]. */
 /* as side-effect, set Gl1_NormShear::strongWeakThresholdForce */
-py::tuple Shop::fabricTensor(bool splitTensor, bool compressionPositive, Real thresholdForce){
+py::tuple Shop::fabricTensor(bool splitTensor, bool revertSign, Real thresholdForce){
 	Scene* scene=Omega::instance().getScene().get();
 	if (!scene->isPeriodic){ throw runtime_error("Can't compute fabric tensor of periodic cell in aperiodic simulation."); }
 	
@@ -561,7 +561,7 @@
 		GenericSpheresContact* geom=YADE_CAST<GenericSpheresContact*>(I->geom.get());
 		NormShearPhys* phys=YADE_CAST<NormShearPhys*>(I->phys.get());
 		const Vector3r& n=geom->normal;
-		Real f=(compressionPositive?-1:1)*phys->normalForce.dot(n);  // FIXME: is it right to consider the sign here?
+		Real f=(revertSign?-1:1)*phys->normalForce.dot(n);  
 		//Real f=phys->normalForce.norm();
 		Fmean+=f;
 	}
@@ -582,13 +582,10 @@
 		GenericSpheresContact* geom=YADE_CAST<GenericSpheresContact*>(I->geom.get());
 		NormShearPhys* phys=YADE_CAST<NormShearPhys*>(I->phys.get());
 		const Vector3r& n=geom->normal;
-		Real  f=(compressionPositive?-1:1)*phys->normalForce.dot(n); 
-		//Real f=phys->normalForce.norm();
-		// FIXME: is it right to check only the normal force against Fmean?
-		// what if we have adhesive contacts?
+		Real  f=(revertSign?-1:1)*phys->normalForce.dot(n); 
 		// slipt the tensor according to the mean contact force or a threshold value if this is given
 		Real Fsplit=(!isnan(thresholdForce))?thresholdForce:Fmean;
-		if (f>Fsplit){ // forces are compared with their sign
+		if (revertSign?(f<Fsplit):(f>Fsplit)){ // reminder: forces are compared with their sign
 			for(int i=0; i<3; i++) for(int j=i; j<3; j++){
 				fabricStrong(i,j)+=n[i]*n[j];
 			}
@@ -615,8 +612,8 @@
 		fabricTot=(1-q)*fabricWeak+q*fabricStrong;
 	}
 	
-	// return fabric tensor or alternatively the two distinct contributions according to strong and weak subnetworks
-	if (!splitTensor){return py::make_tuple(fabric);} // return fabricTot for the 
+	// returns fabric tensor or alternatively the two distinct contributions according to strong and weak subnetworks (or, if thresholdForce is specified, the distinction is made according to that value and not the mean one)
+	if (!splitTensor){return py::make_tuple(fabric);} 
 	else{return py::make_tuple(fabricStrong,fabricWeak);}
 }
 

=== modified file 'pkg/dem/Shop.hpp'
--- pkg/dem/Shop.hpp	2011-02-08 10:17:07 +0000
+++ pkg/dem/Shop.hpp	2011-02-08 16:37:12 +0000
@@ -115,5 +115,5 @@
 		static py::tuple normalShearStressTensors(bool compressionPositive=false);
 		
 		//! Function to compute fabric tensor of periodic cell
-		static py::tuple fabricTensor(bool splitTensor=false, bool compressionPositive=true, Real thresholdForce=NaN);
+		static py::tuple fabricTensor(bool splitTensor=false, bool revertSign=false, Real thresholdForce=NaN);
 };

=== modified file 'py/_utils.cpp'
--- py/_utils.cpp	2011-02-08 10:17:07 +0000
+++ py/_utils.cpp	2011-02-08 16:37:12 +0000
@@ -429,7 +429,7 @@
 Real Shop__getPorosity(Real volume=-1){ return Shop::getPorosity(Omega::instance().getScene(),volume); }
 
 Matrix3r Shop__stressTensorOfPeriodicCell(bool smallStrains=false){return Shop::stressTensorOfPeriodicCell(smallStrains);}
-py::tuple Shop__fabricTensor(bool splitTensor=false, bool compressionPositive=true, Real thresholdForce=NaN){return Shop::fabricTensor(splitTensor,compressionPositive,thresholdForce);}
+py::tuple Shop__fabricTensor(bool splitTensor=false, bool revertSign=false, Real thresholdForce=NaN){return Shop::fabricTensor(splitTensor,revertSign,thresholdForce);}
 py::tuple Shop__normalShearStressTensors(bool compressionPositive=false){return Shop::normalShearStressTensors(compressionPositive);}
 py::tuple Shop__getStressLWForEachBody(bool revertSign=false){return Shop::getStressLWForEachBody(revertSign);}
 
@@ -476,7 +476,7 @@
 	py::def("getViscoelasticFromSpheresInteraction",getViscoelasticFromSpheresInteraction,(py::arg("tc"),py::arg("en"),py::arg("es")),"Compute viscoelastic interaction parameters from analytical solution of a pair spheres collision problem:\n\n\n.. math::\n\t:nowrap:\n\n\n\t\\begin{align*}k_n&=\\frac{m}{t_c^2}\\left(\\pi^2+(\\ln e_n)^2\\right)\\\\  c_n&=-\\frac{2m}{t_c}\\ln e_n \\\\k_t&=\\frac27\\frac{m}{t_c^2}\\left(\\pi^2+(\\ln e_t)^2\\right)  \\\\ c_t=-\\frac27\\frac{m}{t_c}\\ln e_t \\end{align*}\n\n\nwhere $k_n$, $c_n$ are normal elastic and viscous coefficients and $k_t$, $c_t$ shear elastic and viscous coefficients. For details see [Pournin2001]_.\n\n:param float m: sphere mass $m$\n:param float tc: collision time $t_c$\n:param float en: normal restitution coefficient $e_n$\n:param float es: tangential restitution coefficient $e_s$\n:return: dictionary with keys ``kn`` (the value of $k_n$), ``cn`` ($c_n$), ``kt`` ($k_t$), ``ct`` ($c_t$).");
 	py::def("stressTensorOfPeriodicCell",Shop__stressTensorOfPeriodicCell,(py::args("smallStrains")=false),"Compute overall (macroscopic) stress of periodic cell using equation published in [Kuhl2001]_:\n\n.. math:: \\vec{\\sigma}=\\frac{1}{V}\\sum_cl^c[\\vec{N}^cf_N^c+\\vec{T}^{cT}\\cdot\\vec{f}^c_T],\n\nwhere $V$ is volume of the cell, $l^c$ length of interaction $c$, $f^c_N$ normal force and $\\vec{f}^c_T$ shear force. Sumed are values over all interactions $c$. $\\vec{N}^c$ and $\\vec{T}^{cT}$ are projection tensors (see the original publication for more details):\n\n.. math:: \\vec{N}=\\vec{n}\\otimes\\vec{n}\\rightarrow N_{ij}=n_in_j\n\n.. math:: \\vec{T}^T=\\vec{I}_{sym}\\cdot\\vec{n}-\\vec{n}\\otimes\\vec{n}\\otimes\\vec{n}\\rightarrow T^T_{ijk}=\\frac{1}{2}(\\delta_{ik}\\delta_{jl}+\\delta_{il}\\delta_{jk})n_l-n_in_jn_k\n\n.. math:: \\vec{T}^T\\cdot\\vec{f}_T\\equiv T^T_{ijk}f_k=(\\delta_{ik}n_j/2+\\delta_{jk}n_i/2-n_in_jn_k)f_k=n_jf_i/2+n_if_j/2-n_in_jn_kf_k,\n\nwhere $n$ is unit vector oriented along the interaction (:yref:`normal<GenericSpheresContact::normal>`) and $\\delta$ is Kronecker's delta. As $\\vec{n}$ and $\\vec{f}_T$ are perpendicular (therfore $n_if_i=0$) we can write\n\n.. math:: \\sigma_{ij}=\\frac{1}{V}\\sum l[n_in_jf_N+n_jf^T_i/2+n_if^T_j/2]\n\n:param bool smallStrains: if false (large strains), real values of volume and interaction lengths are computed. If true, only :yref:`refLength<Dem3DofGeom::refLength>` of interactions and initial volume are computed (can save some time).\n\n:return: macroscopic stress tensor as Matrix3");
 	py::def("normalShearStressTensors",Shop__normalShearStressTensors,(py::args("compressionPositive")=false),"Compute overall stress tensor of the periodic cell decomposed in 2 parts, one contributed by normal forces, the other by shear forces. The formulation can be found in [Thornton2000]_, eq. (3):\n\n.. math:: \\tens{\\sigma}_{ij}=\\frac{2}{V}\\sum R N \\vec{n}_i \\vec{n}_j+\\frac{2}{V}\\sum R T \\vec{n}_i\\vec{t}_j\n\nwhere $V$ is the cell volume, $R$ is \"contact radius\" (in our implementation, current distance between particle centroids), $\\vec{n}$ is the normal vector, $\\vec{t}$ is a vector perpendicular to $\\vec{n}$, $N$ and $T$ are norms of normal and shear forces.");
-	py::def("fabricTensor",Shop__fabricTensor,(py::args("splitTensor")=false,py::args("compressionPositive")=true,py::args("thresholdForce")=NaN),"Compute the fabric tensor of the periodic cell. The original paper can be found in [Satake1982]_.\n\n:param bool splitTensor: split the fabric tensor into two parts related to the strong and weak contact forces respectively.\n\n:param bool compressionPositive: compressive forces are considered as positive.\n\n:param Real thresholdForce: if the fabric tensor is split into two parts, a threshold value can be specified otherwise the mean contact force is considered by default.");
+	py::def("fabricTensor",Shop__fabricTensor,(py::args("splitTensor")=false,py::args("revertSign")=false,py::args("thresholdForce")=NaN),"Compute the fabric tensor of the periodic cell. The original paper can be found in [Satake1982]_.\n\n:param bool splitTensor: split the fabric tensor into two parts related to the strong and weak contact forces respectively.\n\n:param bool revertSign: it must be set to true if the contact law's convention takes compressive forces as positive.\n\n:param Real thresholdForce: if the fabric tensor is split into two parts, a threshold value can be specified otherwise the mean contact force is considered by default. It is worth to note that this value has a sign and the user needs to set it according to the convention adopted for the contact law. To note that this value could be set to zero if one wanted to make distinction between compressive and tensile forces.");
 	py::def("bodyStressTensors",Shop__getStressLWForEachBody,(py::args("revertSign")=false),"Compute and return a table with per-particle stress tensors. Each tensor represents the average stress in one particle, obtained from the contour integral of applied load as detailed below. This definition is considering each sphere as a continuum. It can be considered exact in the context of spheres at static equilibrium, interacting at contact points with negligible volume changes of the solid phase (this last assumption is not restricting possible deformations and volume changes at the packing scale).\n\nProof:\n\nFirst, we remark the identity:  $\\sigma_{ij}=\\delta_{ij}\\sigma_{ij}=x_{i,j}\\sigma_{ij}=(x_{i}\\sigma_{ij})_{,j}-x_{i}\\sigma_{ij,j}$.\n\nAt equilibrium, the divergence of stress is null: $\\sigma_{ij,j}=\\vec{0}$. Consequently, after divergence theorem: $\\frac{1}{V}\\int_V \\sigma_{ij}dV = \\frac{1}{V}\\int_V (x_{i}\\sigma_{ij})_{,j}dV = \\frac{1}{V}\\int_{\\partial V}x_i.\\sigma_{ij}.\\vec{n_j}.dS = \\frac{1}{V}\\sum_kx_i^k.f_j^k$.\n\nThe last equality is implicitely based on the representation of external loads as Dirac distributions whose zeros are the so-called *contact points*: 0-sized surfaces on which the *contact forces* are applied, located at $x_i$ in the deformed configuration.\n\nA weighted average of per-body stresses will give the average stress inside the solid phase. There is a simple relation between the stress inside the solid phase and the stress in an equivalent continuum in the absence of fluid pressure. For porosity $n$, the relation reads: $\\sigma_{ij}^{equ.}=(1-n)\\sigma_{ij}^{solid}$.\n\n:param bool revertSign: invert the sign of returned tensors components.");
 	py::def("maxOverlapRatio",maxOverlapRatio,"Return maximum overlap ration in interactions (with :yref:`ScGeom`) of two :yref:`spheres<Sphere>`. The ratio is computed as $\\frac{u_N}{2(r_1 r_2)/r_1+r_2}$, where $u_N$ is the current overlap distance and $r_1$, $r_2$ are radii of the two spheres in contact.");
 }