calculate external work in Uniaxial Strainer

Hello,

You tried to modify the source code, that's it ? Did you also modify UniaxialStrainer.cpp in order to update this new variable "dAX" ?

As a side note, did you consider to use TriaxialStressController ? In spite of its name, it should be possible to use it for uniaxial loadings. And TriaxialStressController provides already access to external work:
https://yade-dem.org/doc/yade.wrapper.html#yade.wrapper.TriaxialStressController.externalWork

Jerome

Hi Jabrane,

you came to c++ variables and scopes [2].

to get correct behavior, one option is to change the declaration of dAX
[1]. Currently (Real dAX=....) it creates a new variable, unrelated to
strainer.dAX.
Just delete Real (i.e. assign directly dAX=...) and then dAX refers to the
class member strainer.dAX.

Another option is leave it as it is and add
this->dAX = dAX;
somewhere after dAX is computed.

cheers
Jan

[1]
https://github.com/yade/trunk/blob/master/pkg/dem/UniaxialStrainer.cpp#L95
[2]
http://stackoverflow.com/questions/30494407/declaring-a-local-variable-within-class-scope-with-same-name-as-a-class-attribut

2016-06-14 17:37 GMT+02:00 Yor1 <email address hidden>:

> New question #295301 on Yade:
> https://answers.launchpad.net/yade/+question/295301
>
> Hello Yade users !
>
> I try to calculate the external work on a sample in tensile test.
> In fact, i calculate the external work with this simple formulate
> "Wext_incremental=AppliedForce*displacement_incremental".
> I want to integrate this formulate in the sources UniaxialStrainer.cpp
> In uniaxialStrainer.cpp the incremental displacement is calculated and
> named "dAX".
> In order to verify the behavior of incremental displacement, i define
> "dAX" in uniaxialStrainer.hpp like this
> ((Real,dAX,0,,"Current incremental displacement")) in the line 43 of
> UniaxialStrainer.hpp
>
> My problem is that in the simulation of tensile test , dAX=0 during the
> simulation.
> I call dAX in the python script strainer.dAX
> I don't understand where is the problem.
>
> https://github.com/yade/trunk/blob/master/pkg/dem/UniaxialStrainer.hpp
> https://github.com/yade/trunk/blob/master/pkg/dem/UniaxialStrainer.cpp
>
> Best regards.
> Jabrane.
>
>
> --
> You received this question notification because your team yade-users is
> an answer contact for Yade.
>
> _______________________________________________
> Mailing list: https://launchpad.net/~yade-users
> Post to : <email address hidden>
> Unsubscribe : https://launchpad.net/~yade-users
> More help : https://help.launchpad.net/ListHelp
>

Hello Jerome and Jan

Thank you of your response.

Jerome, I simulate the tensile test with TriaxialStressController using a positive strainRate but i doesn't work.
I get sigma=sigma2=sigma3=0 and that is anormal.

Jan, your solution work but i can't have a definitive opinion until i calculate Wext.
But i have a question how do you know that dAX is a local variable ?

Best regards.
Jabrane.

Hi Jabrane,

But i have a question how do you know that dAX is a local variable ?

it just comes from how C++ works / is defined, nothing else.

"Real dAX" in cpp file inside a function simply declares a local variable.
No matter if the class has a member with the same name.

cheers
Jan

2016-06-14 19:37 GMT+02:00 Yor1 <email address hidden>:

> Question #295301 on Yade changed:
> https://answers.launchpad.net/yade/+question/295301
>
> Yor1 posted a new comment:
> Hello Jerome and Jan
>
> Thank you of your response.
>
> Jerome, I simulate the tensile test with TriaxialStressController using a
> positive strainRate but i doesn't work.
> I get sigma=sigma2=sigma3=0 and that is anormal.
>
> Jan, your solution work but i can't have a definitive opinion until i
> calculate Wext.
> But i have a question how do you know that dAX is a local variable ?
>
> Best regards.
> Jabrane.
>
> --
> You received this question notification because your team yade-users is
> an answer contact for Yade.
>
> _______________________________________________
> Mailing list: https://launchpad.net/~yade-users
> Post to : <email address hidden>
> Unsubscribe : https://launchpad.net/~yade-users
> More help : https://help.launchpad.net/ListHelp
>

Hello Jan,

