# Some questions related to the Triaxial code by Bruno Chareyre

Asked by ehsan benabbas on 2020-01-07

Hi everyone,

I am using Ubuntu 18.04, and Yade 2019-08-08.git-775ae74

I use the Triaxial code by Bruno Chareyre [1] to understand how it works.

My goal is to code a Triaxial test with 20000 particles in different sizes (sand) under different strain and stress paths and get all inter-particle forces, branch vectors, contact normals, and fabric tensor. (with the properties defined in the code)

This is the code I am running. My questions are at the end of this message.

###################################################
############################################################################################################################
######### TRIAXIAL PROBLEM, Y IS THE VERTICAL AXIS, X IS THE RIGHT AXIS, Z IS THE FRONT AXIS #########
############################################################################################################################

import datetime
aa = datetime.datetime.now()
print ('************** START **************')
import numpy as np
import math
from yade import pack, plot, qt, export, utils
from datetime import datetime

######################################################
######### DEFINING VARIABLES #########

print ('============ DEFINING VARIABLES ============')
num_spheres=20000,
compFricDegree = 30,
key='_triax_base_',
unknownOk=True
)

num_spheres=table.num_spheres
key=table.key
targetPorosity = 0.35
compFricDegree = table.compFricDegree
finalFricDegree = 30
rate=-0.02
damp=0.2
stabilityThreshold=0.01
young=5e6
poisson=0.3
mn,mx=Vector3(0,0,0),Vector3(20,20,20)
sigmaIso=-50e3

######################################################
######### DEFINING MATERIALS #########

print ('============ DEFINING MATERIALS ============')
O.materials.append(FrictMat(young=young,poisson=poisson,frictionAngle=0,density=0,label='frictionlesswalls'))

####################################################
######### DEFINING PACKING #########

print ('============ DEFINING PACKING ============')
frictionlesswalls=aabbWalls([mn,mx],thickness=0,material='frictionlesswalls')
wallIds=O.bodies.append(frictionlesswalls)
sp=pack.SpherePack()
clumps=False
if clumps:
volume = (mx[0]-mn[0])*(mx[1]-mn[1])*(mx[2]-mn[2])
sp.makeClumpCloud(mn,mx,[c1],periodic=False)
sp.toSimulation(material='spheres')
O.bodies.updateClumpProperties()
else:
volume = (mx[0]-mn[0])*(mx[1]-mn[1])*(mx[2]-mn[2])
sp.makeCloud(mn,mx,-1,0.3333,num_spheres,False, 0.95,seed=1)

##########################################################
######### DEFINING TRIAXIAL TEST #########

print ('============ DEFINING TRIAXIAL TEST ============')
triax=TriaxialStressController(
maxMultiplier=1.+2e4/young,
finalMaxMultiplier=1.+2e3/young,
thickness = 0,
internalCompaction=True,
)
newton=NewtonIntegrator(damping=damp)

####################################################
######### DEFINING ENGINES #########

print ('============ DEFINING ENGINES ============')
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()]
),
GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
triax,
TriaxialStateRecorder(iterPeriod=100,file='WallStresses'+table.key),
newton
]
Gl1_Sphere.stripes=True
print ('Number of elements: ', len(O.bodies))
print ('Box Volume: ', triax.boxVolume)
print ('Box Volume calculated: ', volume)

###############################################################
######### APPLYING CONFINING PRESSURE #########

print ('============ APPLYING CONFINING PRESSURE ============')
triax.goal1=sigmaIso
triax.goal2=sigmaIso
triax.goal3=sigmaIso
while 1:
O.run(1000, True)
unb=unbalancedForce()
meanS=(triax.stress(triax.wall_right_id)[0]+triax.stress(triax.wall_top_id)[1]+triax.stress(triax.wall_front_id)[2])/3
print ('unbalanced force:',unb,' mean stress engine: ',triax.meanStress,' mean stress (Calculated): ',meanS)
print ('porosity=',triax.porosity)
print ('void ratio=',triax.porosity/(1-triax.porosity))
print ('~~~~~~~~~~~~~ Phase_01: Converging to Isotropic Compression, 50kPa ~~~~~~~~~~~~~')
if unb<stabilityThreshold and abs(sigmaIso-triax.meanStress)/sigmaIso<0.001:
break
O.save('confinedPhase'+key+'.xml')
print ('################## Isotropic phase is finished and saved successfully ##################')
e22Check=triax.strain[1]

######################################################

triax.internalCompaction=False
triax.goal1=sigmaIso
triax.goal2=rate
triax.goal3=sigmaIso
newton.damping=0.1
unb=unbalancedForce()
axialS=triax.stress(triax.wall_top_id)[1]
print ('step=', O.iter, 'unbalanced force:',unb,' sigma2: ',axialS, 'q=', axialS-sigmaIso)
print ('axial deformation (%)', (triax.strain[1]-e22Check)*100)
print ('~~~~~~~~~~~~~ Phase_02: Converging to Deviatoric Compression, Strain Rate ~~~~~~~~~~~~~')
O.save('final.xml')
print ('################## Deviatoric phase is finished and saved successfully ##################')

######################################################
######### DEFINING FUNCTIONS #########

print ('============ DEFINING FUNCTIONS ============')
def history():
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)

########################################################
######### RECORD AND PLOT DATA #########

