yade-dev team mailing list archive

[noreply@xxxxxxxxxxxxx]

------------------------------------------------------------
revno: 3665
committer: Raphael Maurin <raph_maurin@xxxxxxxxxxx>
timestamp: Tue 2015-06-02 19:09:57 +0200
message:
  Law2_ScGeom_ViscElPhys_Basic: add a new formulation + check
  Add a new formulation of the law when imposing in the material young, poisson, and normal restitution coefficient.
  The effective spring constants ks and kn are evaluated "classically" with the young modulus and the poisson's ratio.
  The specificity of the formulation stands in the fact that it is the restitution coefficient of the contact (en = 2*en1*en2/(en1+en2)) that is imposed, and not the damping constant. To achieve this, the damping constant of the contact is evaluated resolving numerically the analytical expression 21 of [Schwager2007] (with a minus sign, there is a mistake in the article) with a Newton-Raphson algorithm (This was made by Francois Kneib). This does not seem to affect the calculation time.
  With this formulation, there is no tangential viscous damping, and it is usually not possible to access value of normal restitution coefficient lower than 0.1 (in that case the numerical resolution does not converged in the maximum 15 iteration allowed in the code)
  Together with the new formulation, add a short check script to verify that the imposed restitution coefficient is the one obtained.
added:
  scripts/checks-and-tests/checks/checkViscElPM2.py
modified:
  pkg/dem/ViscoelasticPM.cpp
  pkg/dem/ViscoelasticPM.hpp


--
lp:yade
https://code.launchpad.net/~yade-pkg/yade/git-trunk

Your team Yade developers is subscribed to branch lp:yade.
To unsubscribe from this branch go to https://code.launchpad.net/~yade-pkg/yade/git-trunk/+edit-subscription

=== modified file 'pkg/dem/ViscoelasticPM.cpp'
--- pkg/dem/ViscoelasticPM.cpp	2015-03-16 10:41:06 +0000
+++ pkg/dem/ViscoelasticPM.cpp	2015-06-02 17:09:57 +0000
@@ -19,6 +19,8 @@
 /* ViscElPhys */
 ViscElPhys::~ViscElPhys(){}
 
