yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #11674
[Branch ~yade-pkg/yade/git-trunk] Rev 3541: Many changes towards a more consistent sign convention throughout the code for stresses and strai...
------------------------------------------------------------
revno: 3541
committer: jduriez <jerome.duriez@xxxxxxxxxxx>
timestamp: Tue 2014-11-18 12:46:27 -0700
message:
Many changes towards a more consistent sign convention throughout the code for stresses and strains. Should fix e.g. https://bugs.launchpad.net/yade/+bug/1381282. The goal is to know directly the meaning of a strain or stress value each time one is encountered. Relying on the classical Continuum Mechanics convention (editorial choice...). A global announcement is planned in a couple of days, in case some errors appear in the meantime.
modified:
examples/FluidCouplingPFV/oedometer.py
examples/clumps/triax-basic-with-clumps.py
examples/test/triax-basic.py
examples/test/triax-cohesive.py
examples/triax-tutorial/script-session1.py
pkg/dem/CapillaryTriaxialTest.hpp
pkg/dem/CohesiveTriaxialTest.hpp
pkg/dem/TriaxialCompressionEngine.cpp
pkg/dem/TriaxialCompressionEngine.hpp
pkg/dem/TriaxialStressController.cpp
pkg/dem/TriaxialStressController.hpp
pkg/dem/TriaxialTest.hpp
py/pack/pack.py
scripts/checks-and-tests/checks/DEM-PFV-check.py
scripts/checks-and-tests/checks/checkTestTriax.py
--
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 'examples/FluidCouplingPFV/oedometer.py'
--- examples/FluidCouplingPFV/oedometer.py 2014-04-01 13:18:38 +0000
+++ examples/FluidCouplingPFV/oedometer.py 2014-11-18 19:46:27 +0000
@@ -64,12 +64,12 @@
newton
]
-triax.goal1=triax.goal2=triax.goal3=10000
+triax.goal1=triax.goal2=triax.goal3=-10000
while 1:
O.run(1000, True)
unb=unbalancedForce()
- if unb<0.001 and abs(10000-triax.meanStress)/10000<0.001:
+ if unb<0.001 and abs(-10000-triax.meanStress)/10000<0.001:
break
setContactFriction(radians(finalFricDegree))
@@ -83,14 +83,14 @@
triax.internalCompaction=False
triax.wall_bottom_activated=False
#load
-triax.goal2=11000; O.run(2000,1)
+triax.goal2=-11000; O.run(2000,1)
#unload
-triax.goal2=10000; O.run(2000,1)
+triax.goal2=-10000; O.run(2000,1)
#load
-triax.goal2=11000; O.run(2000,1)
+triax.goal2=-11000; O.run(2000,1)
e22=triax.strain[1]
#unload
-triax.goal2=10000; O.run(2000,1)
+triax.goal2=-10000; O.run(2000,1)
e22=e22-triax.strain[1]
modulus = 1000./abs(e22)
@@ -120,12 +120,12 @@
newton.damping=0
#we want the theoretical value from Terzaghi's solution
-#keep in mind that we are not in an homogeneous material and the smal strain
+#keep in mind that we are not in an homogeneous material and the small strain
#assumption is not verified => we don't expect perfect match
#there can be also an overshoot of pressure in the very beginning due to dynamic effects
Cv=permeability*modulus/1e4
zeroTime=O.time
-zeroe22 = triax.strain[1]
+zeroe22 = - triax.strain[1]
dryFraction=0.05 #the top layer is affected by drainage on a certain depth, we account for it here
drye22 = 1000/modulus*dryFraction
wetHeight=1*(1-dryFraction)
@@ -137,7 +137,7 @@
U=U-2/M**2*exp(-M**2*Tv)
return U
-triax.goal2=11000
+triax.goal2=-11000
@@ -145,8 +145,8 @@
## a function saving variables
def history():
- plot.addData(e22=triax.strain[1]-zeroe22,e22_theory=drye22+(1-dryFraction)*consolidation((O.time-zeroTime)*Cv/wetHeight**2)*1000./modulus,t=O.time,p=flow.getPorePressure((0.5,0.1,0.5)),s22=triax.stress(3)[1]-10000)
- #plot.addData(e22=triax.strain[1],t=O.time,s22=-triax.stress(2)[1],p=flow.MeasurePorePressure((0.5,0.5,0.5)))
+ plot.addData(e22=-triax.strain[1]-zeroe22,e22_theory=drye22+(1-dryFraction)*consolidation((O.time-zeroTime)*Cv/wetHeight**2)*1000./modulus,t=O.time,p=flow.getPorePressure((0.5,0.1,0.5)),s22=-triax.stress(3)[1]-10000)
+ #plot.addData(e22=-triax.strain[1],t=O.time,s22=-triax.stress(2)[1],p=flow.MeasurePorePressure((0.5,0.5,0.5)))
O.engines=O.engines+[PyRunner(iterPeriod=200,command='history()',label='recorder')]
##make nice animations:
=== modified file 'examples/clumps/triax-basic-with-clumps.py'
--- examples/clumps/triax-basic-with-clumps.py 2013-03-28 11:01:40 +0000
+++ examples/clumps/triax-basic-with-clumps.py 2014-11-18 19:46:27 +0000
@@ -39,8 +39,8 @@
wall_back_id=wallIds[4],
wall_front_id=wallIds[5],
internalCompaction=False,
- sigmaIsoCompaction=50e3,
- sigmaLateralConfinement=50e3,
+ sigmaIsoCompaction=-50e3,
+ sigmaLateralConfinement=-50e3,
strainRate=0.02,
frictionAngleDegree=30,
StabilityCriterion = 0.01,
=== modified file 'examples/test/triax-basic.py'
--- examples/test/triax-basic.py 2013-03-28 11:13:12 +0000
+++ examples/test/triax-basic.py 2014-11-18 19:46:27 +0000
@@ -32,8 +32,8 @@
## define the rest of triax params here
## see in pkg/dem/PreProcessor/TriaxialTest.cpp:524 etc
## which are assigned in the c++ preprocessor actually
- sigmaIsoCompaction=50e3,
- sigmaLateralConfinement=50e3,
+ sigmaIsoCompaction=-50e3,
+ sigmaLateralConfinement=-50e3,
max_vel=10,
strainRate=0.01,
label="triax"
@@ -56,7 +56,7 @@
from yade import plot
O.engines=O.engines[0:5]+[PyRunner(iterPeriod=20,command='history()',label='recorder')]+O.engines[5:7]
def history():
- plot.addData(e11=triax.strain[0], e22=triax.strain[1], e33=triax.strain[2],
+ plot.addData(e11=-triax.strain[0], e22=-triax.strain[1], e33=-triax.strain[2],
s11=-triax.stress(0)[0],
s22=-triax.stress(2)[1],
s33=-triax.stress(4)[2],
=== modified file 'examples/test/triax-cohesive.py'
--- examples/test/triax-cohesive.py 2013-03-28 11:13:12 +0000
+++ examples/test/triax-cohesive.py 2014-11-18 19:46:27 +0000
@@ -30,8 +30,8 @@
wall_back_id=wallIds[4],
wall_front_id=wallIds[5],
internalCompaction=False,
- sigmaIsoCompaction=50e3,
- sigmaLateralConfinement=50e3,
+ sigmaIsoCompaction=-50e3,
+ sigmaLateralConfinement=-50e3,
max_vel=10,
strainRate=0.03,
label="triax"
@@ -60,7 +60,7 @@
from yade import plot
O.engines=O.engines[0:5]+[PyRunner(iterPeriod=20,command='history()',label='recorder')]+O.engines[5:7]
def history():
- plot.addData(e11=O.engines[4].strain[0], e22=O.engines[4].strain[1], e33=O.engines[4].strain[2],
+ plot.addData(e11=-O.engines[4].strain[0], e22=-O.engines[4].strain[1], e33=-O.engines[4].strain[2],
s11=-O.engines[4].stress(0)[0],
s22=-O.engines[4].stress(2)[1],
s33=-O.engines[4].stress(4)[2],
=== modified file 'examples/triax-tutorial/script-session1.py'
--- examples/triax-tutorial/script-session1.py 2014-09-22 23:13:24 +0000
+++ examples/triax-tutorial/script-session1.py 2014-11-18 19:46:27 +0000
@@ -39,7 +39,7 @@
targetPorosity = 0.43 #the porosity we want for the packing
compFricDegree = table.compFricDegree # initial contact friction during the confining phase (will be decreased during the REFD compaction process)
finalFricDegree = 30 # contact friction during the deviatoric loading
-rate=0.02 # loading rate (strain rate)
+rate=-0.02 # loading rate (strain rate)
damp=0.2 # damping coefficient
stabilityThreshold=0.01 # we test unbalancedForce against this value in different loops (see below)
young=5e6 # contact stiffness
@@ -125,14 +125,14 @@
#######################################
#the value of (isotropic) confining stress defines the target stress to be applied in all three directions
-triax.goal1=triax.goal2=triax.goal3=10000
+triax.goal1=triax.goal2=triax.goal3=-10000
#while 1:
#O.run(1000, True)
##the global unbalanced force on dynamic bodies, thus excluding boundaries, which are not at equilibrium
#unb=unbalancedForce()
#print 'unbalanced force:',unb,' mean stress: ',triax.meanStress
- #if unb<stabilityThreshold and abs(10000-triax.meanStress)/10000<0.001:
+ #if unb<stabilityThreshold and abs(-10000-triax.meanStress)/10000<0.001:
#break
#O.save('confinedState'+key+'.yade.gz')
@@ -174,10 +174,10 @@
##set stress control on x and z, we will impose strain rate on y
#triax.stressMask = 5
##now goal2 is the target strain rate
-#triax.goal2=-rate
+#triax.goal2=rate
## we define the lateral stresses during the test, here the same 10kPa as for the initial confinement.
-#triax.goal1=10000
-#triax.goal3=10000
+#triax.goal1=-10000
+#triax.goal3=-10000
##we can change damping here. What is the effect in your opinion?
#newton.damping=0.1
@@ -193,11 +193,11 @@
### a function saving variables
#def history():
- #plot.addData(e11=triax.strain[0], e22=triax.strain[1], e33=triax.strain[2],
+ #plot.addData(e11=-triax.strain[0], e22=-triax.strain[1], e33=-triax.strain[2],
#ev=-triax.strain[0]-triax.strain[1]-triax.strain[2],
- #s11=triax.stress(triax.wall_right_id)[0],
- #s22=triax.stress(triax.wall_top_id)[1],
- #s33=triax.stress(triax.wall_front_id)[2],
+ #s11=-triax.stress(triax.wall_right_id)[0],
+ #s22=-triax.stress(triax.wall_top_id)[1],
+ #s33=-triax.stress(triax.wall_front_id)[2],
#i=O.iter)
#if 1:
=== modified file 'pkg/dem/CapillaryTriaxialTest.hpp'
--- pkg/dem/CapillaryTriaxialTest.hpp 2014-10-15 06:44:01 +0000
+++ pkg/dem/CapillaryTriaxialTest.hpp 2014-11-18 19:46:27 +0000
@@ -109,8 +109,8 @@
((Real,dampingMomentum,0.2,,"Coefficient of Cundal-Non-Viscous damping (applied on on the 3 components of torques)"))
((Real,StabilityCriterion,0.01,,"Value of unbalanced force for which the system is considered stable. Used in conditionals to switch between loading stages."))
((Real,wallOversizeFactor,1.3,,"Make boundaries larger than the packing to make sure spheres don't go out during deformation."))
- ((Real,sigmaIsoCompaction,50000,,"Confining stress during isotropic compaction."))
- ((Real,sigmaLateralConfinement,50000,,"Lateral stress during triaxial loading. An isotropic unloading is performed if the value is not equal to :yref:`CapillaryTriaxialTest::SigmaIsoCompaction`."))
+ ((Real,sigmaIsoCompaction,-50000,,"Confining stress during isotropic compaction (< 0 for real - compressive - compaction)."))
+ ((Real,sigmaLateralConfinement,-50000,,"Lateral stress during triaxial loading (< 0 for classical compressive cases). An isotropic unloading is performed if the value is not equal to :yref:`CapillaryTriaxialTest::SigmaIsoCompaction`."))
((int,timeStepUpdateInterval,50,,"interval for :yref:`GlobalStiffnessTimeStepper`"))
((int,timeStepOutputInterval,50,,"interval for outputing general information on the simulation (stress,unbalanced force,...)"))
=== modified file 'pkg/dem/CohesiveTriaxialTest.hpp'
--- pkg/dem/CohesiveTriaxialTest.hpp 2014-10-15 06:44:01 +0000
+++ pkg/dem/CohesiveTriaxialTest.hpp 2014-11-18 19:46:27 +0000
@@ -97,8 +97,8 @@
((Real,dampingMomentum,0.2,,"Coefficient of Cundal-Non-Viscous damping (applied on on the 3 components of torques)"))
((Real,StabilityCriterion,0.01,,"Value of unbalanced force for which the system is considered stable. Used in conditionals to switch between loading stages."))
((Real,wallOversizeFactor,1.3,,"Make boundaries larger than the packing to make sure spheres don't go out during deformation."))
- ((Real,sigmaIsoCompaction,50000,,"Confining stress during isotropic compaction."))
- ((Real,sigmaLateralConfinement,50000,,"Lateral stress during triaxial loading. An isotropic unloading is performed if the value is not equal to :yref:`TriaxialTest::sigmaIsoCompaction`."))
+ ((Real,sigmaIsoCompaction,-50000,,"Confining stress during isotropic compaction (< 0 for real - compressive - compaction)."))
+ ((Real,sigmaLateralConfinement,-50000,,"Lateral stress during triaxial loading (< 0 for classical compressive cases). An isotropic unloading is performed if the value is not equal to :yref:`TriaxialTest::sigmaIsoCompaction`."))
((int,timeStepUpdateInterval,50,,"interval for :yref:`GlobalStiffnessTimeStepper`"))
((int,wallStiffnessUpdateInterval,10,,"interval for updating the stiffness of sample/boundaries contacts"))
=== modified file 'pkg/dem/TriaxialCompressionEngine.cpp'
--- pkg/dem/TriaxialCompressionEngine.cpp 2014-10-15 06:44:01 +0000
+++ pkg/dem/TriaxialCompressionEngine.cpp 2014-11-18 19:46:27 +0000
@@ -121,7 +121,7 @@
if ( scene->iter % testEquilibriumInterval == 0 )
{
updateParameters ();
- maxStress = max(maxStress,stress[wall_top][1]);
+ maxStress = max(maxStress,-stress[wall_top][1]);
LOG_INFO("UnbalancedForce="<< UnbalancedForce<<", rel stress "<< std::abs ( ( meanStress-sigma_iso ) /sigma_iso ));
}
if ( saveSimulation )
=== modified file 'pkg/dem/TriaxialCompressionEngine.hpp'
--- pkg/dem/TriaxialCompressionEngine.hpp 2014-10-15 06:44:01 +0000
+++ pkg/dem/TriaxialCompressionEngine.hpp 2014-11-18 19:46:27 +0000
@@ -98,9 +98,9 @@
((int,testEquilibriumInterval,20,,"interval of checks for transition between phases, higher than 1 saves computation time."))
((stateNum,currentState,1,,"There are 5 possible states in which TriaxialCompressionEngine can be. See above :yref:`yade.wrapper.TriaxialCompressionEngine` "))
((stateNum,previousState,1,,"Previous state (used to detect manual changes of the state in .xml)"))
- ((Real,sigmaIsoCompaction,1,,"Prescribed isotropic pressure during the compaction phase"))
+ ((Real,sigmaIsoCompaction,1,,"Prescribed isotropic pressure during the compaction phase (< 0 for real - compressive - compaction)"))
((Real,previousSigmaIso,1,,"Previous value of inherited sigma_iso (used to detect manual changes of the confining pressure)"))
- ((Real,sigmaLateralConfinement,1,,"Prescribed confining pressure in the deviatoric loading; might be different from :yref:`TriaxialCompressionEngine::sigmaIsoCompaction`"))
+ ((Real,sigmaLateralConfinement,1,,"Prescribed confining pressure in the deviatoric loading (< 0 for classical compressive cases); might be different from :yref:`TriaxialCompressionEngine::sigmaIsoCompaction`"))
((std::string,Key,"",,"A string appended at the end of all files, use it to name simulations."))
((bool,noFiles,false,,"If true, no files will be generated (\\*.xml, \\*.spheres,...)"))
((Real,frictionAngleDegree,-1,,"Value of friction assigned just before the deviatoric loading"))
@@ -108,7 +108,7 @@
((Real,uniaxialEpsilonCurr,1,,"Current value of axial deformation during confined loading (is reference to strain[1])"))
((Real,fixedPoroCompaction,false,,"A special type of compaction with imposed final porosity :yref:`TriaxialCompressionEngine::fixedPorosity` (WARNING : can give unrealistic results!)"))
((Real,fixedPorosity,0,,"Value of porosity chosen by the user"))
- ((Real,maxStress,0,,"Max value of stress during the simulation (for post-processing)"))
+ ((Real,maxStress,0,,"Max absolute value of axial stress during the simulation (for post-processing)"))
((Real,sigma_iso,0,,"prescribed confining stress (see :yref:TriaxialCompressionEngine::isAxisymetric`)"))
((bool,isAxisymetric,false,,"if true, sigma_iso is assigned to sigma1, 2 and 3 (applies at each iteration and overrides user-set values of s1,2,3)"))
,
=== modified file 'pkg/dem/TriaxialStressController.cpp'
--- pkg/dem/TriaxialStressController.cpp 2014-10-15 06:44:01 +0000
+++ pkg/dem/TriaxialStressController.cpp 2014-11-18 19:46:27 +0000
@@ -75,7 +75,7 @@
}
}
-void TriaxialStressController::controlExternalStress(int wall, Vector3r resultantForce, State* p, Real wall_max_vel)
+void TriaxialStressController::controlExternalStress(int wall, Vector3r resultantForce, State* p, Real wall_max_vel) // controls walls such that Sum Forces from Sample on Wall = resultantForce
{
scene->forces.sync();
Real translation=normal[wall].dot(getForce(scene,wall_id[wall])-resultantForce);
@@ -168,31 +168,31 @@
if (!internalCompaction) {
Vector3r wallForce (0, goal2*width*depth, 0);
if (wall_bottom_activated) {
- if (stressMask & 2) controlExternalStress(wall_bottom, -wallForce, p_bottom, max_vel2);
+ if (stressMask & 2) controlExternalStress(wall_bottom, wallForce, p_bottom, max_vel2);
else p_bottom->vel[1] += (-normal[wall_bottom][1]*0.5*goal2*height -p_bottom->vel[1])*(1-strainDamping);
} else p_bottom->vel=Vector3r::Zero();
if (wall_top_activated) {
- if (stressMask & 2) controlExternalStress(wall_top, wallForce, p_top, max_vel2);
+ if (stressMask & 2) controlExternalStress(wall_top, -wallForce, p_top, max_vel2);
else p_top->vel[1] += (-normal[wall_top][1]*0.5*goal2*height -p_top->vel[1])*(1-strainDamping);
} else p_top->vel=Vector3r::Zero();
wallForce = Vector3r(goal1*height*depth, 0, 0);
if (wall_left_activated) {
- if (stressMask & 1) controlExternalStress(wall_left, -wallForce, p_left, max_vel1);
+ if (stressMask & 1) controlExternalStress(wall_left, wallForce, p_left, max_vel1);
else p_left->vel[0] += (-normal[wall_left][0]*0.5*goal1*width -p_left->vel[0])*(1-strainDamping);
} else p_left->vel=Vector3r::Zero();
if (wall_right_activated) {
- if (stressMask & 1) controlExternalStress(wall_right, wallForce, p_right, max_vel1);
+ if (stressMask & 1) controlExternalStress(wall_right, -wallForce, p_right, max_vel1);
else p_right->vel[0] += (-normal[wall_right][0]*0.5*goal1*width -p_right->vel[0])*(1-strainDamping);
} else p_right->vel=Vector3r::Zero();
wallForce = Vector3r(0, 0, goal3*height*width);
if (wall_back_activated) {
- if (stressMask & 4) controlExternalStress(wall_back, -wallForce, p_back, max_vel3);
+ if (stressMask & 4) controlExternalStress(wall_back, wallForce, p_back, max_vel3);
else p_back->vel[2] += (-normal[wall_back][2]*0.5*goal3*depth -p_back->vel[2])*(1-strainDamping);
} else p_back->vel=Vector3r::Zero();
if (wall_front_activated) {
- if (stressMask & 4) controlExternalStress(wall_front, wallForce, p_front, max_vel3);
+ if (stressMask & 4) controlExternalStress(wall_front, -wallForce, p_front, max_vel3);
else p_front->vel[2] += (-normal[wall_front][2]*0.5*goal3*depth -p_front->vel[2])*(1-strainDamping);
} else p_front->vel=Vector3r::Zero();
}
@@ -200,14 +200,14 @@
{
p_bottom->vel=Vector3r::Zero(); p_top->vel=Vector3r::Zero(); p_left->vel=Vector3r::Zero(); p_right->vel=Vector3r::Zero(); p_back->vel=Vector3r::Zero(); p_front->vel=Vector3r::Zero();
if (isARadiusControlIteration) {
- Real sigma_iso_ = bool(stressMask & 1)*goal1 + bool(stressMask & 2)*goal2 + bool(stressMask & 4)*goal3;
+ Real sigma_iso_ = std::abs( bool(stressMask & 1)*goal1 + bool(stressMask & 2)*goal2 + bool(stressMask & 4)*goal3 );
sigma_iso_ /= bool(stressMask & 1) + bool(stressMask & 2) + bool(stressMask & 4);
- if (sigma_iso_<=meanStress) maxMultiplier = finalMaxMultiplier;
+ if (sigma_iso_<=std::abs(meanStress)) maxMultiplier = finalMaxMultiplier;
if (meanStress==0) previousMultiplier = maxMultiplier;
else {
// previousMultiplier = 1+0.7*(sigma_iso-s)*(previousMultiplier-1.f)/(s-previousStress); // = (Dsigma/apparentModulus)*0.7
// previousMultiplier = std::max(2-maxMultiplier, std::min(previousMultiplier, maxMultiplier));
- previousMultiplier = 1.+(sigma_iso_-meanStress)/sigma_iso_*(maxMultiplier-1.); // = (Dsigma/apparentModulus)*0.7
+ previousMultiplier = 1.+(sigma_iso_-std::abs(meanStress))/sigma_iso_*(maxMultiplier-1.); // = (Dsigma/apparentModulus)*0.7
}
previousStress = meanStress;
//Real apparentModulus = (s-previousStress)/(previousMultiplier-1.f);
@@ -234,24 +234,24 @@
if (height0 == 0) height0 = height;
if (width0 == 0) width0 = width;
if (depth0 == 0) depth0 = depth;
- strain[0] = log(width0/width);
- strain[1] = log(height0/height);
- strain[2] = log(depth0/depth);
+ strain[0] = log(width/width0); // all strain values are positiv for extension
+ strain[1] = log(height/height0);
+ strain[2] = log(depth/depth0);
volumetricStrain=strain[0]+strain[1]+strain[2];
Real invXSurface = 1.f/(height*depth);
Real invYSurface = 1.f/(width*depth);
Real invZSurface = 1.f/(width*height);
- force[wall_bottom]=getForce(scene,wall_id[wall_bottom]); stress[wall_bottom]=force[wall_bottom]*invYSurface;
- force[wall_top]= getForce(scene,wall_id[wall_top]); stress[wall_top]=force[wall_top]*invYSurface;
+ force[wall_bottom]=getForce(scene,wall_id[wall_bottom]); stress[wall_bottom]=force[wall_bottom]*invYSurface; // all stress values are positiv for tension
+ force[wall_top]= getForce(scene,wall_id[wall_top]); stress[wall_top]=-force[wall_top]*invYSurface;
force[wall_left]= getForce(scene,wall_id[wall_left]); stress[wall_left]=force[wall_left]*invXSurface;
- force[wall_right]= getForce(scene,wall_id[wall_right]); stress[wall_right]= force[wall_right]*invXSurface;
- force[wall_front]= getForce(scene,wall_id[wall_front]); stress[wall_front]=force[wall_front]*invZSurface;
+ force[wall_right]= getForce(scene,wall_id[wall_right]); stress[wall_right]= -force[wall_right]*invXSurface;
+ force[wall_front]= getForce(scene,wall_id[wall_front]); stress[wall_front]=-force[wall_front]*invZSurface;
force[wall_back]= getForce(scene,wall_id[wall_back]); stress[wall_back]= force[wall_back]*invZSurface;
- for (int i=0; i<6; i++) meanStress-=stress[i].dot(normal[i]);
- meanStress/=6.;
+ for (int i=0; i<6; i++) meanStress+=stress[i].dot(pow(-1.0,i)*normal[i]); // normal[i] is always inwards
+ meanStress/=6.; // ( sXX(xLeft) + sXX(xRight) + sYY(yBottom) + sYY(yTop) + sZZ(zBack) + sZZ(zFront) ) / 6
}
void TriaxialStressController::controlInternalStress ( Real multiplier )
=== modified file 'pkg/dem/TriaxialStressController.hpp'
--- pkg/dem/TriaxialStressController.hpp 2014-10-15 06:44:01 +0000
+++ pkg/dem/TriaxialStressController.hpp 2014-11-18 19:46:27 +0000
@@ -75,7 +75,7 @@
Vector3r getStrainRate();
YADE_CLASS_BASE_DOC_ATTRS_DEPREC_INIT_CTOR_PY(TriaxialStressController,BoundaryController,
- "An engine maintaining constant stresses or constant strain rates on some boundaries of a parallepipedic packing. The stress/strain control is defined for each axis using :yref:`TriaxialStressController::stressMask` (a bitMask) and target values are defined by goal1,goal2, and goal3. sigmaIso has to be defined during growing phases."
+ "An engine maintaining constant stresses or constant strain rates on some boundaries of a parallepipedic packing. The stress/strain control is defined for each axis using :yref:`TriaxialStressController::stressMask` (a bitMask) and target values are defined by goal1,goal2, and goal3. All strains and stresses values are considered as positiv for extension."
"\n\n.. note::\n\t The algorithms used have been developed initialy for simulations reported in [Chareyre2002a]_ and [Chareyre2005]_. They have been ported to Yade in a second step and used in e.g. [Kozicki2008]_,[Scholtes2009b]_,[Jerier2010b]."
,
((unsigned int,stiffnessUpdateInterval,10,,"target strain rate (./s)"))
@@ -103,9 +103,9 @@
((Real,width0,0,,"Reference size for strain definition. See :yref:`TriaxialStressController::width`"))
((Real,depth0,0,,"Reference size for strain definition. See :yref:`TriaxialStressController::depth`"))
- ((Real,goal1,0,,"prescribed stress/strain rate on axis 1, as defined by :yref:`TriaxialStressController::stressMask` (see also :yref:`TriaxialStressController::isAxisymetric`)"))
- ((Real,goal2,0,,"prescribed stress/strain rate on axis 2, as defined by :yref:`TriaxialStressController::stressMask` (see also :yref:`TriaxialStressController::isAxisymetric`)"))
- ((Real,goal3,0,,"prescribed stress/strain rate on axis 3, as defined by :yref:`TriaxialStressController::stressMask` (see also :yref:`TriaxialStressController::isAxisymetric`)"))
+ ((Real,goal1,0,,"prescribed stress/strain rate (> 0 in extension) on axis 1, as defined by :yref:`TriaxialStressController::stressMask`"))
+ ((Real,goal2,0,,"prescribed stress/strain rate (> 0 in extension) on axis 2, as defined by :yref:`TriaxialStressController::stressMask`"))
+ ((Real,goal3,0,,"prescribed stress/strain rate (> 0 in extension) on axis 3, as defined by :yref:`TriaxialStressController::stressMask`"))
((unsigned int,stressMask,7,,"Bitmask determining if the components of :yref:`TriaxialStressController::goal` are stress (1) or strain (0). 0 for none, 7 for all, 1 for sigma1, etc."))
((Real,maxMultiplier,1.001,,"max multiplier of diameters during internal compaction (initial fast increase - :yref:`TriaxialStressController::finalMaxMultiplier` is used in a second stage)"))
((Real,finalMaxMultiplier,1.00001,,"max multiplier of diameters during internal compaction (secondary precise adjustment - :yref:`TriaxialStressController::maxMultiplier` is used in the initial stage)"))
@@ -113,7 +113,7 @@
((Real,previousStress,0,Attr::readonly,"|yupdate|"))
((Real,previousMultiplier,1,Attr::readonly,"|yupdate|"))
((bool,internalCompaction,true,,"Switch between 'external' (walls) and 'internal' (growth of particles) compaction."))
- ((Real,meanStress,0,Attr::readonly,"Mean stress in the packing. |yupdate|"))
+ ((Real,meanStress,0,Attr::readonly,"Mean stress (> 0 in extension) in the packing. |yupdate|"))
((Real,volumetricStrain,0,Attr::readonly,"Volumetric strain (see :yref:`TriaxialStressController::strain`).|yupdate|"))
((Real,externalWork,0,Attr::readonly,"Energy provided by boundaries.|yupdate|"))
((bool,updatePorosity,false,,"If true porosity calculation will be updated once (will automatically reset to false after one calculation step). Can be used, when volume of particles changes during the simulation (e.g. when particles are erased or when clumps are created)."))
@@ -127,6 +127,7 @@
/* extra initializers */
,
/* constructor */
+ strain = Vector3r::Zero();
first = true;
stiffness.resize(6);
previousTranslation.assign(Vector3r::Zero());
@@ -139,7 +140,7 @@
normal[wall_back].z()=1;
porosity=1;
,
- .def_readonly("strain",&TriaxialStressController::strain,"Current strain in a vector (exx,eyy,ezz). The values reflect true (logarithmic) strain.")
+ .def_readonly("strain",&TriaxialStressController::strain,"Current strain in a vector (exx,eyy,ezz). The values reflect true (logarithmic) strain and are > 0 in extension.")
.def_readonly("strainRate",&TriaxialStressController::getStrainRate,"Current strain rate in a vector d/dt(exx,eyy,ezz).")
.def_readonly("porosity",&TriaxialStressController::porosity,"Porosity of the packing.")
.def_readonly("boxVolume",&TriaxialStressController::boxVolume,"Total packing volume.")
@@ -148,7 +149,7 @@
.def_readonly("max_vel1",&TriaxialStressController::max_vel1,"see :yref:`TriaxialStressController::max_vel` |ycomp|")
.def_readonly("max_vel2",&TriaxialStressController::max_vel2,"see :yref:`TriaxialStressController::max_vel` |ycomp|")
.def_readonly("max_vel3",&TriaxialStressController::max_vel3,"see :yref:`TriaxialStressController::max_vel` |ycomp|")
- .def("stress",&TriaxialStressController::getStress,(boost::python::arg("id")),"Return the mean stress vector acting on boundary 'id', with 'id' between 0 and 5.")
+ .def("stress",&TriaxialStressController::getStress,(boost::python::arg("id")),"Returns the mean stress (> 0 in tension) existing in the sample, as measured along boundary 'id', with 'id' between 0 and 5.")
)
DECLARE_LOGGER;
};
=== modified file 'pkg/dem/TriaxialTest.hpp'
--- pkg/dem/TriaxialTest.hpp 2014-10-15 06:44:01 +0000
+++ pkg/dem/TriaxialTest.hpp 2014-11-18 19:46:27 +0000
@@ -99,8 +99,8 @@
((Real,dampingMomentum,0.2,,"Coefficient of Cundal-Non-Viscous damping (applied on on the 3 components of torques)"))
((Real,StabilityCriterion,0.01,,"Value of unbalanced force for which the system is considered stable. Used in conditionals to switch between loading stages."))
((Real,wallOversizeFactor,1.3,,"Make boundaries larger than the packing to make sure spheres don't go out during deformation."))
- ((Real,sigmaIsoCompaction,50000,,"Confining stress during isotropic compaction."))
- ((Real,sigmaLateralConfinement,50000,,"Lateral stress during triaxial loading. An isotropic unloading is performed if the value is not equal to :yref:`TriaxialTest::sigmaIsoCompaction`."))
+ ((Real,sigmaIsoCompaction,-50000,,"Confining stress during isotropic compaction (< 0 for real - compressive - compaction)."))
+ ((Real,sigmaLateralConfinement,-50000,,"Lateral stress during triaxial loading (< 0 for classical compressive cases). An isotropic unloading is performed if the value is not equal to :yref:`TriaxialTest::sigmaIsoCompaction`."))
((int,timeStepUpdateInterval,50,,"interval for :yref:`GlobalStiffnessTimeStepper`"))
((int,wallStiffnessUpdateInterval,10,,"interval for updating the stiffness of sample/boundaries contacts"))
=== modified file 'py/pack/pack.py'
--- py/pack/pack.py 2014-01-06 14:10:07 +0000
+++ py/pack/pack.py 2014-11-18 19:46:27 +0000
@@ -469,7 +469,7 @@
# upperCorner is just size ratio, if radiusMean is specified
upperCorner=fullDim,
## no need to touch any the following
- noFiles=True,lowerCorner=[0,0,0],sigmaIsoCompaction=1e7,sigmaLateralConfinement=1e5,compactionFrictionDeg=1,StabilityCriterion=.05,strainRate=.2,thickness=-1,maxWallVelocity=.1,wallOversizeFactor=1.5,autoUnload=True,autoCompressionActivation=False).load()
+ noFiles=True,lowerCorner=[0,0,0],sigmaIsoCompaction=-1e7,sigmaLateralConfinement=-1e5,compactionFrictionDeg=1,StabilityCriterion=.05,strainRate=.2,thickness=-1,maxWallVelocity=.1,wallOversizeFactor=1.5,autoUnload=True,autoCompressionActivation=False).load()
while ( numpy.isnan(utils.unbalancedForce()) or utils.unbalancedForce()>0.005 ) :
O.run(100,True)
sp=SpherePack(); sp.fromSimulation()
=== modified file 'scripts/checks-and-tests/checks/DEM-PFV-check.py'
--- scripts/checks-and-tests/checks/DEM-PFV-check.py 2014-08-07 13:38:16 +0000
+++ scripts/checks-and-tests/checks/DEM-PFV-check.py 2014-11-18 19:46:27 +0000
@@ -51,12 +51,12 @@
newton
]
- triax.goal1=triax.goal2=triax.goal3=10000
+ triax.goal1=triax.goal2=triax.goal3=-10000
while 1:
O.run(200, True)
unb=unbalancedForce()
- if unb<0.01 and abs(10000-triax.meanStress)/10000<0.01: break
+ if unb<0.01 and abs(-10000-triax.meanStress)/10000<0.01: break
setContactFriction(radians(finalFricDegree))
@@ -68,11 +68,11 @@
triax.internalCompaction=False
triax.wall_bottom_activated=False
- triax.goal2=11000; O.run(2000,1)
- triax.goal2=10000; O.run(2000,1)
- triax.goal2=11000; O.run(2000,1)
+ triax.goal2=-11000; O.run(2000,1)
+ triax.goal2=-10000; O.run(2000,1)
+ triax.goal2=-11000; O.run(2000,1)
e22=triax.strain[1]
- triax.goal2=10000; O.run(2000,1)
+ triax.goal2=-10000; O.run(2000,1)
e22=e22-triax.strain[1]
modulus = 1000./abs(e22)
@@ -119,7 +119,7 @@
zeroTime=O.time
zeroe22 = triax.strain[1]
- triax.goal2=11000
+ triax.goal2=-11000
O.timingEnabled=1
from yade import timing
@@ -132,7 +132,7 @@
errors+=1
print "The difference is more, than the critical tolerance!"
- target=0.00258113045083
+ target=-0.00258113045083
if abs((triax.strain[1]-zeroe22-target)/target)>toleranceWarning:
print "DEM-PFV: difference in final deformation",triax.strain[1]-zeroe22," vs. target ",target
if (abs((triax.strain[1]-zeroe22-target)/target)>toleranceCritical):
=== modified file 'scripts/checks-and-tests/checks/checkTestTriax.py'
--- scripts/checks-and-tests/checks/checkTestTriax.py 2013-08-29 10:30:31 +0000
+++ scripts/checks-and-tests/checks/checkTestTriax.py 2014-11-18 19:46:27 +0000
@@ -15,13 +15,13 @@
O.load("checkTest.yade")
O.run(2020,True)
if interactive:
- print O.engines[4].stress(0)[1],O.engines[4].stress(1)[1], O.engines[4].stress(2)[0], O.engines[4].stress(3)[0], O.engines[4].stress(4)[2], O.engines[4].stress(5)[2]
+ print O.engines[4].stress(0)[0],O.engines[4].stress(1)[0], O.engines[4].stress(2)[1], O.engines[4].stress(3)[1], O.engines[4].stress(4)[2], O.engines[4].stress(5)[2]
os.system('gnuplot -e "plot \'./WallStresses\' using 1:3; replot'+checksPath+'\'/data/WallStressesCheckTest\' using 1:3; replot '+checksPath+'\'/data/WallStresses\' using 1:4; replot '+checksPath+'\'/data/WallStressesCheckTest\' using 1:4; pause -1"')
-if abs((O.engines[4].stress(3)[1]-107157.2)/107157.2)>tolerance :
+if abs((O.engines[4].stress(3)[1]+107157.2)/107157.2)>tolerance :
print "Triaxial checkTest: difference on peak stress"
errors+=1
-if abs((O.engines[4].stress(1)[0]-50058.7)/50058.7)>tolerance :
+if abs((O.engines[4].stress(1)[0]+50058.7)/50058.7)>tolerance :
print "Triaxial checkTest: difference on confining stress"
errors+=1
Follow ups
