Triaxial test on dry granular material

2) triax.stress is a function [1], so as the error says, it has no __getitem__ (you cannot do triax.stress[0])
Jan

[1] https://yade-dem.org/doc/yade.wrapper.html#yade.wrapper.TriaxialStressController.stress

Hi

Should not the stressMask be equal to 0, for strain controlled?

Cheers
Amiya

To expand on Jan's answer:

Remember to read the documentation for the functions and controllers you plan to use. For example, stressMask [1]. You are telling the walls to reach stress goals with stressMask=7. But then you assign what appears to be a strain_rate to these goals. If I were you, I would spend 5 minutes just reading through all the functionalities of TriaxialStressController.

It's probably also worth your time to work through this triaxial example [2], and modify it as necessary.

[1]https://yade-dem.org/doc/yade.wrapper.html#yade.wrapper.TriaxialStressController.stressMask
[2]https://github.com/yade/trunk/blob/master/examples/triax-tutorial/script-session1.py

Thanks all.

@Amiya & Robert : sorry, the stressMask=7 was a typo. I have read about it in detail but just missed it there ;)

@Jan & Robert: I read the code in [1] and tried to rectify my script. triax is used there in a similar way as here. The problem still persists. Moreover, the simulation now goes blank after some time. The error obtained is "cannot convert float NaN to integer".

Any ideas to correct this now?

Here is the modified script.

from yade import pack,plot
import math

num_spheres=10000 # number of spheres
young=3e9
compFricDegree = 30 # initial contact friction during the confining phase
finalFricDegree = 30 # contact friction during the deviatoric loading
mn,mx=Vector3(0,0,0),Vector3(0.1,0.1,0.1) # corners of the initial packing

O.materials.append(FrictMat(young=young,poisson=0.3,frictionAngle=radians(compFricDegree),density=2600,label='spheres'))
O.materials.append(FrictMat(young=young,poisson=0.5,frictionAngle=0,density=0,label='walls'))
walls=aabbWalls([mn,mx],thickness=0,material='walls')
wallIds=O.bodies.append(walls)

sp=pack.SpherePack()
sp.makeCloud(mn,mx,-1,0.5,num_spheres,False, 0.65,seed=1) #"seed" make the "random" generation always the same
sp.toSimulation(material='spheres')

triax=TriaxialStressController(
 maxMultiplier=1.+2e4/young, # spheres growing factor (fast growth)
 finalMaxMultiplier=1.+2e3/young, # spheres growing factor (slow growth)
 thickness = 0,
 stressMask = 0,
 max_vel = 0.005,
 internalCompaction=True, # If true the confining pressure is generated by growing particles
)

newton=NewtonIntegrator(gravity=(0,0,-9.81),damping=0.3)

O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
  [Ip2_FrictMat_FrictMat_FrictPhys()],
  [Law2_ScGeom_FrictPhys_CundallStrack()],label="iloop"
 ),

 GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),PyRunner(iterPeriod=100,command="addPlotData()"),
 triax,
 newton
]

#loading only along Z-direction (1D consolidation test)
strain_rate=-0.0005
triax.goal1=triax.goal2=triax.goal3=strain_rate #---> Is this correct

O.dt=0.5*PWaveTimeStep()
def addPlotData():
  plot.addData(Sz=-triax.stress[triax.wall_right_id][0], Dz=-triax.strain[0])

plot.plots= {'Dz':('Sz')}
plot.plot()

[1] https://github.com/yade/trunk/blob/master/examples/triax-tutorial/script-session1.py

Correction:

For strain loading in one direction it should be ( for example):

strain_rate=-0.004
triax.goal1=strain_rate
triax.goal2=triax.goal3=0

But the main issue ( as described above) still remains

Hello,

> The problem still persists

again, triax.stress is a function and you cannot do triax.stress[...], first you have to call it with (...), so I guess
Sz=-triax.stress(triax.wall_right_id)[0]
could work (actually it is like that in your example script [1])