print ('============ RECORD AND PLOT DATA ============')
O.run(5000,True)
plot.plots={'e22':('s11','s22','s33',None,'ev')}
plot.labels={'s11':'$\sigma_{11}$' , 's22':'$\sigma_{22}$' , 's33':'$\sigma_{33}$' , 'e22':'$\epsilon_{22}$' , 'ev':'$\epsilon_{V}$'}
plot.plot()
plot.saveDataTxt('results'+key)
plot.saveGnuplot('plotScript'+key)
print ('************** END **************')
bb = datetime.datetime.now()
cc = bb - aa
print( int(cc.total_seconds()))
print(elapsed_time)
###################################################

Question:
1- When the simulation is finished? When I run the code, It takes some minutes to show the ************** END ************** message on the screen, while the Play button has not been chosed in the Yade window yet. If I press the Play button, it takes some hours and then I get run errors.
2- In the code we have defined the strain rate for the deviatoric part, what's the criterion to finish the test? does it work with damping? I have read about it but it's not clear for me yet.
3- I want to get the time of running but I get the error for it. I have used datetime function for it as it is in the code.
4- I don't get the conventional graph of triaxial test. the graph window shows up but it is just a point.
5- If I get the simple triaxial test correctly, how can I get the fabric tensor, inter-particle forces, branch vectors, and contact normals as the raw data?

Thank you.

## Question information

Language:
English Edit question
Status:
Solved
For:
Assignee:
No assignee Edit question
Solved by:
Jan Stránský
Solved:
2020-01-09
Last query:
2020-01-09
2020-01-08
 Jan Stránský (honzik) said on 2020-01-07: #1

Hello,

> then I get run errors.
> I get the error

please be more specific (ideally copy-paste the error here). Especially if it "it takes some hours", which nobody will like to test on his computer.

if the error is something like "elapsed_time is not defined", it is just because it is not defined anywhere in the code.

1)
> When the simulation is finished?

I would say when "************** END **************" is printed..

> > When I run the code, It takes some minutes to show the ************** END ************** message on the screen, while the Play button has not been chosed
>
> while 1:
> O.run(1000, True)

because you have O.run called in the code

2)
> what's the criterion to finish the test?

O.run(5000,True)

> does it work with damping?

newton=NewtonIntegrator(damping=damp)

> I have read about it but it's not clear for me yet.

please be more specific what is not clear, otherwise it is very difficult to give some reasonable suggestion

3)
> I want to get the time of running

you can use O.realtime, which is the time from the script start to the time of calling it

5)
> the fabric tensor

utils.fabricTensor()

>inter-particle forces

for i in O.interactions:
fn = i.phys.normalForce
fs = i.phys.shearForce

> branch vectors, and contact normals

for i in O.interactions:
n = i.geom.normal

does "branch vector" and "contact normal" differ?

> as the raw data

please be more specific on what "raw data" is

cheers
Jan

 ehsan benabbas (ehsanben) said on 2020-01-08: #2

Dear Jan, Thanks for your very much response. Just for clarification

1) I would say when "************** END **************" is printed..
>> So the output is ready and I don't need to choose the Play button, right? it does not make any difference in result. is that correct?

2) newton=NewtonIntegrator(damping=damp)
>> I don't know what this damping does exactly. It damps the rate of strain? I think this one is the criterion of finishing the test.

3) you can use O.realtime, which is the time from the script start to the time of calling it
>> So I just need to put O.realtime at the end of my code. right?

5) does "branch vector" and "contact normal" differ?
>> Yes, branch vectors connect the centroids of 2 particles while contact normals are the exterior unit normal vector at contact. The fabric tensor can be defined based on both with different formulations. Can I get the fabric tensor directy by yade?

Raw data >>> I mean the aforementioned micromechanical parameters in a text file which I can use for post processing.

Is there anyone who can help me to find the answer of question 4?

Thank you very much.

 Jan Stránský (honzik) said on 2020-01-08: #3

> So the output is ready and I don't need to choose the Play button, right? it does not make any difference in result. is that correct?

there is no code to add further results afterwards.
Of course the state of simulation will probably change

> I don't know what this damping does exactly. It damps the rate of strain?

> I think this one is the criterion of finishing the test.

there are a few stop conditions in the script.

For isotropic confinement part, it is based on unbalanced force:
unb=unbalancedForce()
if unb<stabilityThreshold and abs(sigmaIso-triax.meanStress)/sigmaIso<0.001:
break

The deviatoric part is just number of time steps:
O.run(5000,True)

> So I just need to put O.realtime at the end of my code. right?

it is one option, yes. You can of course use time or datetime module or anything else, depending on your needs

> branch vectors connect the centroids of 2 particles while contact normals are the exterior unit normal vector at contact. The fabric tensor can be defined based on both with different formulations.

right, then
for i in O.interactions:
normal = i.geom.normal
b1,b2 = [O.bodies[id] for id in (i.id1,i.id2)]
p1,p2 = [b.state.pos for b in (b1,b2)]
branch = p2 - p1 # maybe multiplied by -1, depending on what you do with it

> Can I get the fabric tensor directy by yade?

or you can compute it "manually"

> a text file which I can use for post processing.

it is still too ambiguous.
In python, you have access to the data and in python it is easy to save data to files. Yade provides in export module utilities to save to various formats, or you can just save it "manually" to the format you need.

cheers
Jan

 ehsan benabbas (ehsanben) said on 2020-01-09: #4

Thanks Jan Stránský, that solved my question.