The behavior of dAX is related to the line in which i put "this->dAX = dAX"
In fact i put "this->dAX = dAX" just after https://github.com/yade/trunk/blob/master/pkg/dem/UniaxialStrainer.cpp#L95
and just after https://github.com/yade/trunk/blob/master/pkg/dem/UniaxialStrainer.cpp#L122 and i get different value of dAX.

The question is where do I put "this->dAX = dAX" to get the real behavior?

Best regards.
Jabrane.

Hi Jabrane,
you should put "this->dAX = dAX" to a place such that dAX is not changed
afterwards. e.g. it is changed at
https://github.com/yade/trunk/blob/master/pkg/dem/UniaxialStrainer.cpp#L106
So after line 121 (in the original file) it should be the correct position
Jan

2016-06-14 19:47 GMT+02:00 Jan Stránský <
<email address hidden>>:

> Question #295301 on Yade changed:
> https://answers.launchpad.net/yade/+question/295301
>
> Jan Stránský proposed the following answer:
> Hi Jabrane,
>
> But i have a question how do you know that dAX is a local variable ?
>
>
> it just comes from how C++ works / is defined, nothing else.
>
> "Real dAX" in cpp file inside a function simply declares a local variable.
> No matter if the class has a member with the same name.
>
> cheers
> Jan
>
>
> 2016-06-14 19:37 GMT+02:00 Yor1 <email address hidden>:
>
> > Question #295301 on Yade changed:
> > https://answers.launchpad.net/yade/+question/295301
> >
> > Yor1 posted a new comment:
> > Hello Jerome and Jan
> >
> > Thank you of your response.
> >
> > Jerome, I simulate the tensile test with TriaxialStressController using a
> > positive strainRate but i doesn't work.
> > I get sigma=sigma2=sigma3=0 and that is anormal.
> >
> > Jan, your solution work but i can't have a definitive opinion until i
> > calculate Wext.
> > But i have a question how do you know that dAX is a local variable ?
> >
> > Best regards.
> > Jabrane.
> >
> > --
> > You received this question notification because your team yade-users is
> > an answer contact for Yade.
> >
> > _______________________________________________
> > Mailing list: https://launchpad.net/~yade-users
> > Post to : <email address hidden>
> > Unsubscribe : https://launchpad.net/~yade-users
> > More help : https://help.launchpad.net/ListHelp
> >
>
> --
> You received this question notification because your team yade-users is
> an answer contact for Yade.
>
> _______________________________________________
> Mailing list: https://launchpad.net/~yade-users
> Post to : <email address hidden>
> Unsubscribe : https://launchpad.net/~yade-users
> More help : https://help.launchpad.net/ListHelp
>

Hi Jan,

Thank you for the answer.
Your alternative work but it doesn't resolve my problem because the external work obtained is underestimated (i compare it with the elastic energy).
I tried a second alternative which based on the difference of the total displacement between iteration "i" and "i+1".
The problem is that in UniaxialStrainer.cpp i can't do this difference because we can't control the
I tried these lines of code but it doesn't work:

Real displacement; Real S;
S=axialLength-originalLength;
displacement=(axialLength-originalLength)-S;

I obtained displacement=0. May be this solution is stupid but i don't know how to do the get the difference in
UniaxialStrainer.cpp

Best regards.
Jabrane.

>
>
> Your alternative work but it doesn't resolve my problem because the
> external work obtained is underestimated (i compare it with the elastic
> energy).
>

Do you have some values? for example
https://github.com/yade/trunk/blob/master/pkg/dem/UniaxialStrainer.cpp#L106
changes the value to half, but if you compute work from it, you should
multiply it by 2 again.. Could you also post a script you use?

> I tried a second alternative which based on the difference of the total
> displacement between iteration "i" and "i+1".
> The problem is that in UniaxialStrainer.cpp i can't do this difference
> because we can't control the
> I tried these lines of code but it doesn't work:
>
> Real displacement; Real S;
> S=axialLength-originalLength;
> displacement=(axialLength-originalLength)-S;

> I obtained displacement=0. May be this solution is stupid but i don't know
> how to do the get the difference in
> UniaxialStrainer.cpp
>

of course :-) displacement=(axialLength-originalLength)-S
= (axialLength-originalLength) - (axialLength-originalLength) = 0 for any
values of axisLength and originalLength

the solution is to store the previous value inside the class

cheers
Jan

Hi Jan,

This is the script which i use:

from yade import ymport, utils , plot
import math

PACKING='X1Y2Z1_2k'
OUT=PACKING+'_tensionTest_r0002_energy'

DAMP=0.4 s
saveData=100
iterMax=40000
saveVTK=1000

