yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #07205
[Branch ~yade-dev/yade/trunk] Rev 2773: - Add formulation for determination of time step for elastic spheres according to Rayleigh wave s...
------------------------------------------------------------
revno: 2773
committer: Chiara Modenese <c.modenese@xxxxxxxxx>
branch nick: yade
timestamp: Mon 2011-02-28 11:32:03 +0000
message:
- Add formulation for determination of time step for elastic spheres according to Rayleigh wave speed.
modified:
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/Shop.cpp'
--- pkg/dem/Shop.cpp 2011-02-27 13:54:43 +0000
+++ pkg/dem/Shop.cpp 2011-02-28 11:32:03 +0000
@@ -346,6 +346,26 @@
}
return dt;
}
+
+/* Detremination of time step as according to Rayleigh wave speed of force propagation (see Thornton 2000, ref. MillerPursey_1955) */
+Real Shop::RayleighWaveTimeStep(const shared_ptr<Scene> _rb){
+ shared_ptr<Scene> rb=(_rb?_rb:Omega::instance().getScene());
+ Real dt=std::numeric_limits<Real>::infinity();
+ FOREACH(const shared_ptr<Body>& b, *rb->bodies){
+ if(!b || !b->material || !b->shape) continue;
+
+ shared_ptr<ElastMat> ebp=dynamic_pointer_cast<ElastMat>(b->material);
+ shared_ptr<Sphere> s=dynamic_pointer_cast<Sphere>(b->shape);
+ if(!ebp || !s) continue;
+
+ Real density=b->state->mass/((4/3.)*Mathr::PI*pow(s->radius,3));
+ Real ShearModulus=ebp->young/(2.*(1+ebp->poisson));
+ Real lambda=0.1631*ebp->poisson+0.876605;
+ dt=min(dt,Mathr::PI*s->radius/lambda*sqrt(density/ShearModulus));
+ }
+ return dt;
+}
+
/* Project 3d point into 2d using spiral projection along given axis;
* the returned tuple is
*
@@ -428,9 +448,9 @@
*/
Vector3r Shop::scalarOnColorScale(Real x, Real xmin, Real xmax){
Real xnorm=min((Real)1.,max((x-xmin)/(xmax-xmin),(Real)0.));
- if(xnorm<.25) return Vector3r(0,4.*xnorm,1);
+ if(xnorm<.25) return Vector3r(0,.4*xnorm,1);
if(xnorm<.5) return Vector3r(0,1,1.-4.*(xnorm-.25));
- if(xnorm<.75) return Vector3r(4.*(xnorm-.5),1.,0);
+ if(xnorm<.75) return Vector3r(4*(xnorm-.5),1.,0);
return Vector3r(1,1-4*(xnorm-.75),0);
}
@@ -533,7 +553,6 @@
// *** Normal stress tensor split into two parts according to subnetworks of strong and weak forces (or other distinction if a threshold value for the force is assigned) ***/
Real Fmean(0); Matrix3r f, fs, fw;
- /// FIXME: this does not compute Fmean as it should?
fabricTensor(Fmean,f,fs,fw,false,compressionPositive,NaN);
Matrix3r sigNStrong(Matrix3r::Zero()), sigNWeak(Matrix3r::Zero());
FOREACH(const shared_ptr<Interaction>& I, *scene->interactions){
=== modified file 'pkg/dem/Shop.hpp'
--- pkg/dem/Shop.hpp 2011-02-09 11:58:56 +0000
+++ pkg/dem/Shop.hpp 2011-02-28 11:32:03 +0000
@@ -58,6 +58,9 @@
//! Estimate timestep based on P-wave propagation speed
static Real PWaveTimeStep(const shared_ptr<Scene> rb=shared_ptr<Scene>());
+
+ //! Estimate timestep based on Rayleigh-wave propagation speed
+ static Real RayleighWaveTimeStep(const shared_ptr<Scene> rb=shared_ptr<Scene>());
//! return 2d coordinates of a 3d point within plane defined by rotation axis and inclination of spiral, wrapped to the 0th period
static boost::tuple<Real, Real, Real> spiralProject(const Vector3r& pt, Real dH_dTheta, int axis=2, Real periodStart=std::numeric_limits<Real>::quiet_NaN(), Real theta0=0);
=== modified file 'py/_utils.cpp'
--- py/_utils.cpp 2011-02-25 18:19:15 +0000
+++ py/_utils.cpp 2011-02-28 11:32:03 +0000
@@ -76,6 +76,7 @@
}
Real PWaveTimeStep(){return Shop::PWaveTimeStep();};
+Real RayleighWaveTimeStep(){return Shop::RayleighWaveTimeStep();};
Real elasticEnergyInAABB(py::tuple Aabb){
Vector3r bbMin=py::extract<Vector3r>(Aabb[0])(), bbMax=py::extract<Vector3r>(Aabb[1])();
@@ -440,10 +441,10 @@
YADE_SET_DOCSTRING_OPTS;
-
py::def("PWaveTimeStep",PWaveTimeStep,"Get timestep accoring to the velocity of P-Wave propagation; computed from sphere radii, rigidities and masses.");
+ py::def("RayleighWaveTimeStep",RayleighWaveTimeStep,"Determination of time step according to Rayleigh wave speed of force propagation.");
py::def("getSpheresVolume",Shop__getSpheresVolume,"Compute the total volume of spheres in the simulation (might crash for now if dynamic bodies are not spheres)");
- py::def("porosity",Shop__getPorosity,(py::arg("volume")=-1),"Compute packing porosity $\\frac{V-V_s}{V}$ where $V$ is overall volume and $V_s$ is volume of spheres.\n\n:param float volume: overall volume which must be specified for aperiodic simulations. For periodic simulations, current volume of the :yref:`Cell` is used.\n");
+ py::def("porosity",Shop__getPorosity,(py::arg("volume")=-1),"Compute packing poro sity $\\frac{V-V_s}{V}$ where $V$ is overall volume and $V_s$ is volume of spheres.\n\n:param float volume: overall volume which must be specified for aperiodic simulations. For periodic simulations, current volume of the :yref:`Cell` is used.\n");
py::def("aabbExtrema",aabbExtrema,(py::arg("cutoff")=0.0,py::arg("centers")=false),"Return coordinates of box enclosing all bodies\n\n:param bool centers: do not take sphere radii in account, only their centroids\n:param float∈〈0…1〉 cutoff: relative dimension by which the box will be cut away at its boundaries.\n\n\n:return: (lower corner, upper corner) as (Vector3,Vector3)\n\n");
py::def("ptInAABB",isInBB,"Return True/False whether the point p is within box given by its min and max corners");
py::def("negPosExtremeIds",negPosExtremeIds,negPosExtremeIds_overloads(py::args("axis","distFactor"),"Return list of ids for spheres (only) that are on extremal ends of the specimen along given axis; distFactor multiplies their radius so that sphere that do not touch the boundary coordinate can also be returned."));
