New question #692524 on Yade:
https://answers.launchpad.net/yade/+question/692524
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
from yade import pack, plot
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'))
MatSand =
O.materials.append(FrictMat(young=young,poisson=0.5,frictionAngle=radians(compFricDegree),density=2600,label='spheres'))
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)
O.bodies.append([sphere(center,rad,material='spheres') for center,rad in sp])
triax=TriaxialStressController(
maxMultiplier=1.+3e7/young,
finalMaxMultiplier=1.+16e4/young,
thickness = 0,
stressMask = 7,
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)[1]
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
setContactFriction(radians(compFricDegree))
print "\r Friction: ",compFricDegree," porosity:",triax.porosity,
sys.stdout.flush()
O.run(500,1)
O.save('sample.yade.gz')
print "### Compacted state saved ###"
######script2#######
from yade import pack, plot
O.load('sample.yade.gz')
confinement=100e3
triax=TriaxialStressController(
thickness = 0,
stressMask = 5,
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,
]
triax.stressMask=7
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
--
You received this question notification because your team yade-users is
an answer contact for Yade.
_______________________________________________
Mailing list: https://launchpad.net/~yade-users
Post to : [email protected]
Unsubscribe : https://launchpad.net/~yade-users
More help : https://help.launchpad.net/ListHelp