> Moreover, the simulation now goes blank after some time. The error obtained is "cannot convert float NaN to integer".

please in this situation always attach full error message and preferably as specific information as possible (e.g. instead of "after some time" something like "at iteration 1234" etc.)

cheers
Jan

Thanks for the help Jan.
I am trying the same: Sz=-triax.stress(triax.wall_right_id)[0] but it is not working with my script. I have been trying all possibilities to correct it but the same result everytime.
I am not able to figure out another way to do the task.

The following is the full message:
TypeError Traceback (most recent call last)
/usr/bin/yade in <module>()

/usr/bin/yade in addPlotData()
     47 O.dt=0.5*PWaveTimeStep()
     48 def addPlotData():
---> 49 plot.addData(i=O.iter,Sz=-triax.stress[triax.wall_right_id][0], Dz=-triax.strain[0])
     50
     51 plot.plots= {'Dz':('Sz')}

TypeError: 'instancemethod' object has no attribute '__getitem__'

regarding the other error: Lets solve this first ;)

Hello

"I am trying the same: Sz=-triax.stress(triax.wall_right_id)[0]"

It actually seems you're not....

See "triax.stress[triax.wall_right_id][0]" in your error message with [..] around "triax.wall_right_id" and not parenthesis (..) as kindly adviced several times by Jan...

Thanks Jerome for pointing that out!

I am a seriously bad programmer. should be more cautious with my typing of commands.

I ran the plot but it does not update.

Jan, Jerome.

I also ran the example script to check on the plot (https://github.com/yade/trunk/blob/master/examples/triax-tutorial/script-session1.py)

The strain there remains zero as the stress increases and later no change.
In my theoretical understanding, the stress strain curve should follow a proportional (linear) curve. Reach a maximum and then attain a strain hardening phase more or less. Do we have to make some suitable changes?

the script I am trying only gives me a strain hardening (constant zero stress - strain goes on increasing). what could be the error?

> I ran the plot but it does not update.

It updates, but the values of Dz ans Sz are 0 or NaN. You can check by
print plot.data

>I also ran the example script to check on the plot (https://github.com/yade/trunk/blob/master/examples/triax-tutorial/script-session1.py)
>The strain there remains zero as the stress increases and later no change.
In that script, internalCompaction=True is used meaning that strain is indeed not changing and stress is caused by particle growth.
You can also see the "movie" of the simulation.

cheers
Jan

Jan,
But do you not think that it is not modeling the behavior of the granular material appropriately in that case.

It should generally be like this : the stain increases with the stress proportionally, peaks, followed by a slight dip and then strain hardening occurs.

we are not able to achieve that here.

Hello,

From my point of view, the "we" is incorrect here since it suggest you're facing a YADE limitation, whereas there has been a great deal of published scientific papers showing appropriate stress-strain curves and done with YADE...

I thus would like to warmly advice you take the time to carefully go through all attributes of TriaxialStressController, and the example script at https://github.com/yade/trunk/blob/master/examples/triax-tutorial/script-session1.py.

You will see this internalCompaction attribute is intended for the first ("consolidation") step of triaxial tests imposed as an isotropic compression. During the deviatoric loading part where "the stain increases with the stress proportionally [...]", any YADE user should pay attention this attribute is switched/set to false, as in the above-mentioned script, at the 1st line of the "DEVIATORIC LOADING" section (commented in the state of the script on github, but this is actually no reason for not reading it) :

https://github.com/yade/trunk/blob/master/examples/triax-tutorial/script-session1.py#L169

Sorry, for using "we" before. I did not intend it to use it seriously though.
But your response was perfect Jerome. I included the deviatoric component of the script into the current code.

Earlier, I was trying it without the deviatoric part not noticing its significance to the initial set-up here (even more than that in a real lab experiment ;)). Now I can understand the functioning of the code better.

the full script is as follows. another question at this stage: What command to use in order to confine the spheres within the parallelopiped container. The spheres are moving out during the loading with my script. I am looking for a solution in the documentation while you answer ;)