strainRate=0.002

intR=1.5028
DENS=4000
YOUNG=65e9
FRICT=10
ALPHA=0.4
TENS=8e6
COH=160e6

def sphereMat(): return JCFpmMat(type=1,density=DENS,young=YOUNG,poisson = ALPHA,frictionAngle=radians(FRICT),tensileStrength=TENS,cohesion=COH)

O.bodies.append(ymport.text(PACKING+'.spheres',scale=1.,shift=Vector3(0,0,0),material=sphereMat))

dim=utils.aabbExtrema()
xinf=dim[0][0]
xsup=dim[1][0]
X=xsup-xinf
yinf=dim[0][1]
ysup=dim[1][1]
Y=ysup-yinf
zinf=dim[0][2]
zsup=dim[1][2]
Z=zsup-zinf

bb=utils.uniaxialTestFeatures()
negIds,posIds,axis,crossSectionArea=bb['negIds'],bb['posIds'],bb['axis'],bb['area']

O.engines=[
 ForceResetter(),
        InsertionSortCollider([Bo1_Box_Aabb(),Bo1_Sphere_Aabb(aabbEnlargeFactor=intR,label='Saabb')]),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=intR,label='SSgeom'),Ig2_Box_Sphere_ScGeom()],
  [Ip2_JCFpmMat_JCFpmMat_JCFpmPhys(cohesiveTresholdIteration=1,label='interactionPhys')],
  [Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM(label='interactionLaw')]
 ),
 UniaxialStrainer(strainRate=strainRate,axis=axis,asymmetry=0,posIds=posIds,negIds=negIds,crossSectionArea=crossSectionArea,blockDisplacements=1,blockRotations=1,setSpeeds=0,stopStrain=0.1,label='strainer'),
 GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=10,timestepSafetyCoefficient=0.4, defaultDt=0.1*utils.PWaveTimeStep()),
 NewtonIntegrator(damping=DAMP,label='newton'),
 PyRunner(iterPeriod=int(saveData),initRun=True,command='recorder()',label='data'),
        #VTKRecorder(iterPeriod=int(saveVTK),initRun=True,fileName=OUT+'-',recorders=['spheres','jcfpm','cracks'],Key=OUT,label='vtk')
]

O.step();

SSgeom.interactionDetectionFactor=-1.
Saabb.aabbEnlargeFactor=-1.

layerSize=0.2
for o in O.bodies:
  if isinstance(o.shape,Sphere):
    if ( o.state.pos[axis]<(dim[0][axis]+layerSize*(dim[1][axis]-dim[0][axis])) ) or ( o.state.pos[axis]>(dim[1][axis]-layerSize*(dim[1][axis]-dim[0][axis])) ) :
      o.shape.color=(1,1,1)

O.run(iterMax)

Hi Jabrane,
thanks for the script. However, I was mainly interested in how you compute
the elastic energy and the work :-)
Jan

2016-06-16 9:28 GMT+02:00 Yor1 <email address hidden>:

