# How to reduce unb while keeping state constant

Hi,
I read from the book <Particulate Discrete Element Modelling> in which the author suggests that in servo-controlled simulations once the desired specific stress state is attained, it is good practice to halt all boundary motion and run through a number of DEM calculation cycles to ensure that the system is in equilibrium at that required stress level.

So I want to halt the 6 walls in the triaxial compression test and then run some steps. Here I tied in 2 ways:

The first approach I came up with was fixing all the walls (please see script2-approach 1)

After I O.run(10000), I found that: 1)the unbalancedForce decreased which is expected, 2) triax.porosity is the same which is good. 3) However, the stress on each wall seems to become Vector3(0,0,0), which means the stress state of the sampel is not remained.

The second approach I tried was to keep the stress of all the walls constant as confining pressure (please see script2-approach 2).
After I O.run(10000), I found that: 1) unb decreased which is expected, 2) triax.porosity changed which is not desired because it means the dense state of the sample changed. 3) the stress on each wall get more precisely close to the confining stress which is good.

Both the two ways couldn't keep the current state , do you have any ideas to let simulations run but keep current state?
---
Here is the MWE (modified from Bruno's script) if you want to reproduce this. The MWE is in two-steps: 1st for making sample and 2nd for reloading sample and checking.
######script1#######
from __future__ import division
num_spheres=7000#
targetPorosity = 0.44
compFricDegree = 30
finalFricDegree = 18
rate=-0.01
damp=0.6
stabilityThreshold=0.001
young=2e8
confinement=100e3

mn,mx=Vector3(0,0,0),Vector3(0.07,0.14,0.07)

MatWall=O.materials.append(FrictMat(young=young,poisson=0.5,frictionAngle=0,density=0,label='walls'))
walls=aabbWalls([mn,mx],thickness=0,material='walls')
wallIds=O.bodies.append(walls)

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

triax=TriaxialStressController(
maxMultiplier=1.+3e7/young,
finalMaxMultiplier=1.+16e4/young,
thickness = 0,
internalCompaction=True,
)
newton=NewtonIntegrator(damping=damp)

O.engines=[
ForceResetter(),
InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
InteractionLoop(
[Ig2_Sphere_Sphere_ScGeom6D(),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'),
newton,
]

Gl1_Sphere.stripes=0
triax.goal1=triax.goal2=triax.goal3=-confinement

while 1:
O.run(1000, True)
#the global unbalanced force on dynamic bodies, thus excluding boundaries, which are not at equilibrium
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.0001:
break

print "### Isotropic state saved ###"

import sys
while triax.porosity>targetPorosity:
compFricDegree = 0.95*compFricDegree
print "\r Friction: ",compFricDegree," porosity:",triax.porosity,
sys.stdout.flush()
O.run(500,1)

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

######script2#######

confinement=100e3
triax=TriaxialStressController(
thickness = 0,
internalCompaction=False,
)
newton=NewtonIntegrator(damping=0.4)

O.engines=[
ForceResetter(),
InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
InteractionLoop(
[Ig2_Sphere_Sphere_ScGeom6D(),Ig2_Box_Sphere_ScGeom()],
[Ip2_FrictMat_FrictMat_FrictPhys()],
[Law2_ScGeom6D_CohFrictPhys_CohesionMoment(useIncrementalForm=True,always_use_moment_law=True),Law2_ScGeom_FrictPhys_CundallStrack()]
),
GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
triax,
TriaxialStateRecorder(iterPeriod=1000,file='WallStresses'),
newton,
]
O.step()

## Try approach 1
triax.wall_top_activated=False
triax.wall_bottom_activated=False
triax.wall_front_activated=False
triax.wall_back_activated=False
triax.wall_left_activated=False
triax.wall_right_activated=False

## Try approach 2
# triax.goal1=triax.goal2=triax.goal3=-confinement

print "before run:",triax.stress(triax.wall_front_id),triax.stress(triax.wall_top_id),triax.stress(triax.wall_left_id)
print "porosity:",triax.porosity,"unb:",unbalancedForce()
O.run(10000,True)
print "after run:",triax.stress(triax.wall_front_id),triax.stress(triax.wall_top_id),triax.stress(triax.wall_left_id)
print "porosity:",triax.porosity,"unb:",unbalancedForce()

Thanks
Leonard

## Question information

Language:
English Edit question
Status:
Solved
For:
Assignee:
No assignee Edit question
Solved by:
Jan Stránský
Solved:
2020-08-22
Last query:
2020-08-22
2020-08-22 Jan Stránský (honzik) said on 2020-08-22: #1

Hello,

I think what you require is impossible.

All the three quantities (stress, unbalanced force, porosity) are not independent.
fixed porosity: higher unbalanced force -> higher stress
fixed unbalanced force: lower porosity -> higher stress
fixed stress: higher unbalanced force -> lower porosity
So you cannot fix both stress and porosity and let unbalanced force disappear...

cheers
Jan

 Leonard (z2521899293) said on 2020-08-22: #2

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