yade-dev team mailing list archive

[noreply@xxxxxxxxxxxxx]

------------------------------------------------------------
revno: 2508
committer: Sergei D. <sega@think>
branch nick: trunk
timestamp: Fri 2010-10-22 18:30:53 +0400
message:
  Viscoelastic basic particle model: move mass of particles from material to interaction
modified:
  pkg/dem/Shop.cpp
  pkg/dem/Shop.hpp
  pkg/dem/ViscoelasticPM.cpp
  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/Shop.cpp'
--- pkg/dem/Shop.cpp	2010-10-19 14:57:27 +0000
+++ pkg/dem/Shop.cpp	2010-10-22 14:30:53 +0000
@@ -1158,12 +1158,12 @@
 	return v0+((v2-v0)*(v1-v0).norm()+(v1-v0)*(v2-v0).norm())/((v1-v0).norm()+(v2-v1).norm()+(v0-v2).norm());
 }
 
-void Shop::getViscoelasticFromSpheresInteraction( Real m, Real tc, Real en, Real es, shared_ptr<ViscElMat> b)
+void Shop::getViscoelasticFromSpheresInteraction( Real tc, Real en, Real es, shared_ptr<ViscElMat> b)
 {
-    b->kn = m/tc/tc * ( Mathr::PI*Mathr::PI + pow(log(en),2) );
-    b->cn = -2.0*m/tc * log(en);
-    b->ks = 2.0/7.0 * m/tc/tc * ( Mathr::PI*Mathr::PI + pow(log(es),2) );
-    b->cs = -2.0/7.0 * m/tc * log(es);
+    b->kn = 1/tc/tc * ( Mathr::PI*Mathr::PI + pow(log(en),2) );
+    b->cn = -2.0 /tc * log(en);
+    b->ks = 2.0/7.0 /tc/tc * ( Mathr::PI*Mathr::PI + pow(log(es),2) );
+    b->cs = -2.0/7.0 /tc * log(es);
 
     if (abs(b->cn) <= Mathr::ZERO_TOLERANCE ) b->cn=0;
     if (abs(b->cs) <= Mathr::ZERO_TOLERANCE ) b->cs=0;

=== modified file 'pkg/dem/Shop.hpp'
--- pkg/dem/Shop.hpp	2010-10-13 16:23:08 +0000
+++ pkg/dem/Shop.hpp	2010-10-22 14:30:53 +0000
@@ -94,9 +94,9 @@
 		static Vector3r inscribedCircleCenter(const Vector3r& v0, const Vector3r& v1, const Vector3r& v2);
 
 		/// Get viscoelastic parameters kn,cn,ks,cs from analytical solution of
-		/// a problem of interaction of pair spheres with mass m, collision
+		/// a problem of interaction of pair spheres with mass 1, collision
 		/// time tc and restitution coefficients en,es.
-	    static void getViscoelasticFromSpheresInteraction(Real m, Real tc, Real en, Real es, shared_ptr<ViscElMat> b);
+	    static void getViscoelasticFromSpheresInteraction(Real tc, Real en, Real es, shared_ptr<ViscElMat> b);
 
 		//! Get unbalanced force of the whole simulation
 		static Real unbalancedForce(bool useMaxForce=false, Scene* _rb=NULL);

=== modified file 'pkg/dem/ViscoelasticPM.cpp'
--- pkg/dem/ViscoelasticPM.cpp	2010-10-13 16:23:08 +0000
+++ pkg/dem/ViscoelasticPM.cpp	2010-10-22 14:30:53 +0000
@@ -15,17 +15,22 @@
 
 /* Ip2_ViscElMat_ViscElMat_ViscElPhys */
 void Ip2_ViscElMat_ViscElMat_ViscElPhys::go(const shared_ptr<Material>& b1, const shared_ptr<Material>& b2, const shared_ptr<Interaction>& interaction) {
+	// no updates of an existing contact 
 	if(interaction->phys) return;
-	ViscElMat* sdec1 = static_cast<ViscElMat*>(b1.get());
-	ViscElMat* sdec2 = static_cast<ViscElMat*>(b2.get());
+	ViscElMat* mat1 = static_cast<ViscElMat*>(b1.get());
+	ViscElMat* mat2 = static_cast<ViscElMat*>(b2.get());
+	const Real mass1 = Body::byId(interaction->getId1())->state->mass;
+	const Real mass2 = Body::byId(interaction->getId2())->state->mass;
+	const Real kn1 = mat1->kn*mass1; const Real cn1 = mat1->cn*mass1;
+	const Real ks1 = mat1->ks*mass1; const Real cs1 = mat1->cs*mass1;
+	const Real kn2 = mat2->kn*mass2; const Real cn2 = mat2->cn*mass2;
+	const Real ks2 = mat2->ks*mass2; const Real cs2 = mat2->cs*mass2;
 	ViscElPhys* phys = new ViscElPhys();
-	phys->kn = sdec1->kn * sdec2->kn / (sdec1->kn + sdec2->kn);
-	phys->ks = sdec1->ks * sdec2->ks / (sdec1->ks + sdec2->ks);
-	phys->cn = ( (sdec1->cn==0) ? 0 : 1/sdec1->cn ) + ( (sdec2->cn==0) ? 0 : 1/sdec2->cn );
-	phys->cs = ( (sdec1->cs==0) ? 0 : 1/sdec1->cs ) + ( (sdec2->cs==0) ? 0 : 1/sdec2->cs );
-	if (phys->cn) phys->cn = 1/phys->cn;
-	if (phys->cs) phys->cs = 1/phys->cs;
-	phys->tangensOfFrictionAngle = std::tan(std::min(sdec1->frictionAngle, sdec2->frictionAngle)); 
+	phys->kn = 1/(kn1?1/kn1:0 + kn2?1/kn2:0);
+	phys->ks = 1/(ks1?1/ks1:0 + ks2?1/ks2:0);
+	phys->cn = cn1?1/cn1:0 + cn2?1/cn2:0; phys->cn = phys->cn?1/phys->cn:0;
+	phys->cs = cs1?1/cs1:0 + cs2?1/cs2:0; phys->cs = phys->cs?1/phys->cs:0;
+	phys->tangensOfFrictionAngle = std::tan(std::min(mat1->frictionAngle, mat2->frictionAngle)); 
 	phys->shearForce = Vector3r(0,0,0);
 	phys->prevNormal = Vector3r(0,0,0);
 	interaction->phys = shared_ptr<ViscElPhys>(phys);

=== modified file 'py/_utils.cpp'
--- py/_utils.cpp	2010-10-01 22:30:00 +0000
+++ py/_utils.cpp	2010-10-22 14:30:53 +0000
@@ -168,10 +168,10 @@
 {
 	return Shop::inscribedCircleCenter(v0,v1,v2);
 }
-py::dict getViscoelasticFromSpheresInteraction(Real m, Real tc, Real en, Real es)
+py::dict getViscoelasticFromSpheresInteraction(Real tc, Real en, Real es)
 {
 	shared_ptr<ViscElMat> b = shared_ptr<ViscElMat>(new ViscElMat());
-	Shop::getViscoelasticFromSpheresInteraction(m,tc,en,es,b);
+	Shop::getViscoelasticFromSpheresInteraction(tc,en,es,b);
 	py::dict d;
 	d["kn"]=b->kn;
 	d["cn"]=b->cn;
@@ -465,7 +465,7 @@
 	py::def("wireNone",wireNone,"Set :yref:`Shape::wire` on all bodies to False, rendering them as solids.");
 	py::def("wireNoSpheres",wireNoSpheres,"Set :yref:`Shape::wire` to True on non-spherical bodies (:yref:`Facets<Facet>`, :yref:`Walls<Wall>`).");
 	py::def("flipCell",&Shop::flipCell,(py::arg("flip")=Matrix3r(Matrix3r::Zero())),"Flip periodic cell so that angles between $R^3$ axes and transformed axes are as small as possible. This function relies on the fact that periodic cell defines by repetition or its corners regular grid of points in $R^3$; however, all cells generating identical grid are equivalent and can be flipped one over another. This necessiatates adjustment of :yref:`Interaction.cellDist` for interactions that cross boundary and didn't before (or vice versa), and re-initialization of collider. The *flip* argument can be used to specify desired flip: integers, each column for one axis; if zero matrix, best fit (minimizing the angles) is computed automatically.\n\nIn c++, this function is accessible as ``Shop::flipCell``.\n\n.. warning::\n\t This function is currently broken and should not be used.");
-	py::def("getViscoelasticFromSpheresInteraction",getViscoelasticFromSpheresInteraction,(py::arg("m"),py::arg("tc"),py::arg("en"),py::arg("es")),"Get viscoelastic interaction parameters from analytical solution of a pair spheres collision problem. \n\n:Parameters:\n\t`m` : float\n\t\tsphere mass\n\t`tc` : float\n\t\tcollision time\n\t`en` : float\n\t\tnormal restitution coefficient\n\t`es` : float\n\t\ttangential restitution coefficient.\n\n\n:return: \n\tdict with keys:\n\n\tkn : float\n\t\tnormal elastic coefficient computed as:\n\n.. math:: k_n=\\frac{m}{t_c^2}\\left(\\pi^2+(\\ln e_n)^2\\right)\n\n\tcn : float\n\t\tnormal viscous coefficient computed as:\n\n.. math:: c_n=-\\frac{2m}{t_c}\\ln e_n\n\n\n\tkt : float\n\t\ttangential elastic coefficient computed as:\n\n.. math:: k_t=\\frac27\\frac{m}{t_c^2}\\left(\\pi^2+(\\ln e_t)^2\\right)\n\n\tct : float\n\t\ttangential viscous coefficient computed as:\n\n.. math:: c_t=-\\frac27\\frac{m}{t_c}\\ln e_t.\n\n\nFor details see [Pournin2001]_. ");
+	py::def("getViscoelasticFromSpheresInteraction",getViscoelasticFromSpheresInteraction,(py::arg("tc"),py::arg("en"),py::arg("es")),"Get viscoelastic interaction parameters from analytical solution of a pair spheres collision problem. \n\n:Parameters:\n\t`m` : float\n\t\tsphere mass\n\t`tc` : float\n\t\tcollision time\n\t`en` : float\n\t\tnormal restitution coefficient\n\t`es` : float\n\t\ttangential restitution coefficient.\n\n\n:return: \n\tdict with keys:\n\n\tkn : float\n\t\tnormal elastic coefficient computed as:\n\n.. math:: k_n=\\frac{m}{t_c^2}\\left(\\pi^2+(\\ln e_n)^2\\right)\n\n\tcn : float\n\t\tnormal viscous coefficient computed as:\n\n.. math:: c_n=-\\frac{2m}{t_c}\\ln e_n\n\n\n\tkt : float\n\t\ttangential elastic coefficient computed as:\n\n.. math:: k_t=\\frac27\\frac{m}{t_c^2}\\left(\\pi^2+(\\ln e_t)^2\\right)\n\n\tct : float\n\t\ttangential viscous coefficient computed as:\n\n.. math:: c_t=-\\frac27\\frac{m}{t_c}\\ln e_t.\n\n\nFor details see [Pournin2001]_. ");
 	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("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.");