from yade import pack,plot,qt
import math

num_spheres=10000 # number of spheres
young=3e10
compFricDegree = 30 # initial contact friction during the confining phase
finalFricDegree = 30 # contact friction during the deviatoric loading
mn,mx=Vector3(0,0,0),Vector3(0.1,0.1,0.1) # corners of the initial packing

O.materials.append(FrictMat(young=young,poisson=0.3,frictionAngle=radians(compFricDegree),density=2600,label='spheres'))
O.materials.append(FrictMat(young=young,poisson=0.5,frictionAngle=0,density=0,label='walls'))
walls=aabbWalls([mn,mx],thickness=0,material='walls')
wallIds=O.bodies.append(walls)

sp=pack.SpherePack()
sp.makeCloud(mn,mx,-1,0.5,num_spheres,False,0.40,seed=1) #"seed" make the "random" generation always the same
sp.toSimulation(material='spheres')

newton=NewtonIntegrator(damping=0.3)

O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb(),Bo1_Wall_Aabb()]),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
  [Ip2_FrictMat_FrictMat_FrictPhys()],
  [Law2_ScGeom_FrictPhys_CundallStrack()],label='iloop'
 ),
 TriaxialStressController(label='triax',
 maxMultiplier=1.+2e4/young, # spheres growing factor (fast growth)
 finalMaxMultiplier=1.+2e3/young, # spheres growing factor (slow growth)
 thickness = 0,
 stressMask =7,
 max_vel = 0.005,
 internalCompaction=True, # If true the confining pressure is generated by growing particles
),GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=50,timestepSafetyCoefficient=0.8),PyRunner(iterPeriod=100,command='addPlotData()'),
 newton
]

#O.dt=0.005*PWaveTimeStep()
#strain_rate=-0.05
triax.goal1=triax.goal2=triax.goal3=-5000

#Deviatoric loading
##We move to deviatoric loading, let us turn internal compaction off to keep particles sizes constant
triax.internalCompaction=False

## Change contact friction (remember that decreasing it would generate instantaneous instabilities)
setContactFriction(radians(finalFricDegree))

##set stress control on y and z, we will impose strain rate on x
triax.stressMask = 6
##now goal2 is the target strain rate
triax.goal2=-0.01
## we define the lateral stresses during the test, here the same 10kPa as for the initial confinement.
triax.goal1=-10000
triax.goal3=-10000
##Save temporary state in live memory. This state will be reloaded from the interface with the "reload" button.
O.saveTmp()

def addPlotData():
  plot.addData(i=O.iter,Sz=-triax.stress(triax.wall_right_id)[0], Dz=-triax.strain[0])

plot.plots= {'Dz':('Sz')}
plot.plot()

I tried changing the young's modulus to a higher value ( e.g., young=5e10) and to play with the time step ( O.dt=0.5*PWavetimeStep()). The particles still move out. What should I try to prevent this from happening? It affects the results.

> triax.stressMask = 6
> triax.goal2=-0.01
> triax.goal1=-10000
> triax.goal3=-10000

Are you sure these values fit together? I would **guess** for stressMask=6 (=0b110) it should be
triax.goal1=-0.01
triax.goal2=-10000
triax.goal3=-10000

Jan

Yes Jan,my bad.

However, it still does not give the correct (expected) results upon inducing the change required.
Am I doing something else wrong?

> it still does not give the correct (expected) results

what results it gives and what are correct (expected) results? :-)
You are using the "most simple" Law2_ScGeom_FrictPhys_CundallStrack, do you know about any "correct (expected)" results with this model?
Maybe with different material parameters it gives expected results?

cheers
Jan

Yes, I mean the strain just goes on increasing with a constant zero stress (please try running to cross-check). It should be increasing with the stress, peak and then go into a strain hardening mode.

