Some questions related to the Triaxial code by Bruno Chareyre

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.

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

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 ============')
nRead=readParamsFromTable(
 num_spheres=20000,
 compFricDegree = 30,
 key='_triax_base_',
 unknownOk=True
)

from yade.params import table

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=radians(compFricDegree),density=2600,label='spheres'))
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])
    mean_rad = pow(0.09*volume/num_spheres,0.3333)
    c1=pack.SpherePack([((-0.2*mean_rad,0,0),0.5*mean_rad),((0.2*mean_rad,0,0),0.5*mean_rad)])
    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)
    O.bodies.append([sphere(center,rad,material='spheres') for center,rad in sp])

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

print ('============ DEFINING TRIAXIAL TEST ============')
triax=TriaxialStressController(
 maxMultiplier=1.+2e4/young,
 finalMaxMultiplier=1.+2e3/young,
 thickness = 0,
 stressMask = 7,
 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
if nRead==0: yade.qt.Controller(), yade.qt.View()
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.stressmask=7
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]

######################################################
######### DEVIATORIC LOADING #########

print ('============ APPLYING DEVIATORIC LOADING ============')
triax.internalCompaction=False
setContactFriction(radians(finalFricDegree))
triax.stressMask = 5
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():
 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)

########################################################
######### 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.