> Question #295301 on Yade changed:
> https://answers.launchpad.net/yade/+question/295301
>
> Yor1 posted a new comment:
> Hi Jan,
>
> This is the script which i use:
>
> from yade import ymport, utils , plot
> import math
>
> PACKING='X1Y2Z1_2k'
> OUT=PACKING+'_tensionTest_r0002_energy'
>
> DAMP=0.4 s
> saveData=100
> iterMax=40000
> saveVTK=1000
>
> strainRate=0.002
>
> intR=1.5028
> DENS=4000
> YOUNG=65e9
> FRICT=10
> ALPHA=0.4
> TENS=8e6
> COH=160e6
>
> def sphereMat(): return JCFpmMat(type=1,density=DENS,young=YOUNG,poisson
> = ALPHA,frictionAngle=radians(FRICT),tensileStrength=TENS,cohesion=COH)
>
>
> O.bodies.append(ymport.text(PACKING+'.spheres',scale=1.,shift=Vector3(0,0,0),material=sphereMat))
>
> dim=utils.aabbExtrema()
> xinf=dim[0][0]
> xsup=dim[1][0]
> X=xsup-xinf
> yinf=dim[0][1]
> ysup=dim[1][1]
> Y=ysup-yinf
> zinf=dim[0][2]
> zsup=dim[1][2]
> Z=zsup-zinf
>
> bb=utils.uniaxialTestFeatures()
>
> negIds,posIds,axis,crossSectionArea=bb['negIds'],bb['posIds'],bb['axis'],bb['area']
>
> O.engines=[
> ForceResetter(),
>
> InsertionSortCollider([Bo1_Box_Aabb(),Bo1_Sphere_Aabb(aabbEnlargeFactor=intR,label='Saabb')]),
> InteractionLoop(
>
> [Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=intR,label='SSgeom'),Ig2_Box_Sphere_ScGeom()],
>
> [Ip2_JCFpmMat_JCFpmMat_JCFpmPhys(cohesiveTresholdIteration=1,label='interactionPhys')],
>
> [Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM(label='interactionLaw')]
> ),
>
> UniaxialStrainer(strainRate=strainRate,axis=axis,asymmetry=0,posIds=posIds,negIds=negIds,crossSectionArea=crossSectionArea,blockDisplacements=1,blockRotations=1,setSpeeds=0,stopStrain=0.1,label='strainer'),
>
> GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=10,timestepSafetyCoefficient=0.4,
> defaultDt=0.1*utils.PWaveTimeStep()),
> NewtonIntegrator(damping=DAMP,label='newton'),
>
> PyRunner(iterPeriod=int(saveData),initRun=True,command='recorder()',label='data'),
>
> #VTKRecorder(iterPeriod=int(saveVTK),initRun=True,fileName=OUT+'-',recorders=['spheres','jcfpm','cracks'],Key=OUT,label='vtk')
> ]
>
> O.step();
>
> SSgeom.interactionDetectionFactor=-1.
> Saabb.aabbEnlargeFactor=-1.
>
> layerSize=0.2
> for o in O.bodies:
> if isinstance(o.shape,Sphere):
> if (
> o.state.pos[axis]<(dim[0][axis]+layerSize*(dim[1][axis]-dim[0][axis])) ) or
> ( o.state.pos[axis]>(dim[1][axis]-layerSize*(dim[1][axis]-dim[0][axis])) ) :
> o.shape.color=(1,1,1)
>
> O.run(iterMax)
>
> --
> You received this question notification because your team yade-users is
> an answer contact for Yade.
>
> _______________________________________________
> Mailing list: https://launchpad.net/~yade-users
> Post to : <email address hidden>
> Unsubscribe : https://launchpad.net/~yade-users
> More help : https://help.launchpad.net/ListHelp
>

Hi Jan,

I compute the elastic energy in the recorder (in the script) with these lines:

for i in O.interactions:
      if not i.isReal : continue
      E+=0.5*(i.phys.normalForce.squaredNorm()/i.phys.kn + i.phys.shearForce.squaredNorm()/i.phys.ks)

With these lines i guarantee that i calculate the elastic energy in all contacts (contact between spheres, contact between sphere and wall, etc ...)

For the external work i try to compute it in the sources in https://github.com/yade/trunk/blob/master/pkg/dem/UniaxialStrainer.cpp#L106

if(asymmetry==0){
dAX*=.5;displacement=2*dAX;externalWork+=displacement(sumPosForces+sumNegForces)/2;externalWork+=externalWork;}

Best regards.
Jabrane.

Hi Jabrane,

dAX*=.5;
> displacement=2*dAX;
> externalWork += displacement(sumPosForces+sumNegForces)/2;
> externalWork += externalWork;

a few points:
1) displacement(sumPosForces+sumNegForces)/2 should give compilation error.
I consider there should be multiplication displacement*(
sumPosForces+sumNegForces)..

2) also from your previous post with S=axialLength-originalLength;
displacement=(axialLength-originalLength)-S;, think twice what your
commands do

externalWork += displacement(sumPosForces+sumNegForces)/2; // you add
dspl*force to externalWork, so far OK
externalWork += externalWork; // this is equal to externalWork =
2*externalWork.. I don's see any reason for this

3) sumPosForces and sumNegForces are updated at the very end of action
function by computeAxialForce, and only once per stressUpdateInterval, so
to compute it more accurately, you should put your computation after this
update and use some low value (1 to get really exact results) for
stressUpdateInterval.

please consider to update you computations and let us know about new results
cheers
Jan

2016-06-16 14:13 GMT+02:00 Yor1 <email address hidden>:

> Question #295301 on Yade changed:
> https://answers.launchpad.net/yade/+question/295301
>
> Yor1 posted a new comment:
> Hi Jan,
>
> I compute the elastic energy in the recorder (in the script) with these
> lines:
>
> for i in O.interactions:
> if not i.isReal : continue
> E+=0.5*(i.phys.normalForce.squaredNorm()/i.phys.kn +
> i.phys.shearForce.squaredNorm()/i.phys.ks)
>
> With these lines i guarantee that i calculate the elastic energy in all
> contacts (contact between spheres, contact between sphere and wall, etc
> ...)
>
> For the external work i try to compute it in the sources in
> https://github.com/yade/trunk/blob/master/pkg/dem/UniaxialStrainer.cpp#L106
>
> if(asymmetry==0){
>
> dAX*=.5;displacement=2*dAX;externalWork+=displacement(sumPosForces+sumNegForces)/2;externalWork+=externalWork;}
>
> Best regards.
> Jabrane.
>
> --
> You received this question notification because your team yade-users is
> an answer contact for Yade.
>
> _______________________________________________
> Mailing list: https://launchpad.net/~yade-users
> Post to : <email address hidden>
> Unsubscribe : https://launchpad.net/~yade-users
> More help : https://help.launchpad.net/ListHelp
>

Hi Jan,

I calculate the displacement and the external work with these lines:

if(asymmetry==0){ dAX*=.5;displacement=2*max(negIds.size(),posIds.size())*dAX;}
externalWork += displacement*(sumPosForces+sumNegForces)/2;

The eslastic energy and the external work are equal in the pre-pic stage and that is reasonable for me.
Now i have to modify the calculation of the damping energy. In fact with NewtonIntegrator, the damped energy is calculated in all the particles. But in my case, i have to calculate the damping energy in the particles that don't belong to posIds and negIds.

Best regards
Jabrane

>
> displacement=2*max(negIds.size(),posIds.size())*dAX;

I don't like that the work depends on number of particles, it should not..
Jan

2016-06-16 18:37 GMT+02:00 Yor1 <email address hidden>:

> Question #295301 on Yade changed:
> https://answers.launchpad.net/yade/+question/295301
>
> Yor1 posted a new comment:
> Hi Jan,
>
> I calculate the displacement and the external work with these lines:
>
> if(asymmetry==0){
> dAX*=.5;displacement=2*max(negIds.size(),posIds.size())*dAX;}
> externalWork += displacement*(sumPosForces+sumNegForces)/2;
>
>
> The eslastic energy and the external work are equal in the pre-pic stage
> and that is reasonable for me.
> Now i have to modify the calculation of the damping energy. In fact with
> NewtonIntegrator, the damped energy is calculated in all the particles. But
> in my case, i have to calculate the damping energy in the particles that
> don't belong to posIds and negIds.
>
> Best regards
> Jabrane
>
> --
> You received this question notification because your team yade-users is
> an answer contact for Yade.
>
> _______________________________________________
> Mailing list: https://launchpad.net/~yade-users
> Post to : <email address hidden>
> Unsubscribe : https://launchpad.net/~yade-users
> More help : https://help.launchpad.net/ListHelp
>

I change the manner that i calculate the displacement

if(asymmetry==0){
   dAX*=.5;
   for(size_t i=0; i<max(posIds.size(),negIds.size()); i++){
    displacement+=2*dAX;
   }
 }

I'm inspiring from the calculation of the displacement for "asymmetry!=1" and "assymetry!=-1".
For example, the displacement in the case of "asymmetry!=1" is calculated like this:

if(asymmetry!=1){
  for(size_t i=0; i<negIds.size(); i++){
   negCoords[i]-=dAX;
   axisVel(negIds[i]) = -dAX/scene->dt; // update current position

  }

 }

Jabrane.

>
>
>
> if(asymmetry!=1){
> for(size_t i=0; i<negIds.size(); i++){
> negCoords[i]-=dAX;
> axisVel(negIds[i]) = -dAX/scene->dt; // update
> current position
>
> }
>
> }
>
>
there is no computation of displacement here.. just each **coordinate**
(from 0 to negIds.size()) is updated by displacement dAX, so .. Really try
to understand what the code does before being inspired, otherwise it might
be "dangerous" :-)

dAX (or maximum 2*dAX) is the actual change in length of the specimen. If
you compute external work like this (multiplying it by number of fixed
particles) and it equals the elastic energy, there must be a problem
somewhere..

Jan

Hello Jan !

I modified the computation of the displacement and the external work by these formula:

displacement += strainRate*originalLength*scene->dt;
externalWork += displacement*(sumPosForces+sumNegForces)/2;

And it works and that resolves my problem.

Jabrane.