yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #06959
[Branch ~yade-dev/yade/trunk] Rev 2709: - Add calculation of fabric tensor for periodic cell as according to Satake (bib is updated too)....
------------------------------------------------------------
revno: 2709
committer: Chiara Modenese <c.modenese@xxxxxxxxx>
branch nick: yade
timestamp: Fri 2011-02-04 18:00:01 +0000
message:
- Add calculation of fabric tensor for periodic cell as according to Satake (bib is updated too). The tensor can also be spit into two different parts each of them related to strong and weak contact forces respectively.
- Small changes in HM.
modified:
doc/references.bib
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 'doc/references.bib'
--- doc/references.bib 2011-01-21 19:55:58 +0000
+++ doc/references.bib 2011-02-04 18:00:01 +0000
@@ -6,7 +6,7 @@
1. Keep entries in the form Author2008 (Author is the first author), Author2008b if repeated
-2. Try to fill mandatory fields for given type of citations (http://en.wikipedia.org/wiki/Bibtex\\\#Entry\\\_Types)
+2. Try to fill mandatory fields for given type of citations (http://en.wikipedia.org/wiki/Bibtex\\\\\#Entry\\\\\_Types)
3. Do not use {\'i} funny escapes for accents, put everything in straight utf-8
@@ -24,7 +24,7 @@
year = "2008"
}
-@article{ Bertrand2005,
+@Article{ Bertrand2005,
title = "Modelling a geo-composite cell using discrete analysis",
author = "D. Bertrand and F. Nicot and P. Gotteland and S. Lambert",
journal = "Computers and Geotechnics",
@@ -57,11 +57,11 @@
doi = "10.1061/(ASCE)0733-9399(2005)131:7(689)"
}
-@Article{Duriez2010,
+@Article{ Duriez2010,
author = "J. Duriez and F.Darve and F.-V.Donze",
title = "A discrete modeling-based constitutive relation for infilled rock joints",
year = "2010",
- journal = "International Journal of Rock Mechanics & Mining Sciences",
+ journal = "International Journal of Rock Mechanics \& Mining Sciences",
note = "in press"
}
@@ -403,3 +403,12 @@
Pages = "5--19",
Year = "2002"
}
+
+@InProceedings{ Satake1982,
+ title = "Fabric tensor in granular materials.",
+ author = "M. Satake",
+ booktitle = "Proc., IUTAM Symp. on Deformation and Failure of Granular materials, Delft, The Netherlands",
+ year = "1982",
+ pages = "63--68"
+}
+
=== modified file 'pkg/dem/HertzMindlin.cpp'
--- pkg/dem/HertzMindlin.cpp 2011-01-29 22:47:18 +0000
+++ pkg/dem/HertzMindlin.cpp 2011-02-04 18:00:01 +0000
@@ -239,7 +239,8 @@
const shared_ptr<Body>& b2=Body::byId(id2,scene);
bool useDamping=(phys->betan!=0. || phys->betas!=0. || phys->alpha!=0.);
- if(phys->alpha==0. && !LinDamp) throw std::invalid_argument("Law2_ScGeom_MindlinPhys_Mindlin: viscous contact damping for non-linear elastic law is specified, specify en and es in Ip2.");
+ bool LinDamp=true;
+ if (phys->alpha!=0.) {LinDamp=false;} // use non linear damping
// tangential and normal stiffness coefficients, recomputed from betan,betas at every step
Real cn=0, cs=0;
=== modified file 'pkg/dem/HertzMindlin.hpp'
--- pkg/dem/HertzMindlin.hpp 2011-01-26 18:19:52 +0000
+++ pkg/dem/HertzMindlin.hpp 2011-02-04 18:00:01 +0000
@@ -136,7 +136,7 @@
((bool,includeAdhesion,false,,"bool to include the adhesion force following the DMT formulation. If true, also the normal elastic energy takes into account the adhesion effect."))
((bool,calcEnergy,false,,"bool to calculate energy terms (shear potential energy, dissipation of energy due to friction and dissipation of energy due to normal and tangential damping)"))
((bool,includeMoment,false,,"bool to consider rolling resistance (if :yref:`Ip2_FrictMat_FrictMat_MindlinPhys::eta` is 0.0, no plastic condition is applied.)"))
- ((bool,LinDamp,true,,"bool to activate linear viscous damping (if false, en and es have to be defined in place of betan and betas)"))
+ //((bool,LinDamp,true,,"bool to activate linear viscous damping (if false, en and es have to be defined in place of betan and betas)"))
((OpenMPAccumulator<Real>,frictionDissipation,,Attr::noSave,"Energy dissipation due to sliding"))
((OpenMPAccumulator<Real>,shearEnergy,,Attr::noSave,"Shear elastic potential energy"))
=== modified file 'pkg/dem/Shop.cpp'
--- pkg/dem/Shop.cpp 2011-01-29 22:47:18 +0000
+++ pkg/dem/Shop.cpp 2011-02-04 18:00:01 +0000
@@ -527,6 +527,70 @@
return py::make_tuple(sigN,sigT);
}
+/* Return the fabric tensor as according to [Satake1982]. */
+py::tuple Shop::fabricTensor(bool splitTensor){
+ Scene* scene=Omega::instance().getScene().get();
+ if (!scene->isPeriodic){ throw runtime_error("Can't compute fabric tensor of periodic cell in aperiodic simulation."); }
+ Matrix3r fabric(Matrix3r::Zero());
+ int count=0;
+ FOREACH(const shared_ptr<Interaction>& I, *scene->interactions){
+ if(!I->isReal()) continue;
+ GenericSpheresContact* geom=YADE_CAST<GenericSpheresContact*>(I->geom.get());
+ const Vector3r& n=geom->normal;
+ for(int i=0; i<3; i++) for(int j=i; j<3; j++){
+ fabric(i,j)+=n[i]*n[j];
+ }
+ count++;
+ }
+ // fill terms under the diagonal
+ fabric(1,0)=fabric(0,1); fabric(2,0)=fabric(0,2); fabric(2,1)=fabric(1,2);
+ fabric/=count;
+
+ // calculate average contact force
+ Real Fmean(0);
+ FOREACH(const shared_ptr<Interaction>& I, *scene->interactions){
+ if(!I->isReal()) continue;
+ GenericSpheresContact* geom=YADE_CAST<GenericSpheresContact*>(I->geom.get());
+ NormShearPhys* phys=YADE_CAST<NormShearPhys*>(I->phys.get());
+ Vector3r Fcontact=phys->shearForce+phys->normalForce;
+ Real f=Fcontact.norm();
+ Fmean+=f;
+ }
+ Fmean/=count;
+
+ // evaluate two different parts of the fabric tensor
+ // make distinction between strong and weak network of contact forces
+ Matrix3r fabricStrong(Matrix3r::Zero()), fabricWeak(Matrix3r::Zero());
+ int nStrong(0), nWeak(0);
+ FOREACH(const shared_ptr<Interaction>& I, *scene->interactions){
+ if(!I->isReal()) continue;
+ GenericSpheresContact* geom=YADE_CAST<GenericSpheresContact*>(I->geom.get());
+ NormShearPhys* phys=YADE_CAST<NormShearPhys*>(I->phys.get());
+ Vector3r Fcontact=phys->shearForce+phys->normalForce;
+ Real f=Fcontact.norm();
+ const Vector3r& n=geom->normal;
+ if (f>Fmean){
+ for(int i=0; i<3; i++) for(int j=i; j<3; j++){
+ fabricStrong(i,j)+=n[i]*n[j];
+ }
+ nStrong++;
+ }
+ else{
+ for(int i=0; i<3; i++) for(int j=i; j<3; j++){
+ fabricWeak(i,j)+=n[i]*n[j];
+ }
+ nWeak++;
+ }
+ }
+ // fill terms under the diagonal
+ fabricStrong(1,0)=fabricStrong(0,1); fabricStrong(2,0)=fabricStrong(0,2); fabricStrong(2,1)=fabricStrong(1,2);
+ fabricWeak(1,0)=fabricWeak(0,1); fabricWeak(2,0)=fabricWeak(0,2); fabricWeak(2,1)=fabricWeak(1,2);
+ fabricStrong/=nStrong;
+ fabricWeak/=nWeak;
+
+ if (!splitTensor){return py::make_tuple(fabric);}
+ else{return py::make_tuple(fabricStrong,fabricWeak);}
+}
Matrix3r Shop::stressTensorOfPeriodicCell(bool smallStrains){
Scene* scene=Omega::instance().getScene().get();
=== modified file 'pkg/dem/Shop.hpp'
--- pkg/dem/Shop.hpp 2011-01-17 15:44:53 +0000
+++ pkg/dem/Shop.hpp 2011-02-04 18:00:01 +0000
@@ -113,4 +113,7 @@
//! Compute overall ("macroscopic") stress of periodic cell, returning 2 tensors
//! (contribution of normal and shear forces)
static py::tuple normalShearStressTensors(bool compressionPositive=false);
+
+ //! Function to compute fabric tensor of periodic cell
+ static py::tuple fabricTensor(bool splitTensor=false);
};
=== modified file 'py/_utils.cpp'
--- py/_utils.cpp 2011-01-29 22:47:18 +0000
+++ py/_utils.cpp 2011-02-04 18:00:01 +0000
@@ -429,6 +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){return Shop::fabricTensor(splitTensor);}
py::tuple Shop__normalShearStressTensors(bool compressionPositive=false){return Shop::normalShearStressTensors(compressionPositive);}
py::tuple Shop__getStressLWForEachBody(bool revertSign=false){return Shop::getStressLWForEachBody(revertSign);}
@@ -475,6 +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),"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.");
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.");
}
