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.

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

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
O.materials.append(FrictMat(young=young,poisson=poisson,frictionAngle=radians(compFricDegree),density=particleDensity,label='spheres'))
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])
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 : #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 : #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?
###
from yade import export
export.text(fileName) # [1]
###

cheers
Jan

[1] https://yade-dem.org/doc/yade.export.html#yade.export.text

ehsan benabbas (ehsanben) said : #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.

>> from yade import export
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 : #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)

any answer please?

Best Jan Stránský (honzik) said : #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 : #6

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