I think the values I am trying here seem fine.

Is it possible to see the plot produced by the original script in [1] so I could know where I am going wrong?

[1] https://github.com/yade/trunk/blob/master/examples/triax-tutorial/script-session1.py.

> please try running to cross-check

to be sure I am trying the same code, please put here the current version

> Is it possible to see the plot produced by the original script in [1] so I could know where I am going wrong?

sure, just uncomment the section
#####################################################
### Example of how to record and plot data ###
#####################################################

cheers
Jan

Yes, the one I am trying is as follows: Perhaps I might be deviating but it looks syntactically alright :)
#triaxial test on dry sand

from yade import pack,plot,qt
import math

num_spheres=10000 # number of spheres
young=5e10
compFricDegree = 30 # initial contact friction during the confining phase
finalFricDegree = 30 # contact friction during the deviatoric loading
mn,mx=Vector3(0,0,0),Vector3(0.1,0.1,0.1) # corners of the initial packing

O.materials.append(FrictMat(young=young,poisson=0.3,frictionAngle=radians(compFricDegree),density=2600,label='spheres'))
O.materials.append(FrictMat(young=young,poisson=0.5,frictionAngle=0,density=0,label='walls'))
walls=aabbWalls([mn,mx],thickness=0,material='walls')
wallIds=O.bodies.append(walls)

sp=pack.SpherePack()
sp.makeCloud(mn,mx,-1,0.5,num_spheres,False,0.55,seed=1) #"seed" make the "random" generation always the same
sp.toSimulation(material='spheres')

newton=NewtonIntegrator(damping=0.3)
O.dt=0.1*PWaveTimeStep()

O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb(),Bo1_Wall_Aabb()]),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom(), Ig2_Wall_Sphere_ScGeom()],
  [Ip2_FrictMat_FrictMat_FrictPhys()],
  [Law2_ScGeom_FrictPhys_CundallStrack(),]
 ),
 TriaxialStressController(label='triax',
 maxMultiplier=1.+2e4/young, # spheres growing factor (fast growth)
 finalMaxMultiplier=1.+2e3/young, # spheres growing factor (slow growth)
 thickness = 0,
 stressMask =7,
 max_vel = 0.005,
 internalCompaction=True, # If true the confining pressure is generated by growing particles
),GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=50,timestepSafetyCoefficient=0.8),PyRunner(iterPeriod=100,command='addPlotData()'),
 newton
]

#O.dt=0.005*PWaveTimeStep()
#strain_rate=-0.05
triax.goal1=triax.goal2=triax.goal3=-5000

#Deviatoric loading
##We move to deviatoric loading, let us turn internal compaction off to keep particles sizes constant
triax.internalCompaction=False

## Change contact friction (remember that decreasing it would generate instantaneous instabilities)
setContactFriction(radians(finalFricDegree))

##set stress control on y and z, we will impose strain rate on x
triax.stressMask = 6
##now goal2 is the target strain rate
triax.goal1=-0.001
## we define the lateral stresses during the test, here the same 10kPa as for the initial confinement.
triax.goal2=-10000
triax.goal3=-10000
##Save temporary state in live memory. This state will be reloaded from the interface with the "reload" button.
O.saveTmp()

def addPlotData():
  plot.addData(i=O.iter,Sz=-triax.stress(triax.wall_right_id)[0], Dz=-triax.strain[0],Tx=math.sqrt(O.time))

plot.plots= {'Dz':('Sz'),'Tx':('Dz')}
plot.plot()

Yeah, I am also trying to run the example script but getting this error:

Traceback (most recent call last):
  File "/usr/bin/yade", line 182, in runScript
    execfile(script,globals())
  File "newtriaxtest.py", line 184
    O.engines=O.engines[0:5]+[PyRunner(iterPeriod=20,command='history()',label='recorder')]+O.engines[5:7]
                                                                                                         ^
IndentationError: unindent does not match any outer indentation level

