Question #698594 on Yade changed:
https://answers.launchpad.net/yade/+question/698594
Status: Answered => Open
Leonard is still having a problem:
Hello Jan and Luc,
Thanks very much for your reply.
Here I made two MWE for rigid wall case and periodic boundary case
respectively. It may take around 8 min for each of the simulation, or
you can also find the results of stress-strain and volumetric strain
evolution at [1].
>same strain rate (I think)
Because in rigid wall which uses TriaxialStressController, I can specify the
loading rate by triax.goal. But as far as I understand, it seems that for
PeriTriaxController, it can only specify the loading rate by maxStrainRate. So
I made triax.strain=1 (rigid case) and triax.maxStrainRate = (1., 1., 1.) in
periodic case. Do you know how to specify the strain rate in
PeriTriaxController?
>the two cases use ... the same of sample size
Now I know they are not the same. Just explain here why I though they were same
before. If you run the periodic MWE below, it prints out the cell size just
before the deviatoric loading, this size is almost the same as the size in
rigid case.
>initial porostiy
Sorry for not clearly describe it. Yes, the initial porosity I mentioned is the
porosity at the begining of the deviatoric loading. In the following two MWE,
they are both 0.358 getting by function of porosity(). But note that, in rigid
boundary case, I also print the porosity at the same state by using
triax.porosity, which retures 0.33.
>coordination number
Thanks Luc, that is a good point to mention. In the two MWE, I print the
coordination number (CN) just before the deviatoric loading. And the CN of
rigid boundary is 7.28, while the CN in periodic is 6.3. Which means there are
more contacts in rigid wall case, indicating a denser state in general. Then my
question is why they have same inputtings (porosity, particle number...) but
they have different CN? And another question is how can I make the same sample
under the two boundary conditions? It would be conflict if I control them as
the same from porosity and from CN.
################################### rigidWall case ##########################
from __future__ import division
from yade import pack, plot
import numpy as np
num_spheres=7000# number of spheres
targetPorosity = 0.33
compFricDegree = 30
finalFricDegree = 35
damp=0.6
stabilityThreshold=0.001
confinement=100e3
mn,mx=Vector3(0,0,0),Vector3(0.07,0.14,0.07) #determine the size of the sample
young=5e6
## create materials for spheres and plates
MatWall=O.materials.append(FrictMat(young=young,poisson=0.3,frictionAngle=0,density=0,label='walls'))
MatSand =
O.materials.append(FrictMat(young=young,poisson=0.3,frictionAngle=radians(30),\
density=2650.0,label='sand'))
walls=aabbWalls([mn,mx],thickness=0,material='walls')
wallIds=O.bodies.append(walls)
sp=pack.SpherePack()
sp.makeCloud(mn,mx,-1,0,num_spheres,False, 0.95,seed=1)
O.bodies.append([sphere(center,rad,material='sand') for center,rad in sp])
Gl1_Sphere.quality=3
triax=TriaxialStressController(
maxMultiplier=1.+2e4/young,
finalMaxMultiplier=1.+2e3/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,
PyRunner(iterPeriod=100,command='history()',label='recorder'),
PyRunner(command='stopIfDamaged()',iterPeriod=200,label="endSimulation"),
]
recorder.dead=True
endSimulation.dead=True
Gl1_Sphere.stripes=0
triax.goal1=triax.goal2=triax.goal3=-confinement
def getCN():
# get coordination number
CN=0
for i in O.bodies:
if isinstance(i.shape,Sphere):
CN +=len(i.intrs())
Z = CN/ 7000.0
return Z
while 1:
O.run(1000, True)
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 "### state 1 Isotropic completed ###"
import sys
while triax.porosity - targetPorosity>0.0001:
compFricDegree = 0.95*compFricDegree
setContactFriction(radians(compFricDegree))
print "\r Friction: ",compFricDegree," porosity:",triax.porosity,
sys.stdout.flush()
O.run(500,1)
print "### state 2 Reach target porosity completed ###"
print "particle radii is: ", O.bodies[2000].shape.radius ## here get the
particle radius which will be used as input in makecloud in periodic boundary
case.
print 'triax.porosity at 100 kPa',triax.porosity,'porosity()',porosity()
print ("coordination number at 100 kPa is",getCN())
topStress_initial=-triax.stress(triax.wall_top_id)[1]
bottomStress_initial=-triax.stress(triax.wall_bottom_id)[1]
# start deviatoric loading
triax.internalCompaction=False
triax.wall_bottom_activated=False
recorder.dead=False
endSimulation.dead=False
triax.stressMask=5
triax.goal2=-1
triax.goal1=-confinement
triax.goal3=-confinement
setContactFriction(radians(finalFricDegree))
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],
q=-triax.stress(triax.wall_top_id)[1]-100000,)
def stopIfDamaged():
if -triax.strain[1]>0.4:
O.pause()
plot.saveDataTxt('results_rigidWall.txt')
plot.plots={'e22':('q',),'e22 ':('ev')}
plot.plot()
###################################### end of the 1st MWE #############
#################### periodic case #######################3
from __future__ import print_function
sigmaIso = -100e3
from yade import pack, qt, plot
O.periodic = True
O.cell.setBox(.09,.18,.09)
young=5e6
compFricDegree=23 # using this value can make the porosity at 100 kpa be 0.358,
which is the same as that of rigid wall
finalFricDegree = 35
MatSand =
O.materials.append(FrictMat(young=young,poisson=0.3,frictionAngle=radians(30),density=2650.0,label='sand'))
sp = pack.SpherePack()
mn,mx=Vector3(0,0,0),Vector3(0.1,0.2,0.1)
sp.makeCloud(mn,mx,0.0025,0,7000,False, 0.95,seed=1)
sp.toSimulation(material="sand")
# here set contact friction angle to reach a low porosity at target confining
pressure.
setContactFriction(radians(compFricDegree))
O.engines = [
ForceResetter(),
InsertionSortCollider([Bo1_Sphere_Aabb()]),
InteractionLoop([Ig2_Sphere_Sphere_ScGeom()],
[Ip2_FrictMat_FrictMat_FrictPhys()], [Law2_ScGeom_FrictPhys_CundallStrack()]),
PeriTriaxController(
label='triax',
goal=(sigmaIso, sigmaIso, sigmaIso),
stressMask=7,
dynCell=True,
maxStrainRate=(10, 10, 10),
maxUnbalanced=.1,
relStressTol=1e-3,
doneHook='compactionFinished()'
),
NewtonIntegrator(damping=.6),
PyRunner(command='addPlotData()', iterPeriod=100,label="recorder"),
]
O.dt = .5 * PWaveTimeStep()
recorder.dead=True
def addPlotData():
plot.addData(
unbalanced=unbalancedForce(),
i=O.iter,
sxx=-triax.stress[0],
syy=-triax.stress[1],q=-triax.stress[1]-100000,
szz=-triax.stress[2],
exx=-triax.strain[0],ev=triax.strain[0]+triax.strain[1]+triax.strain[2],
eyy=-triax.strain[1],
ezz=-triax.strain[2],
)
O.trackEnergy = True
plot.plots = {
'eyy': ('q',),
'eyy ': ('ev'),}
plot.plot()
def compactionFinished():
print ("porosity at 100 kPa is",porosity())
print ("coordination number at 100 kPa is",getCN())
print ("O.cell.size is", O.cell.size)
setContactFriction(radians(finalFricDegree))
O.pause()
O.cell.trsf = Matrix3.Identity
triax.goal = (sigmaIso,-0.4, sigmaIso)
triax.stressMask = 5
triax.maxStrainRate = (1., 1., 1.)
triax.doneHook = 'triaxFinished()'
triax.maxUnbalanced = 10
O.step() # if I didnot use O.step(), the initial eyy will not start
from 0.
recorder.dead=False
def triaxFinished():
print('Finished')
O.pause()
plot.saveDataTxt('results_periodic.txt')
def getCN():
CN=0
for i in O.bodies:
if isinstance(i.shape,Sphere):
CN +=len(i.intrs())
Z = CN/ 7000.0
return Z
##################################################
Thanks
Leonard
[1]https://we.tl/t-xlkvIPB9zx
ps: both of the MWE are modified from YADE trunk.
--
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
