# Not rational result in the Triaxial code

Asked by ehsan benabbas on 2020-01-21

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 run triaxial simulation on my specimen.

Unfortunately, the result is not rational and thus my questions are:

1- Why the "q" (deviatoric stress) is positive? (at the end of the deviatoric part)
2- Why axial deformation (axial strain, volumetric strain) is zero? (at the end of both the confining and deviatoric parts)
3- Why sigma2 is less than SigmaIso in the code? (at the end of the deviatoric part)
4- How can I save the packing (particles' coordinates as X,Y, Z and their radius as R) in a text file after make the packing?

%%%%%%%%%%%%%%%%%%%%%%%%

I only made the following changes in [1]:

num_spheres=20000,
targetPorosity = 0.4
young=200000000
poisson=1
mn,mx=Vector3(0,0,0),Vector3(0.02,0.02,0.02)
sigmaIso=-500000
particleDensity=2000
sp=pack.SpherePack()
clumps=False
volume = (mx[0]-mn[0])*(mx[1]-mn[1])*(mx[2]-mn[2])
sp.makeCloud(mn,mx,0.0002,0.5,num_spheres,False, 0.95,seed=1)
triax.goal1 = sigmaIso
triax.goal2 = sigmaIso
triax.goal3 = sigmaIso
f unb<stabilityThreshold and abs(sigmaIso-triax.meanStress)/sigmaIso<0.001:
triax.goal1 = sigmaIso
triax.goal2 = rate
triax.goal3 = sigmaIso

%%%%%%%%%%%%%%%%%%%%%%%%

And the output is:

ehsan@ehsan:~/Desktop$/home/ehsan/yade/install/bin/yade-2019-08-08.git-775ae74 Q.py Welcome to Yade 2019-08-08.git-775ae74 Using python version: 3.6.9 (default, Nov 7 2019, 10:44:02) [GCC 8.3.0] TCP python prompt on localhost:9000, auth cookie `cyueka' XMLRPC info provider on http://localhost:21000 Running script Q.py ************** START ************** ============ DEFINING VARIABLES ============ ============ DEFINING FUNCTIONS ============ ============ DEFINING MATERIALS ============ ============ DEFINING PACKING ============ ============ DEFINING TRIAXIAL TEST ============ ============ DEFINING ENGINES ============ The constructor with a shareWidget is deprecated, use the regular contructor instead. Number of elements: 20006 Box Volume: 3.714574414115203e-81 Box Volume calculated: 8.000000000000001e-06 ============ APPLYING CONFINING PRESSURE ============ ~~~~~~~~~~~~~ Phase_01: Converging to Isotropic Compression, 50kPa ~~~~~~~~~~~~~ unbalanced force: nan mean stress engine: 0.0 mean stress (Calculated): -0.0 porosity= 0.8920863623534938 void ratio= 8.266669364586804 . . . ~~~~~~~~~~~~~ Phase_01: Converging to Isotropic Compression, 50kPa ~~~~~~~~~~~~~ unbalanced force: 5e-05 mean stress engine: 0.0 mean stress (Calculated): -0.0 porosity= 0.874623180910073 void ratio= 6.975956059969474 Axial Strain 0.0 ################## Isotropic phase is finished and saved successfully ################## ============ REACHING TARGET POROSITY ============ Friction: 28.5 porosity: 0.874623180910073 Friction: 27.075 porosity: 0.8708050821138507 . . . Friction: 2.0832852056170093 porosity: 0.42156339492430284 Friction: 1.9791209453361587 porosity: 0.40486546436292276 Final Number of elements: 20006 Final Box Volume: 8.000000000000001e-06 Final Box Volume calculated: 8.000000000000001e-06 Final porosity= 0.38946593512762423 Final void ratio= 0.6379102453669592 ################## Target porosity is reached and compacted state saved successfully ################## ============ APPLYING DEVIATORIC LOADING ============ step= 59000 unbalanced force: 0.22470167357987816 sigma2: -121121.56365292308 q= 378878.4363470769 Axial Deformation (%) 0.0 ~~~~~~~~~~~~~ Phase_02: Converging to Deviatoric Compression, Strain Rate ~~~~~~~~~~~~~ ################## Deviatoric phase is finished and saved successfully ################## ============ RECORD AND PLOT DATA ============ /home/ehsan/yade/install/lib/x86_64-linux-gnu/yade-2019-08-08.git-775ae74/py/yade/plot.py:444: MatplotlibDeprecationWarning: The 'verts' kwarg was deprecated in Matplotlib 3.0 and will be removed in 3.2. Use 'marker' instead. scatter=pylab.scatter(scatterPt[0] if not math.isnan(scatterPt[0]) else 0,scatterPt[1] if not math.isnan(scatterPt[1]) else 0,s=scatterSize,color=line.get_color(),**scatterMarkerKw) ************** END ************** Analysis has been taken for 661.858 seconds or 11.030966666666666 minutes [[ ^L clears screen, ^U kills line. F12 controller, F11 3D view (press "h" in 3D view for help), F10 both, F9 generator, F8 plot. ]] In [1]: %%%%%%%%%%%%%%%%%%%%%%%% I will copy the whole code in the first comment. Thank you for your help. Ehsan ## Question information Language: English Edit question Status: Solved For: Yade Edit question Assignee: No assignee Edit question Solved by: Jan Stránský Solved: 2020-01-23 Last query: 2020-01-23 Last reply: 2020-01-22  ehsan benabbas (ehsanben) said on 2020-01-21: #1 The code is: ############################################################################################################################ ######### TRIAXIAL PROBLEM, Y IS THE VERTICAL AXIS, X IS THE RIGHT AXIS, Z IS THE FRONT AXIS ######### ############################################################################################################################ print ('************** START **************') import numpy as np import datetime import time import math from yade import 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.4 compFricDegree = table.compFricDegree finalFricDegree = 30 rate=-0.02 damp=0.2 thick=0 stabilityThreshold=0.01 young=200000000 poisson=1 mn,mx=Vector3(0,0,0),Vector3(0.02,0.02,0.02) sigmaIso=-500000 particleDensity=2000 ###################################################### ################# DEFINING FUNCTIONS ################# print ('============ DEFINING FUNCTIONS ============') from yade import plot 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, triax.porosity) ###################################################### ################# DEFINING MATERIALS ################# print ('============ DEFINING MATERIALS ============') O.materials.append(FrictMat(young=young,poisson=poisson,frictionAngle=radians(compFricDegree),density=particleDensity,label='spheres')) O.materials.append(FrictMat(young=young,poisson=poisson,frictionAngle=0,density=0,label='walls')) #################################################### ################# DEFINING PACKING ################# print ('============ DEFINING PACKING ============') walls=aabbWalls([mn,mx],thickness=thick,material='walls') wallIds=O.bodies.append(walls) from yade import pack sp=pack.SpherePack() clumps=False volume = (mx[0]-mn[0])*(mx[1]-mn[1])*(mx[2]-mn[2]) sp.makeCloud(mn,mx,0.0002,0.5,num_spheres,False, 0.95,seed=1) O.bodies.append([sphere(center,rad,material='spheres') for center,rad in sp]) O.save('packing'+key+'.xml') ########################################################## ################# DEFINING TRIAXIAL TEST ################# print ('============ DEFINING TRIAXIAL TEST ============') triax=TriaxialStressController( maxMultiplier=1.+2e4/young, finalMaxMultiplier=1.+2e3/young, thickness = thick, stressMask = 7, internalCompaction=True, ) #################################################### ################# 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, PyRunner(iterPeriod=20,command='history()',label='recorder'), TriaxialStateRecorder(iterPeriod=100,file='WallStresses'+table.key), NewtonIntegrator(damping=damp,label="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 ('~~~~~~~~~~~~~ Phase_01: Converging to Isotropic Compression, 50kPa ~~~~~~~~~~~~~') 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)) if unb<stabilityThreshold and abs(sigmaIso-triax.meanStress)/sigmaIso<0.001: break O.save('confinedPhase'+key+'.xml') e22Check=triax.strain[1] print ('Axial Strain',e22Check) print ('################## Isotropic phase is finished and saved successfully ##################') ############################################################# ################# REACHING TARGET POROSITY ################## print ('============ REACHING TARGET POROSITY ============') import sys #this is only for the flush() below while triax.porosity>targetPorosity: compFricDegree = 0.95*compFricDegree setContactFriction(radians(compFricDegree)) print ("\r Friction: ",compFricDegree," porosity:",triax.porosity) sys.stdout.flush() O.run(1000,1) print ('Final Number of elements: ', len(O.bodies)) print ('Final Box Volume: ', triax.boxVolume) print ('Final Box Volume calculated: ', volume) print ('Final porosity=',triax.porosity) print ('Final void ratio=',triax.porosity/(1-triax.porosity)) O.save('compactedState'+key+'.yade.gz') print ('################## Target porosity is reached and compacted state saved successfully ##################') ###################################################### ################# 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 ##################') ######################################################## ################# RECORD AND PLOT DATA ################# print ('============ RECORD AND PLOT DATA ============') if 1: O.engines=O.engines[0:5]+[PyRunner(iterPeriod=1,command='history()',label='recorder')]+O.engines[5:7] else: O.engines[4]=PyRunner(iterPeriod=1,command='history()',label='recorder') O.run(100,True) plot.plots={'e22':('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)

#####################################################
################# RECORD Micro DATA #################

data = []
for i in O.interactions:
fn = i.phys.normalForce
fs = i.phys.shearForce
cp = i.geom.contactPoint
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
cp,normal,branch,fn,fs = [tuple(v) for v in (cp,normal,branch,fn,fs)] # Vector3 -> tuple
d = dict(cp=cp,normal=normal,branch=branch,fn=fn,fs=fs)
data.append(d)
# new data contains the information, you can save it e.g. as JSON
import json
with open("interactions.json","w") as f:
json.dump(data,f)

#######################################
################# END #################

print ('************** END **************')
O.realtime
print ('Analysis has been taken for',O.realtime, 'seconds or', O.realtime/60, 'minutes')

 Jan Stránský (honzik) said on 2020-01-21: #2

Hello,

please provide a working code. I have tried it, resulting in
File "e.py", line 55
triax.porosity)
^
SyntaxError: positional argument follows keyword argument

> the result is not rational

the results of the original script without changes are rational or not?

> 4- How can I save the packing (particles' coordinates as X,Y, Z and their radius as R) in a text file after make the packing?
###
export.text(fileName) # [1]
###

cheers
Jan

 ehsan benabbas (ehsanben) said on 2020-01-21: #3

Hello,

>> File "e.py", line 55
triax.porosity)

Yes, you are right. Please change it to this:
n = triax.porosity)

>> the results of the original script without changes are rational or not?
Actually the part that gets the axial deformation and deviatoric stress (q) is not from the code, I added those lines. But yes, even on the main code the axial deformation is still 0 (on the line I added to be printed) and the "q" is positive and sigma2 is less than confining pressure. Those lines can be found in the "DEVIATORIC LOADING" section of the code. However, these parameters are shown before pressing the "Play" button.

export.text(fileName) # [1]

That work's perfectly. Thank you so much. This is exactly like what I asked you about micro variables saving in a text file. I am looking for something like this, clear and nice.

Regards,
Ehsan

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

I still need an answer for these 3 questions:

1- Why the "q" (deviatoric stress) is positive? (at the end of the deviatoric part)
2- Why axial deformation (axial strain, volumetric strain) is zero? (at the end of both the confining and deviatoric parts)
3- Why sigma2 is less than SigmaIso in the code? (at the end of the deviatoric part)

 Jan Stránský (honzik) said on 2020-01-22: #5

0)
you have no running in the deviatoric part, so "at the end of deviatoric part" = "at the end of reaching target porosity"

1+3) (same problem)
sigma2 = axialS, q = axialS-sigmaIso

because sigma2=axialS has no relation to sigmaIso (see below)

> if unb<stabilityThreshold and abs(sigmaIso-triax.meanStress)/sigmaIso<0.001:

since sigmaIso is negative, the fraction in the second condition is always negative and therefore the second condition is always met.
So the loop is stopped based only on unbalanced force= with arbitrary mean stress.
In your case triax.meanStress = 0 (according to the output)

2)
> Why axial deformation ... is zero?

confining - because internalCompaction=True (default value)
deviatoric - you have no running, the strain remains 0

cheers
Jan

 ehsan benabbas (ehsanben) said on 2020-01-23: #6

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