yade-users team mailing list archive
-
yade-users team
-
Mailing list archive
-
Message #11780
[Question #271071]: triaxial stress test simulation based on tutorial does not work correctly
New question #271071 on Yade:
https://answers.launchpad.net/yade/+question/271071
Hi All,
I have tried to remodel below paper:
Numerical simulations of triaxial test with sand using DEM, Widuli\'nski, Kozicki, Tejchman, Arch. of Hydroeng. and Env. Mech. \textbf{29} (2009)
Based on triaxial tutorial script. However it does not work and as you can see below the Void ratio and porosity will be "nan"
.
.
.
.
mean stress of engine 8546192948.55
unbalanced force: 0.0108401100247 mean stress: 7400563170.4
void ratio= -0.968143717452 porosity= -30.3909822492
mean stress of engine 7640368975.62
unbalanced force: 0.000639395631987 mean stress: 8796238714.24
void ratio= -0.994549113884 porosity= -182.45641033
mean stress of engine 8640926640.02
unbalanced force: 1.39839137595e-05 mean stress: 7687490942.18
void ratio= -0.998338565076 porosity= -600.889358287
mean stress of engine 7686078969.72
unbalanced force: 1.44304333932 mean stress: 3.04920990888e+12
void ratio= -0.99999998077 porosity= -52002133.7905
mean stress of engine 3.04920991136e+12
unbalanced force: nan mean stress: nan
void ratio= nan porosity= nan
mean stress of engine nan
.
.
.
I just tried to set all of the inputs based on Kozicki paper, to be honest do not have any idea why it does not work. would you please kindly instruct me.
Thanks heaps in advance.
here is the script:
from yade import pack
############################################
### DEFINING VARIABLES AND MATERIALS ###
############################################
# The following 5 lines will be used later for batch execution
nRead=readParamsFromTable(
num_spheres=1000, # number of spheres
compFricDegree = 30, # contact friction during the confining phase
key='_triax_base_', # put you simulation's name here
unknownOk=True,
isoStress=100000, # stress for the isotropic compression phase (1)
conStress=100000 # confinement stress, for the deviatoric loading session
)
from yade.params import table
num_spheres=table.num_spheres # number of spheres
key=table.key
targetPorosity = 0.53 # the porosity we want for the packing
compFricDegree = table.compFricDegree # initial contact friction during the confining phase (will be decreased during the REFD compaction process)
finalFricDegree = 30 # contact friction during the deviatoric loading
rate=-0.1 # loading rate (strain rate)
damp=0.2 # damping coefficient
stabilityThreshold=0.1 # we test unbalancedForce against this value in different loops (see below)
young=30e9 # contact stiffness
mn,mx=Vector3(0,0,0),Vector3(0.1,0.1,0.1) # corners of the initial packing
thick = 0.01 # thickness of the plates
Spoisson=0.3 # poisson ratio of shperes
Wpoisson=0.3 # poisson ratio of Walls
Sdensity=2600 # density of shperes
Wdensity=0 # density of walls
stop_strain= -0.3 # 30%
#kn=7.5e4 # normal stiffness: Kn=2.*young*R1*young*R2/(young*R1+young*R2) ( R1 and R2 are the radii of the two particles in contact)
#ks=2.25e4 # shear stiffness: Ks/Kn=poisson
###################################################
### Create Materials for Spheres and Walls ###
###################################################
O.materials.append(FrictMat(young=young,poisson=Spoisson,frictionAngle=radians(compFricDegree),density=Sdensity,label='spheres'))
O.materials.append(FrictMat(young=young,poisson=Wpoisson,frictionAngle=0,density=Wdensity,label='walls'))
###########################################
### Create Walls Around the Packing ###
###########################################
walls=aabbWalls([mn,mx],thickness=0,material='walls')
wallIds=O.bodies.append(walls)
################################################################################
### use a SpherePack object to generate a random loose particles packing ###
################################################################################
sp=pack.SpherePack()
psdSizes=[0.002,0.003,0.004,0.005,0.006,0.007,0.008] # (sizes or radii of the grains vary from 2mm to 9.5mm)
psdCumm=[0.14,0.28,0.34,0.50,0.65,0.85,1.00] # for the code do not use percentage
sp.makeCloud(mn,mx,-1,0,num_spheres,False, 0.95,psdSizes,psdCumm,False,seed=1) #"seed" make the "random" generation always the same
sp.toSimulation(material='spheres')
############################
### DEFINING ENGINES ###
############################
## TriaxialStressController will be used to control stress and strain. It controls particles size and plates positions.
## this control of boundary conditions was used for instance in http://dx.doi.org/10.1016/j.ijengsci.2008.07.002
triax=TriaxialStressController(
maxMultiplier=1.+2e4/young, # spheres growing factor (fast growth)
finalMaxMultiplier=1.+2e3/young, # spheres growing factor (slow growth)
thickness = 0,
## switch stress/strain control using a bitmask. What is a bitmask, huh?!
## Say x=1 if stess is controlled on x, else x=0. Same for for y and z, which are 1 or 0.
## Then an integer uniquely defining the combination of all these tests is: mask = x*1 + y*2 + z*4
## to put it differently, the mask is the integer whose binary representation is xyz, i.e.
## "100" (1) means "x", "110" (3) means "x and y", "111" (7) means "x and y and z", etc.
stressMask = 7,
internalCompaction=True, # If true the confining pressure is generated by growing particles
)
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()]
),
## 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 #Display spheres with 2 colors for seeing rotations better
if nRead==0: yade.qt.Controller(), yade.qt.View()
#######################################
### APPLYING CONFINING PRESSURE ###
#######################################
#the value of (isotropic) confining stress defines the target stress to be applied in all three directions
triax.goal1=triax.goal2=triax.goal3=-100000
while 1:
O.run(1000, True)
#the global unbalanced force on dynamic bodies, thus excluding boundaries, which are not at equilibrium
unb=unbalancedForce()
#note: triax.stress(k) returns a stress vector, so we need to keep only the normal component
meanS=(triax.stress(triax.wall_right_id)[0]+triax.stress(triax.wall_top_id)[1]+triax.stress(triax.wall_front_id)[2])/3
print 'unbalanced force:',unb,' mean stress: ',meanS
print 'void ratio=',triax.porosity/(1-triax.porosity), 'porosity=', triax.porosity
print 'mean stress of engine', triax.meanStress
if unb<stabilityThreshold and abs(meanS-table.isoStress)/table.isoStress<0.01: #0.001
break
O.save('confinedState'+key+'.xml')
print "### Isotropic state saved ###"
print 'current porosity=',triax.porosity
print 'current void ratio=',triax.porosity/(1-triax.porosity)
print 'step that starts the deviatoric loading ', O.iter
from yade import plot
from pprint import pprint
plot.reset()
e22Check=triax.strain[1]
plot.addData(CheckpointStep=O.iter,Porosity=triax.porosity,e22=triax.strain[1])
plot.saveDataTxt('checkpoint.txt',vars=('CheckpointStep','Porosity','e22'))
O.save('confinedState'+key+'.yade.gz')
print "### Isotropic state saved ###"
###################################################
### REACHING A SPECIFIED POROSITY PRECISELY ###
###################################################
### reach a prescribed value of porosity with the REFD algorithm
### ( http://dx.doi.org/10.2516/ogst/2012032 and
### http://www.geosyntheticssociety.org/Resources/Archive/GI/src/V9I2/GI-V9-N2-Paper1.pdf)
import sys #this is only for the flush() below
while triax.porosity>targetPorosity:
## we decrease friction value and apply it to all the bodies and contacts
compFricDegree = 0.95*compFricDegree
setContactFriction(radians(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(50,1)
O.save('compactedState'+key+'.yade.gz')
print "### Compacted state saved ###"
from pprint import pprint
plot.reset()
e22Check=triax.strain[1]
plot.addData(CheckpointStep=O.iter,Porosity=triax.porosity,e22=triax.strain[1])
plot.saveDataTxt('checkpoint.txt',vars=('CheckpointStep','Porosity','e22'))
##############################
### DEVIATORIC LOADING ###
##############################
##We move to deviatoric loading, let us turn internal compaction off to keep particles sizes constant
triax.internalCompaction=False
## Change contact friction (remember that decreasing it would generate instantaneous instabilities)
setContactFriction(radians(finalFricDegree))
##set stress control on x and z, we will impose strain rate on y
triax.stressMask = 5
##now goal2 is the target strain rate
triax.goal2=rate
## we define the lateral stresses during the test, here the same 10kPa as for the initial confinement.
triax.goal1=-100000
triax.goal3=-100000
newton.damping=0.3 # change damping here ( Just in case)
##Save temporary state in live memory. This state will be reloaded from the interface with the "reload" button.
O.saveTmp()
# Termination condition
def check_stop():
if triax.strain[1] <= stop_strain:
print('Reached target strain, stopping')
O.pause()
#if triax.strain[1] < e22Check # stop condition: the simulation will stop at epsilon = 0.3( az tamrin code)
O.run(50); O.wait()
#####################################################
### Example of how to record and plot data ###
#####################################################
from yade import plot
### a function saving variables
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)
if 1:
## include a periodic engine calling that function in the simulation loop
O.engines=O.engines[0:5]+[PyRunner(iterPeriod=20,command='history()',label='recorder')]+O.engines[5:7]
##O.engines.insert(4,PyRunner(iterPeriod=20,command='history()',label='recorder'))
else:
## With the line above, we are recording some variables twice. We could in fact replace the previous
## TriaxialRecorder
## by our periodic engine. Uncomment the following line:
O.engines[4]=PyRunner(iterPeriod=20,command='history()',label='recorder')
O.run(100,True)
### declare what is to plot. "None" is for separating y and y2 axis
plot.plots={'i':('e11','e22','e33',None,'s11','s22','s33')}
### the traditional triaxial curves would be more like this:
plot.plots={'e22':('s11','s22','s33',None,'ev')}
## display on the screen (doesn't work on VMware image it seems)
plot.plot()
--
You received this question notification because your team yade-users is
an answer contact for Yade.
