Question #482413 on Yade changed:
https://answers.launchpad.net/yade/+question/482413
Status: Answered => Open
Wei Cao is still having a problem:
Hi, thanks again. I modified my code to make sure that the sample is at
equilibrium before applying the deviatic stress, by adding:
{
while 1:
O.run(1000, True)
##the global unbalanced force on dynamic bodies, thus excluding boundaries,
which are not at equilibrium
unb=unbalancedForce()
print 'unbalanced force:',unb,' mean stress: ',triax.meanStress
if unb<stabilityThreshold and abs(-100000-triax.meanStress)/100000<0.001:
break
}
Could you helpt me to look at the code? I only want to know why the strain
control give me so strange results. It really bothers me for a long time.
Thank you！
Here is my code:
##############################################################################################################################
#Authors: Luc Sibille luc.sibi...@3sr-grenoble.fr
# Franck Lomine franck.lom...@insa-rennes.fr
#
# Script to simulate buoyancy in a granular assembly with the DEM-LBM
coupling.
# This script was written for a practical work in the OZ ALERT school in
Grenoble 2011.
# Script updated in 2014
##############################################################################################################################
from yade import pack, timing, utils, qt, plot
############################################
### DEFINING VARIABLES AND MATERIALS ###
############################################
# The following 5 lines will be used later for batch execution
nRead=readParamsFromTable(
num_spheres=10000,# number of spheres
compFricDegree = 35, # contact friction during the confining phase
key='_triax_base_', # put you simulation's name here
unknownOk=True
)
from yade.params import table
num_spheres=table.num_spheres# number of spheres
key=table.key
# targetPorosity = 0.5 #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 = 35 # contact friction during the deviatoric loading
rate=-0.01 # loading rate (strain rate)
damp=0.7 # damping coefficient
stabilityThreshold=0.01 # we test unbalancedForce against this value in
different loops (see below)
young=5e8 # contact stiffness
mn,mx=Vector3(0,0,0),Vector3(1,1,1) # corners of the initial packing
## create materials for spheres
O.materials.append(FrictMat(young=young,poisson=0.3,frictionAngle=radians(compFricDegree),density=2600,label='spheres'))
## control to open the file
Not_Open_flag = 1
if Not_Open_flag:
stress_info =
open("/home/caowei/Documents/myYade/trunk/examples/stress_info.txt",'a')
contact_info =
open("/home/caowei/Documents/myYade/trunk/examples/contactOri.txt",'a')
particle_info =
open("/home/caowei/Documents/myYade/trunk/examples/particle.txt",'a')
Not_Open_flag = 0
sigmaIso = -1e5
rMean = 0.035 #mean radius (m) of solid particles
Gravity=-0.1 #gravity (m.s-2) that will be applied to solid particles
#definition of the simulation domain (position of walls), the domain is a 3D
domain but particles will be generated on a unique plane at z=0, be carefull
that wall normal to z directions do not touch particles.
lowerCornerW = (0.001,0.001,0.001)
upperCornerW = (0.701,0.701,0.701)
#definitions of the domain for particles: positions in z direction is set to 0
to force a 2D granular assembly
#lowerCornerS = (0.00001,0.00001,0)
#upperCornerS = (0.00701,0.00501,0)
lowerCornerS = (0.001,0.001,0.001)
upperCornerS = (0.701,0.701,0.701)
nbSpheres = 5000 #this not the real number of particles but number of
times where we try to insert a new particle, the real number of
particles can be much lower!
def addPlotData():
plot.addData(unbalanced=unbalancedForce(),i=O.iter,
sxx=triax.stress(triax.wall_right_id)[0],syy=triax.stress(triax.wall_top_id)[1],szz=triax.stress(triax.wall_front_id)[2],
exx=triax.strain[0],eyy=triax.strain[1],ezz=triax.strain[2],
# save all available energy data
Etot=O.energy.total(),**O.energy
)
stress_info.write("%6.2f " %(triax.stress(triax.wall_right_id)[0]))
stress_info.write("%6.2f " %(triax.stress(triax.wall_top_id)[1]))
stress_info.write("%6.2f " %(triax.stress(triax.wall_front_id)[2]))
#print("stress:")
#print(triax.stress(triax.wall_right_id)[0])
#print(triax.stress(triax.wall_top_id)[1])
#print(triax.stress(triax.wall_front_id)[2])
#print("strain:")
#print(triax.strain[0])
#print(triax.strain[1])
#print(triax.strain[2])
for i in range(3):
stress_info.write("%6.2f " %(triax.strain[i]))
stress_info.write("\n")
if O.iter % 20000 == 0:
for i in O.interactions:
contact_info.write("%d " %(i.id1))
contact_info.write("%d " %(i.id2))
contact_info.write("%.4f " %(i.phys.normalForce[0]))
contact_info.write("%.4f " %(i.phys.normalForce[1]))
contact_info.write("%.4f " %(i.phys.normalForce[2]))
contact_info.write("%.4f " %(i.phys.shearForce[0]))
contact_info.write("%.4f " %(i.phys.shearForce[1]))
contact_info.write("%.4f " %(i.phys.shearForce[2]))
contact_info.write("\n")
for i in range(len(O.bodies)):
if i > 5:
particle_info.write("%d " %(i))
if O.bodies[i].isClump:
particle_info.write("\n")
else:
particle_info.write("%f " %(O.bodies[i].state.pos[0]))
particle_info.write("%f " %(O.bodies[i].state.pos[1]))
particle_info.write("%f " %(O.bodies[i].state.pos[2]))
particle_info.write("%f " %(O.bodies[i].shape.radius))
particle_info.write("%d " %(O.bodies[i].clumpId))
particle_info.write("\n")
# enable energy tracking in the code
O.trackEnergy=True
# define what to plot
plot.plots={'i':('unbalanced',),'i ':('sxx','syy','szz'),'
i':('exx','eyy','ezz'),
# energy plot
' i ':(O.energy.keys,None,'Etot'),
}
# show the plot
plot.plot()
def compactionFinished():
# set the current cell configuration to be the reference one
# O.cell.trsf=Matrix3.Identity
triax.depth0 = triax.depth
triax.height0 = triax.height
triax.width0 = triax.width
triax.strain[0] = triax.strain[1] = triax.strain[2] = 0
# change control type: keep constant confinement in x,y, 20% compression in z
triax.stressMask = 0
triax.goal1 = 0.001
triax.goal2 = 0.001
triax.goal3 = -0.002
count = 0
countN = 2
while count <= countN:
for i in range(100):
O.step()
if triax.stress(triax.wall_front_id)[2] < -1.1e5:
if triax.goal3 < 0:
count = count + 1
triax.goal1 = -0.001
triax.goal2 = -0.001
triax.goal3 = 0.002
if triax.stress(triax.wall_front_id)[2] > -0.9e5:
if triax.goal3 > 0:
count = count + 1
triax.goal1 = 0.001
triax.goal2 = 0.001
triax.goal3 = -0.002
stress_info.close()
contact_info.close()
particle_info.close()
print 'Finished'
O.pause()
# allow faster deformation along x,y to better maintain stresses
# triax.maxStrainRate=(0.05,0.05,.1)
# next time, call triaxFinished instead of compactionFinished
# triax.doneHook = 'triaxFinished()'
# do not wait for stabilization before calling triaxFinished
# triax.maxUnbalanced =1 0
# O.run(10000, True)
#def triaxFinished():
#stress_info.close()
#contact_info.close()
#particle_info.close()
#print 'Finished'
#O.pause()
newton=NewtonIntegrator(damping=damp)
triax = TriaxialStressController(
## 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
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=False, # If true the confining pressure is generated by
growing particles
# call this function when goal is reached and the packing is stable
)
## Build of the engine vector
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()]
),
triax,
newton,
## 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),
##
NewtonIntegrator(damping=0.2,gravity=(Gravity,0.,0.),label="ENewton"),
PyRunner(command='addPlotData()',iterPeriod=100),
]
###############################################################################################################
# Creation of 6 walls at the limit of the simulation domain
# defintiion of the material for walls
O.materials.append(FrictMat(young=50e6,poisson=.5,frictionAngle=0.0,density=3000,label='walls'))
#lowerCornerW = (0.001,0.001,0.001)
#upperCornerW = (0.701,0.701,0.701)
mn,mx=Vector3(-0.1,-0.1,-0.1),Vector3(0.8,0.8,0.8) # corners of the
initial packing
## create materials for spheres and plates
O.materials.append(FrictMat(young=young,poisson=0.3,frictionAngle=0.0,density=0,label='wallsss'))
## create walls around the packing
walls=aabbWalls([mn,mx],thickness=0,material='wallsss')
wallIds=O.bodies.append(walls)
###############################################################################################################
# Creation of the assembly of particles
# defintiion of the material for particles
O.materials.append(FrictMat(young=50e6,poisson=.5,frictionAngle=35*3.1415926535/180,density=3000,label='spheres'))
## use a SpherePack object to generate a random loose particles packing
sp=pack.SpherePack()
#sp.makeCloud(lowerCornerS,upperCornerS,-1,0.33,nbSpheres,False, 0.5,seed=1)
#"seed" make the "random" generation always the same
#sp.makeCloud(lowerCornerS,upperCornerS,-1,0.33,nbSpheres,False, 0.5,seed=1)
#sp.makeCloud(lowerCornerS,upperCornerS,rMean=.03,seed=1)
c1=pack.SpherePack([((0,0,0),.03),((.03,0,0),.03)])
sp.makeClumpCloud(lowerCornerS,upperCornerS,[c1],periodic=True,num=5000)
sp.toSimulation(material='spheres')
O.dt=.5*PWaveTimeStep()
# for the samples of this size (diametre is around 1.8cm), the O.dt is 8.8e-5
s, which is about a hundred times of
# the small sample, and this leads to much faster calculation speed.
#Display spheres with 2 colors for seeing rotations better
Gl1_Sphere.stripes=0
if nRead==0: yade.qt.Controller()
#######################################
### 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 = sigmaIso
print"___________________"
print"PHASE 1"
print "Settlement of particle under gravity only, fluid flow d"
triax.goal1 = triax.goal2 = triax.goal3 = sigmaIso / 10 * 1
O.run(2000, True)
triax.goal1 = triax.goal2 = triax.goal3 = sigmaIso / 10 * 2
O.run(2000, True)
triax.goal1 = triax.goal2 = triax.goal3 = sigmaIso / 10 * 4
O.run(2000, True)
triax.goal1 = triax.goal2 = triax.goal3 = sigmaIso / 10 * 7
O.run(2000, True)
triax.goal1 = triax.goal2 = triax.goal3 = sigmaIso / 10 * 10
O.run(2000, True)
stabilityThreshold = 0.001
while 1:
O.run(1000, True)
##the global unbalanced force on dynamic bodies, thus excluding boundaries,
which are not at equilibrium
unb=unbalancedForce()
print 'unbalanced force:',unb,' mean stress: ',triax.meanStress
if unb<stabilityThreshold and abs(-100000-triax.meanStress)/100000<0.001:
break
print(O.iter)
#compactionFinished()
triax.depth0 = triax.depth
triax.height0 = triax.height
triax.width0 = triax.width
# triax.strain[0] = triax.strain[1] = triax.strain[2] = 0
# change control type: keep constant confinement in x,y, 20% compression in z
triax.stressMask = 0
triax.goal1 = 0.0001
triax.goal2 = 0.0001
triax.goal3 = -0.0002
--
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 : yade-users@lists.launchpad.net
Unsubscribe : https://launchpad.net/~yade-users
More help : https://help.launchpad.net/ListHelp