The example script is running fine :)

I am not able to figure out the problem with my script ;)

There is zero stress because there are no interactions in your model. You can try
aabbExtrema()
and compare it to
O.bodies[triax.wall_right_id].state.pos
There is simply a gap between the wall and spheres.

Anyway, you created spheres by makeCloud, which makes a loose packing, so you would get very low stress anyway. Consider to use some "precompaction" either using randomDensePack or the TriaxialStressController itself.

Concerning the error, have you uncommented "if 1:" the line above?

> File "newtriaxtest.py", line 184

if you are referring to the original script, it is better to use that one also for testing, e.g. here the error message, line 184 point to nowhere in the original script..

cheers
Jan

thanks Jan.
1) I am trying with the randomDense pack option. could you please check this script.( It takes a lot of time for the periodic compression though ;))

#triaxial test on dry sand

from yade import pack,plot,qt
import math

num_spheres=5000 # number of spheres
young=5e10
compFricDegree = 30 # initial contact friction during the confining phase
finalFricDegree = 30 # contact friction during the deviatoric loading
mn,mx=Vector3(0,0,0),Vector3(0.1,0.1,0.1) # corners of the initial packing

O.materials.append(FrictMat(young=young,poisson=0.3,frictionAngle=radians(compFricDegree),density=2600,label='spheres'))
O.materials.append(FrictMat(young=young,poisson=0.5,frictionAngle=0,density=0,label='walls'))
walls=aabbWalls([mn,mx],thickness=0,material='walls')
wallIds=O.bodies.append(walls)

pred=pack.inAlignedBox(mn,mx)
sp=pack.randomDensePack(pred,radius=0.5,spheresInCell=num_spheres)
O.bodies.append(sp)

newton=NewtonIntegrator(damping=0.1)
O.dt=0.1*PWaveTimeStep()

O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb(),Bo1_Wall_Aabb()]),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom(), Ig2_Wall_Sphere_ScGeom()],
  [Ip2_FrictMat_FrictMat_FrictPhys()],
  [Law2_ScGeom_FrictPhys_CundallStrack(),]
 ),
 TriaxialStressController(label='triax',
 maxMultiplier=1.+2e4/young, # spheres growing factor (fast growth)
 finalMaxMultiplier=1.+2e3/young, # spheres growing factor (slow growth)
 thickness = 0,
 stressMask =7,
 max_vel = 0.005,
 internalCompaction=True, # If true the confining pressure is generated by growing particles
),GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=50,timestepSafetyCoefficient=0.8),PyRunner(iterPeriod=100,command='addPlotData()'),
 newton
]

#O.dt=0.005*PWaveTimeStep()
#strain_rate=-0.05
triax.goal1=triax.goal2=triax.goal3=-5000

#Deviatoric loading
##We move to deviatoric loading, let us turn internal compaction off to keep particles sizes constant
triax.internalCompaction=False

## Change contact friction (remember that decreasing it would generate instantaneous instabilities)
setContactFriction(radians(finalFricDegree))

##set stress control on y and z, we will impose strain rate on x
triax.stressMask = 6
##now goal1 is the target strain rate
triax.goal1=-0.01
## we define the lateral stresses during the test, here the same 10kPa as for the initial confinement.
triax.goal2=-10000
triax.goal3=-10000
##Save temporary state in live memory. This state will be reloaded from the interface with the "reload" button.
O.saveTmp()

def addPlotData():
  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],
      i=O.iter)

plot.plots={'i':('e11','e22','e33',None,'s11','s22','s33')}
plot.plots={'e11':('s11'),'e22':('s22'),'e33':('s33')} #,'s22','s33')}
plot.plot()

I stopped the compilation in between because I realized I have not specified the material properties for the spheres added to the box through the following segment:

pred=pack.inAlignedBox(mn,mx)
sp=pack.randomDensePack(pred,radius=0.5,spheresInCell=num_spheres)
O.bodies.append(sp)

