yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #11807
[Branch ~yade-pkg/yade/git-trunk] Rev 3573: Add new function to compute stress tensor depth profile.
------------------------------------------------------------
revno: 3573
committer: Raphael Maurin <raph_maurin@xxxxxxxxxxx>
timestamp: Fri 2015-02-06 18:07:21 +0100
message:
Add new function to compute stress tensor depth profile.
Compute and return the stress tensor for each cell of height dz, between two limit points given by the user.
The stress tensor computed include contributions from both particles inertia and Love-Weber stress tensor. This last take into account only the part of the branch vector contained in the cell.
modified:
pkg/dem/Shop.hpp
pkg/dem/Shop_02.cpp
py/_utils.cpp
--
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/Shop.hpp'
--- pkg/dem/Shop.hpp 2014-10-15 06:44:01 +0000
+++ pkg/dem/Shop.hpp 2015-02-06 17:07:21 +0000
@@ -121,6 +121,9 @@
static Matrix3r getStress(Real volume=0);
static Matrix3r getCapillaryStress(Real volume=0, bool mindlin=false);
static Matrix3r stressTensorOfPeriodicCell() { LOG_WARN("Shop::stressTensorOfPeriodicCelli is DEPRECATED: use getStress instead"); return Shop::getStress(); }
+
+ static vector<Matrix3r> getStressProfile(Real volume, int nCell, Real dz, Real zRef, vector<Real> vPartAverageX, vector<Real> vPartAverageY, vector<Real> vPartAverageZ);
+
//! Compute overall ("macroscopic") stress of periodic cell, returning 2 tensors
//! (contribution of normal and shear forces)
static py::tuple normalShearStressTensors(bool compressionPositive=false, bool splitNormalTensor=false, Real thresholdForce=NaN);
=== modified file 'pkg/dem/Shop_02.cpp'
--- pkg/dem/Shop_02.cpp 2014-10-15 15:51:34 +0000
+++ pkg/dem/Shop_02.cpp 2015-02-06 17:07:21 +0000
@@ -372,6 +372,106 @@
return stressTensor/volume;
}
+
+
+vector<Matrix3r> Shop::getStressProfile(Real volume, int nCell, Real dz, Real zRef, vector<Real> vPartAverageX, vector<Real> vPartAverageY, vector<Real> vPartAverageZ){
+ int minZ;
+ int maxZ;
+ Real minPosZ;
+ Real maxPosZ;
+ Scene* scene=Omega::instance().getScene().get();
+ vector<Matrix3r> stressTensorProfile(nCell,Matrix3r::Zero());
+ const bool isPeriodic = scene->isPeriodic;
+
+ //
+ //Dynamic contribution to the stress tensor
+ //
+ FOREACH(const shared_ptr<Body>& b,*Omega::instance().getScene()->bodies){
+ int Np = floor((b->state->pos[2]-zRef)/dz); //Define the layer number with 0 corresponding to zRef
+ if ((Np>=0)&&(Np<nCell)){ //To avoid non defined vPartAverage
+ //Velocity fluctuation wrt the average field
+ Vector3r vFluct = b->state->vel - Vector3r(vPartAverageX[Np],vPartAverageY[Np],vPartAverageZ[Np]);
+ //Classical dynamical expression of the stress tensor
+ stressTensorProfile[Np]+= -1/volume*b->state->mass*vFluct*vFluct.transpose();
+ }
+ }
+
+ //
+ //Love Weber contribution (same as getStress(), but with different layers)
+ //
+ FOREACH(const shared_ptr<Interaction>&I, *scene->interactions){ // loop over the interactions
+ if (!I->isReal()) continue;
+ shared_ptr<Body> b1 = Body::byId(I->getId1(),scene);
+ shared_ptr<Body> b2 = Body::byId(I->getId2(),scene);
+
+ if ((b1->state->blockedDOFs!=State::DOF_ALL)||(b2->state->blockedDOFs!=State::DOF_ALL)){// to remove annoying contribution from the fixed particles
+ //Layers in which the particle center is contained
+ int Np1 = floor((b1->state->pos[2] - zRef)/dz);
+ int Np2 = floor((b2->state->pos[2] - zRef)/dz);
+ //Vector between the two centers, from 2 to 1
+ Vector3r branch = b1->state->pos -b2->state->pos;
+ if (isPeriodic) branch -= scene->cell->hSize*I->cellDist.cast<Real>();//to handle periodicity
+
+ //Contact vector (from 1 to 2)
+ NormShearPhys* nsi=YADE_CAST<NormShearPhys*> ( I->phys.get() );
+ Vector3r fContact = nsi->normalForce + nsi->shearForce;
+
+ //The contribution to the stress tensor is taken such that only the part of the branch vector
+ //inside the cell is taken into account
+ //If the whole vector is in the cell, add the whole contribution to the cell
+ if (Np1==Np2){
+ if ((Np1>=0) && (Np1<nCell)){
+ stressTensorProfile[Np1]+= 1/volume*fContact*branch.transpose();
+ }
+ }
+ //Otherwise, find out the cell crossed by the branch vector and assign it the contribution from this part.
+ else {
+ //Find which one is above the other to prepare the loop
+ if (Np1>Np2){
+ minZ = Np2;
+ minPosZ = b2->state->pos[2]-zRef;
+ maxZ = Np1;
+ maxPosZ = b1->state->pos[2]-zRef;
+ }
+ else if (Np2>Np1) {
+ minZ = Np1;
+ minPosZ = b1->state->pos[2]-zRef;
+ maxZ = Np2;
+ maxPosZ = b2->state->pos[2]-zRef;
+ }
+
+ Real branchS_x = pow(branch[0],2.0);
+ Real branchS_y = pow(branch[1],2.0);
+ Real branchS_z = pow(branch[2],2.0);
+
+ //Normalize the branch vector
+ branch/=sqrt(branchS_x + branchS_y+branchS_z);
+
+ //Loop over the cell containing the branch vector
+ int numLayer = minZ;
+ while (numLayer<=maxZ){
+ if ((numLayer>= 0)&&(numLayer<nCell)){
+ //Evaluate the branch height inside the cell
+ Real deltaZ = dz;
+ if (numLayer==minZ) deltaZ = dz - (minPosZ - minZ*dz);
+ else if (numLayer==maxZ) deltaZ = maxPosZ - maxZ*dz;
+ //From it, trigonometry gives us the vector contained in the cell
+ Vector3r branchVectCell = deltaZ*sqrt(1.0 + 1.0/branchS_z*(branchS_x + branchS_y))*branch;
+ //Add the contribution to the stress tensor
+ stressTensorProfile[numLayer]+= 1.0/volume*fContact*branchVectCell.transpose();
+ }
+ //Increment the layer/cell number
+ numLayer+=1;
+
+ }
+ }
+ }
+ }
+ return stressTensorProfile;
+}
+
+
+
Matrix3r Shop::getCapillaryStress(Real volume, bool mindlin){
Scene* scene=Omega::instance().getScene().get();
if (volume==0) volume = scene->isPeriodic?scene->cell->hSize.determinant():1;
=== modified file 'py/_utils.cpp'
--- py/_utils.cpp 2014-12-01 18:29:59 +0000
+++ py/_utils.cpp 2015-02-06 17:07:21 +0000
@@ -455,6 +455,7 @@
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,"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_{ik}\\sigma_{kj}=x_{i,k}\\sigma_{kj}=(x_{i}\\sigma_{kj})_{,k}-x_{i}\\sigma_{kj,k}$.\n\nAt equilibrium, the divergence of stress is null: $\\sigma_{kj,k}=\\vec{0}$. Consequently, after divergence theorem: $\\frac{1}{V}\\int_V \\sigma_{ij}dV = \\frac{1}{V}\\int_V (x_{i}\\sigma_{kj})_{,k}dV = \\frac{1}{V}\\int_{\\partial V}x_i\\sigma_{kj}n_kdS = \\frac{1}{V}\\sum_bx_i^bf_j^b$.\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\nThis last relation may not be very useful if porosity is not homogeneous. If it happens, one can define the equivalent bulk stress a the particles scale by assigning a volume to each particle. This volume can be obtained from :yref:`TesselationWrapper` (see e.g. [Catalano2014a]_)");
py::def("getStress",Shop::getStress,(py::args("volume")=0),"Compute and return Love-Weber stress tensor:\n\n $\\sigma_{ij}=\\frac{1}{V}\\sum_b f_i^b l_j^b$, where the sum is over all interactions, with $f$ the contact force and $l$ the branch vector (joining centers of the bodies). Stress is negativ for repulsive contact forces, i.e. compression. $V$ can be passed to the function. If it is not, it will be equal to the volume of the cell in periodic cases, or to the one deduced from utils.aabbDim() in non-periodic cases.");
+ py::def("getStressProfile",Shop::getStressProfile,(py::args("volume"),py::args("nCell"),py::args("dz"),py::args("vPartAverageX"),py::args("vPartAverageY"),py::args("vPartAverageZ"),py::args("zRef")),"Compute and return the stress tensor depth profile, including the contribution from Love-Weber stress tensor and the dynamic stress tensor taking into account the effect of particles inertia. For each defined cell z, the stress tensor reads: \n\n $\\sigma_{ij}^z= \\frac{1}{V}\\sum_c f_i^c l_j^{c,z} - \\frac{1}{V}\\sum_p m^p u'^p_i u'^p_j$,\n\n where the first sum is made over the contacts which are contained or cross the cell z, f^c is the contact force from particle 1 to particle 2, and l^{c,z} is the part of the branch vector from particle 2 to particle 1, contained in the cell. The second sum is made over the particles, and u'^p is the velocity fluctuations of the particle p with respect to the spatial averaged particle velocity at this point (given as input parameters). The expression of the stress tensor is the same as the one given in getStress plus the inertial contribution. Apart from that, the main difference with getStress stands in the fact that it gives a depth profile of stress tensor, i.e. from the reference horizontal plane at elevation zRef (input parameter) until the plane of elevation zRef+nCell*dz (input parameters), it is computing the stress tensor for each cell of height dz. For the love-Weber stress contribution, the branch vector taken into account in the calculations is only the part of the branch vector contained in the cell considered.\n To validate the formulation, it has been checked that activating only the Love-Weber stress tensor, and suming all the contributions at the different altitude, we recover the same stress tensor as when using getStress. For my own use, I have troubles with strong overlap between fixed object, so that I made a condition to exclude the contribution to the stress tensor of the fixed objects, this can be desactivated easily if needed (and should be desactivated for the comparison with getStress)." );
py::def("getCapillaryStress",Shop::getCapillaryStress,(py::args("volume")=0,py::args("mindlin")=false),"Compute and return Love-Weber capillary stress tensor:\n\n $\\sigma^{cap}_{ij}=\\frac{1}{V}\\sum_b l_i^b f^{cap,b}_j$, where the sum is over all interactions, with $l$ the branch vector (joining centers of the bodies) and $f^{cap}$ is the capillary force. $V$ can be passed to the function. If it is not, it will be equal to one in non-periodic cases, or equal to the volume of the cell in periodic cases. Only the CapillaryPhys interaction type is supported presently. Using this function with physics MindlinCapillaryPhys needs to pass True as second argument.");
py::def("getBodyIdsContacts",Shop__getBodyIdsContacts,(py::args("bodyID")=0),"Get a list of body-ids, which contacts the given body.");
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.");
