# yade-users team mailing list archive

## [Question #684319]: How to reach target confining stress in triaxial test

New question #684319 on Yade:

Hi,
I am using [1] for simulating a triaxial test. The following MWE is from [1] and only some parameters are modified (damp, targetPorosity, stabilityThreshold).
My question is, after the confining process is completed, the confining pressure doesn't reach our desired pressure which is 10kPa, the mean stress and  other 3 stresses are:
meanStress:  9997.34994518; top: 9649.20616727 ;right: 9842.19217112 ;front :10510.4130631.
It is very nice for meanStress, but the value of other 3 directions are far from the target goal, and this results in the stress-strain curve doesn't start from origin point.
Here is the MWE, it takes about 30 seconds to run it.
#####

num_spheres=2000,
compFricDegree = 30,
key='_triax_base_',
unknownOk=True
)

num_spheres=table.num_spheres
key=table.key
targetPorosity = 0.7
compFricDegree = table.compFricDegree
finalFricDegree = 30
rate=-0.01
damp=0.4
stabilityThreshold=0.001
young=5e7

confinement=10e3

mn,mx=Vector3(0,0,0),Vector3(0.7,1.4,0.7)

O.materials.append(FrictMat(young=young,poisson=0.5,frictionAngle=0,density=0,label='walls'))

## create walls around the packing
walls=aabbWalls([mn,mx],thickness=0,material='walls')
wallIds=O.bodies.append(walls)

sp=pack.SpherePack()
sp.makeCloud(mn,mx,-1,0.3333,num_spheres,False, 0.95,seed=1)

Gl1_Sphere.quality=3

triax=TriaxialStressController(

maxMultiplier=1.+4e-3,
finalMaxMultiplier=1.+4e-4,
thickness = 0,

internalCompaction=True,
)
newton=NewtonIntegrator(damping=damp)

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()]
),
## We will use the global stiffness of each body to determine an optimal timestep (see https://yade-dem.org/w/images/1/1b/Chareyre&Villard2005_licensed.pdf)
GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
triax,
TriaxialStateRecorder(iterPeriod=100,file='WallStresses'+table.key),
newton
]

Gl1_Sphere.stripes=0

triax.goal1=triax.goal2=triax.goal3=-confinement

while 1:
O.run(1000, True)
unb=unbalancedForce()
print 'unbF:',unb,' meanStress: ',-triax.meanStress,'top:',-triax.stress(triax.wall_top_id)[1]
if unb<stabilityThreshold and abs(-confinement-triax.meanStress)/confinement<0.001:
break

print "###      Isotropic state saved      ###"
##### below is what I tried, but it seems doesn't work
# while 1:
#   O.run(100, True)
#   print ' meanStress: ',-triax.meanStress,'top:',-triax.stress(triax.wall_top_id)[1],'front:',-triax.stress(triax.wall_front_id)[2],\
# 	  'right:',-triax.stress(triax.wall_right_id)[0]
#   if abs(-confinement-triax.meanStress)/confinement<0.001 and \
# 		  abs(-confinement-triax.stress(triax.wall_right_id)[0])/confinement<0.01 and \
# 		  abs(-confinement-triax.stress(triax.wall_top_id)[1])/confinement<0.01 and \
# 		  abs(-confinement-triax.stress(triax.wall_front_id)[2])/confinement<0.01:
# 	break

import sys
while triax.porosity>targetPorosity:
# we decrease friction value and apply it to all the bodies and contacts
compFricDegree = 0.95*compFricDegree
print "\r Friction: ",compFricDegree," porosity:",triax.porosity,
sys.stdout.flush()
# while we run steps, triax will tend to grow particles as the packing
# keeps shrinking as a consequence of decreasing friction. Consequently
# porosity will decrease
O.run(500,1)

print "###    Compacted state saved      ###"

print 'meanStress: ',-triax.meanStress, 'top', -triax.stress(triax.wall_top_id)[1], 'right',-triax.stress(triax.wall_right_id)[0],'front',-triax.stress(triax.wall_front_id)[2]

triax.internalCompaction=False
triax.goal2=rate
triax.goal1=-confinement
triax.goal3=-confinement

newton.damping=0.4
O.saveTmp()

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],
dev=-triax.stress(triax.wall_top_id)[1]-10000,
i=O.iter)

def stopIfDamaged():
if -triax.strain[1]>0.1:
O.pause()
plot.saveDataTxt('Frict_10kpa_num{}_CaAtGap'.format(num_spheres))

if 1:
O.engines=O.engines[0:5]+[PyRunner(iterPeriod=20,command='history()',label='recorder')]+[PyRunner(iterPeriod=20,command='stopIfDamaged()',label='checkDamage')]+O.engines[5:7]
plot.plots={'e22':('s11','s22','s33'),'e22 ':('dev')}
plot.labels={'e22':'$\epsilon_{a}$','s11':'$\sigma_{11}$','s22':'$\sigma_{22}$','s33':'$\sigma_{33}$'}
plot.plot()
###########
The commented part of code is what I tried to let the stress reach our target value. but it seems doesn't work since the top, front and right stress doesn't change even after run it for several hours.

Could you please give some instruction about how to reach target confining stress?
Thanks
Leonard