How do I assign the label 'spheres' (which contains the material properties) to the above spherical particles?

> could you please check this script

> mn,mx=Vector3(0,0,0),Vector3(0.1,0.1,0.1)
> pred=pack.inAlignedBox(mn,mx)
> sp=pack.randomDensePack(pred,radius=0.5,...)

using spheres (radius=0.5) larger than predicate (0.1) results in zero spheres..

> It takes a lot of time for the periodic compression though

use lower value of spheresInCell, e.g. 300 should be ok

> How do I assign the label 'spheres' (which contains the material properties) to the above spherical particles?

e.g.
sp=pack.randomDensePack(...,material='spheres')

Jan

Yes, Jan. Sorry another typo from me (as usual :P ). It should have been 0.05.

But I have been trying the script with the corrections and it works fine now. My only trouble is with the stress-strain plot. It has a lot of fluctuations. Perhaps I need to find better set of data for my case. Any suggestions with that ?

Also, is there a possibility to insert a certain porosity with the RandomDensePack ?

Another issue is that when I increase the number of sphere to say 5000., increase box dimensions to (1,1,1) and particle size to 0.005, I get the following error:

Running script drysandoedo.py
No suitable packing in database found, running PERIODIC compression
/usr/lib/x86_64-linux-gnu/yade/py/yade/pack.py:296: FutureWarning: The default behavior will change; specify returnSpherePack=True for the new behavior, and False to get rid of this warning (your code will break in the future, however). The returned SpherePack object can be added to the simulation using SpherePack.toSimulation()
  warnings.warn('The default behavior will change; specify returnSpherePack=True for the new behavior, and False to get rid of this warning (your code will break in the future, however). The returned SpherePack object can be added to the simulation using SpherePack.toSimulation()',category=FutureWarning)
[[ ^L clears screen, ^U kills line. F12 controller, F11 3d view (use h-key for showing help), F10 both, F9 generator, F8 plot. ]]
[0;34mYade [[1;34m1[0;34m]: [0mICE default IO error handler doing an exit(), pid = 6298, errno = 32

> My only trouble is with the stress-strain plot. It has a lot of fluctuations. Perhaps I need to find better set of data for my case. Any suggestions with that ?

it depends on "a lot", certain fluctuations are expected since the simulation is dynamic. I would say strain rate would influence this significantly (at least for higher values).
Also compare the fluctuations to the maximum stress reached, maybe in the beginning of simulation there are seemingly a lot of fluctuations, but they are negligible with respect to desired stress.

> Also, is there a possibility to insert a certain porosity with the RandomDensePack ?

no, randomDensePack is predefined function, but you can quite easily create your own function based on the code of randomDensePack. e.g. playinh with material friction angle (which is hardcoded in randomDensePack), you can modify resulting porosity

> I get the following error:

it is no error, just info / warning

> No suitable packing in database found, running PERIODIC compression

info what is happening

> utureWarning: The default behavior will change; specify returnSpherePack=True ...

you can change your code to (be careful with the material):

sp=pack.randomDensePack(..., returnSpherePack=True) # material has no effect in this case
sp.toSimulation(material='cpm') # put material specification here

> mICE default IO error handler doing an exit()

this is after you exit Yade? If so, I would not bother with it..

Jan

Will try making the changes :)

Regarding the last part Jan.

It exits the simulation. I cannot see the set-up anymore. I get a new line on the emulator. E.g. in my case - zg@ZG-Yade ~/scriptf $

I have no idea, try [1,2], googling or ask a separate question if nothing helps..
Jan

[1] https://stackoverflow.com/questions/20621070/ice-default-io-error-handler-doing-an-exit-pid-11281-errno-4
[2] https://stackoverflow.com/questions/39270988/ice-default-io-error-handler-doing-an-exit-pid-errno-32-when-running

I will check that Jan but after 50 days. Going on vacation now after 2 years.

Please do not close the question ;)