# How to reach target confining stress in triaxial test

Hi,
I am using  for simulating a triaxial test. The following MWE is from  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)
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),'front:',-triax.stress(triax.wall_front_id),\
# 'right:',-triax.stress(triax.wall_right_id)
# if abs(-confinement-triax.meanStress)/confinement<0.001 and \
# abs(-confinement-triax.stress(triax.wall_right_id))/confinement<0.01 and \
# abs(-confinement-triax.stress(triax.wall_top_id))/confinement<0.01 and \
# abs(-confinement-triax.stress(triax.wall_front_id))/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), 'right',-triax.stress(triax.wall_right_id),'front',-triax.stress(triax.wall_front_id)

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

newton.damping=0.4
O.saveTmp()

def history():
ev=-triax.strain-triax.strain-triax.strain,
s11=-triax.stress(triax.wall_right_id),
s22=-triax.stress(triax.wall_top_id),
s33=-triax.stress(triax.wall_front_id),
dev=-triax.stress(triax.wall_top_id)-10000,
i=O.iter)

def stopIfDamaged():
if -triax.strain>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

## Question information

Language:
English Edit question
Status:
Solved
For:
Assignee:
No assignee Edit question
Solved by:
Leonard
Solved:
2019-10-07
Last query:
2019-10-07
 Jérôme Duriez (jduriez) said on 2019-09-27: #1

Hi,

Just a small hint: your stress state being not isotropic, you could investigate what may be not isotropic in your script (the packing ? the movements of the surrounding walls ? ... ?)

 Leonard (z2521899293) said on 2019-09-27: #2

Hi Jérôme,
Thanks for your hint, I will dive into it following your suggestion.
If I have solved this problem (or not), I will come back here to share my code.

Cheers,
Leonard

 Leonard (z2521899293) said on 2019-10-07: #3

Hi,
For those who may have the similar questions, one possible way is:
First, using the strategy in , that is the confining pressure is generated by growing particles, to reach a confinement which is roughly close to target confinement. Like the MWE in this question.
Then, set internalCompaction=False, to control the 6 walls moving to reach exact confinement.

Cheers,
Leonard