+Real Ip2_ViscElMat_ViscElMat_ViscElPhys::epsilon = 1.0e-15;
+
 /* 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 
@@ -227,6 +229,11 @@
 			cs1 = mat1->cs;
 			cs2 = mat2->cs;
 		}
+		else if( isfinite(mat1->en) and isfinite(mat1->et)) {
+			const Real En = (en) ? (*en)(mat1->id,mat2->id) : (mat1->en+mat2->en)/2.0;
+			cn1 = cn2 = 2.0*find_cn_from_en(En, massR,contactParameterCalculation(kn1,kn2),interaction);
+			cs1 = cs2 = 0;
+		}
 	}
 	
 	const Real mR1 = mat1->mR;      const Real mR2 = mat2->mR; 
@@ -276,3 +283,40 @@
   else return 0;
 }
 
+Real find_cn_from_en(const Real& en, const Real& m, const Real& kn, const shared_ptr<Interaction>& interaction){
+	Real eps = Ip2_ViscElMat_ViscElMat_ViscElPhys::epsilon;
+	Real cn = eps ; //initial small value
+	Real en_temp=get_en_from_cn(cn,m ,kn);
+	int i =0;
+	Real error = 1.0/eps;
+	while (error > 1.0e-2){
+		if(i>15){
+			cerr<<"Warning in ViscoelasticPM.cpp : Newton-Raphson algorithm did not converged within 15 iterations for contact between "<<interaction->id1<<" and "<<interaction->id2<<". Continue with values : cn="<<cn<<" en="<<en_temp<<endl;
+			break;
+		}
+		i++;
+		Real deriv=(get_en_from_cn(cn-eps,m ,kn)-get_en_from_cn(cn+eps,m ,kn))/(-2.*eps);
+		cn = cn - (en_temp-en)/deriv;
+		en_temp=get_en_from_cn(cn,m ,kn);
+		error = fabs(en_temp-en)/en;
+	}
+// 	cout<<"i="<<i<<" error="<<error<<endl;
+	return cn;
+}
+
+Real get_en_from_cn(const Real& cn, const Real& m, const Real& kn){
+	Real beta = 0.5*cn/m;
+	Real omega0 = sqrt(kn/m);
+	Real omega = sqrt(omega0*omega0-beta*beta);
+	Real Omega = sqrt(beta*beta-omega0*omega0);
+	if ( beta < omega0/sqrt(2.) )
+		return exp(-beta/omega*(Mathr::PI-atan(2.*beta*omega/(omega*omega-beta*beta))));
+	else if ( beta > omega0/sqrt(2.) and  beta < omega0)
+		return exp(-beta/omega*atan(-2.*beta*omega/(omega*omega-beta*beta)));
+	else if ( beta > omega0 )
+		return exp(-beta/Omega*log((beta+Omega)/(beta-Omega)));
+	else if ( beta == omega0/sqrt(2.) or beta == omega0 )
+		return get_en_from_cn(cn + Ip2_ViscElMat_ViscElMat_ViscElPhys::epsilon, m, kn);
+	else return 0;
+}
+

=== modified file 'pkg/dem/ViscoelasticPM.hpp'
--- pkg/dem/ViscoelasticPM.hpp	2015-04-29 17:12:57 +0000
+++ pkg/dem/ViscoelasticPM.hpp	2015-06-02 17:09:57 +0000
@@ -76,6 +76,7 @@
 // Uses the rule of consecutively connection.
 class Ip2_ViscElMat_ViscElMat_ViscElPhys: public IPhysFunctor {
 	public :
+		static Real epsilon;
 		virtual void go(const shared_ptr<Material>& b1,
 					const shared_ptr<Material>& b2,
 					const shared_ptr<Interaction>& interaction);
@@ -104,3 +105,6 @@
 
 Real contactParameterCalculation(const Real& l1,const Real& l2);
 bool computeForceTorqueViscEl(shared_ptr<IGeom>& _geom, shared_ptr<IPhys>& _phys, Interaction* I, Vector3r & force, Vector3r & torque1, Vector3r & torque2);
+
+Real get_en_from_cn(const Real& cn, const Real& m, const Real& kn);
+Real find_cn_from_en(const Real& en, const Real& m, const Real& kn, const shared_ptr<Interaction>& interaction);

=== added file 'scripts/checks-and-tests/checks/checkViscElPM2.py'
--- scripts/checks-and-tests/checks/checkViscElPM2.py	1970-01-01 00:00:00 +0000
+++ scripts/checks-and-tests/checks/checkViscElPM2.py	2015-06-02 17:09:57 +0000
@@ -0,0 +1,76 @@
+#!/usr/bin/env python
+# encoding: utf-8
+
+################################################################################
+# CONSTITUTIVE LAW TESTING: Law2_ScGeom_ViscElPhys_Basic()
+#
+# Two spheres of different diameter and young modulus, with velocities 
+# (v,0,0) and (-v,0,0) collide.
+# The script checks if the restitution coefficient applied is the one obtained
+# Law2_ScGeom_ViscElPhys_Basic() compute the stiffness and damping coefficient
+# differently depending on the input. This test is then different from 
+# checkViscElPM.py and checkViscElEng.py
+#
+
+from yade import plot,ymport
+import math
+
+
+################################################################################
+# MATERIAL
+diameter1 = 1e-2 #[m] diameter of the first particle
+diameter2 = 1e-3    #[m] diameter of the second particle
+rho1 = 1500 # [kg/m^3] density of the first particle
+rho2 = 2500 # [kg/m^3] density of the second particle
+youngMod = 1e5 # [kg/m/s^2] young modulus of the particles
+poissonRatio = 0.5 # [-] poisson's ratio of the particles
+mu  = 0.4 # [-]      friction coefficient of both particle
+en  = 0.5  # [-]     normal restitution coefficient of both particle
+et  = 1.  # [-]      tangential restitution coefficient of both particle, no damping
+
+frictionAngle = math.atan(mu)
+
+O.materials.append(ViscElMat(en=en,et=1.,young = youngMod,poisson = poissonRatio,frictionAngle=frictionAngle,density=rho1, label='Mat1'))
+O.materials.append(ViscElMat(en=en,et=1.,young = youngMod,poisson = poissonRatio,frictionAngle=frictionAngle,density=rho2, label='Mat2'))
+
+################################################################################
+v = 1e-1
+# GEOMETRY
+O.bodies.append(sphere((0,0,0),diameter1/2.0, material = 'Mat1'))
+O.bodies[0].state.vel[2] = v
+O.bodies.append(sphere((0,0,(diameter1+diameter2)/2.+ 0.01*diameter1),diameter2/2.0, material = 'Mat2'))
+O.bodies[-1].state.vel[2] = -v
+
+
+################################################################################
+# SIMULATION
+
+O.dt = 1e-6 # [s]   fixed time step
+
+O.engines=[
+  ForceResetter(),
+  InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb(),Bo1_Wall_Aabb()]),
+  InteractionLoop(
+    [Ig2_Sphere_Sphere_ScGeom(),
+     Ig2_Facet_Sphere_ScGeom(),
+     Ig2_Wall_Sphere_ScGeom()],
+    [Ip2_ViscElMat_ViscElMat_ViscElPhys()],
+    [Law2_ScGeom_ViscElPhys_Basic()]
+  ),
+  NewtonIntegrator(damping=0.),
+  PyRunner(command = 'check()',iterPeriod = int(0.05/O.dt)-1)
+]
+
+
+################################################################################
+# RUN
+O.saveTmp()
+O.run(int(0.05/O.dt))
+################################################################################
+def check():
+	# Compare imposed restitution coefficient and obtained one
+	enMeasured = (O.bodies[1].state.vel[2]-O.bodies[0].state.vel[2])/(2*v)
+	tolerance = 1e-2
+	if (abs(enMeasured -en)/en)<tolerance:
+		resultStatus+=1
+