New question #700213 on Yade:
https://answers.launchpad.net/yade/+question/700213
Hello,
I'm trying to do a uniaxial compression of a periodic box with a 3x3x3 grid of
particles.
I'm implementing periodic boundary conditions for level-set particles. I had
some unexpected results (level set spheres gave a much lower density). Now I
want to do a unit test by comparing a very simple simulation using analytical
spheres and level-set spheres.
Whenever I run the simulation (see code below), the pressure increases to the
target pressure properly. However, hereafter the simulation keeps iterating but
does not initiate the decompression stage defined by compactionFinished(). I
have attempted to change maxUnbalanced and relStressTol to see if the issue is
there but had no success.
Would any of you happen to know why the next stage is not initiated by
PeriTriaxController?
With kind regards,
Danny
###########################################################
# encoding: utf-8
from __future__ import print_function
sigmaIso = -1e6
ydim = 1.0
xdim = 1.0
zdim = 1.0
r = 1.0/3.0/2.0
#import matplotlib
#matplotlib.use('Agg')
# generate loose packing
from yade import pack, qt, plot
O.periodic = True
O.bodies.append(
[
# In xy-plane
sphere(center=(r, r, r), radius=r, fixed=True),
sphere(center=(3*r, r, r), radius=r, fixed=True),
sphere(center=(5*r, r, r), radius=r, fixed=True),
sphere(center=(r, 3*r, r), radius=r, fixed=True),
sphere(center=(3*r, 3*r, r), radius=r, fixed=True),
sphere(center=(5*r, 3*r, r), radius=r, fixed=True),
sphere(center=(r, 5*r, r), radius=r, fixed=True),
sphere(center=(3*r, 5*r, r), radius=r, fixed=True),
sphere(center=(5*r, 5*r, r), radius=r, fixed=True),
# Repeat in z-direction
sphere(center=(r, r, 3*r), radius=r, fixed=True),
sphere(center=(3*r, r, 3*r), radius=r, fixed=True),
sphere(center=(5*r, r, 3*r), radius=r, fixed=True),
sphere(center=(r, 3*r, 3*r), radius=r, fixed=True),
sphere(center=(3*r, 3*r, 3*r), radius=r, fixed=True),
sphere(center=(5*r, 3*r, 3*r), radius=r, fixed=True),
sphere(center=(r, 5*r, 3*r), radius=r, fixed=True),
sphere(center=(3*r, 5*r, 3*r), radius=r, fixed=True),
sphere(center=(5*r, 5*r, 3*r), radius=r, fixed=True),
# Repeat in z-direction
sphere(center=(r, r, 5*r), radius=r, fixed=True),
sphere(center=(3*r, r, 5*r), radius=r, fixed=True),
sphere(center=(5*r, r, 5*r), radius=r, fixed=True),
sphere(center=(r, 3*r, 5*r), radius=r, fixed=True),
sphere(center=(3*r, 3*r, 5*r), radius=r, fixed=True),
sphere(center=(5*r, 3*r, 5*r), radius=r, fixed=True),
sphere(center=(r, 5*r, 5*r), radius=r, fixed=True),
sphere(center=(3*r, 5*r, 5*r), radius=r, fixed=True),
sphere(center=(5*r, 5*r, 5*r), radius=r, fixed=True)
]
)
O.cell.hSize = Matrix3(xdim, 0.0, 0.0, 0.0, ydim, 0.0, 0.0, 0.0, zdim)
O.engines = [
ForceResetter(),
InsertionSortCollider([Bo1_Sphere_Aabb()]),
InteractionLoop([Ig2_Sphere_Sphere_ScGeom()],
[Ip2_FrictMat_FrictMat_FrictPhys()], [Law2_ScGeom_FrictPhys_CundallStrack()]),
PeriTriaxController(
label='triax',
# specify target values and whether they are strains or stresses
goal=(0, 0, sigmaIso),
stressMask=4,
# type of servo-control
dynCell=True,
maxStrainRate=(0, 0, 0.005),
# wait until the unbalanced force goes below this value
maxUnbalanced=.05,
relStressTol=1e-4,
# call this function when goal is reached and the packing is
stable
doneHook='compactionFinished()'
),
NewtonIntegrator(damping=.2),
PyRunner(command='addPlotData()', iterPeriod=100),
]
O.dt = .25 * PWaveTimeStep()
def compactionFinished():
# set the current cell configuration to be the reference one
O.cell.trsf = Matrix3.Identity
# change control type: keep constant confinement in x,y, 20%
compression in z
triax.goal = (0, 0, .2)
triax.stressMask = 3
# allow faster deformation along x,y to better maintain stresses
triax.maxStrainRate = (0, 0, 0.005)
# next time, call triaxFinished instead of compactionFinished
triax.doneHook = 'triaxFinished()'
# do not wait for stabilization before calling triaxFinished
triax.maxUnbalanced = 10
def triaxFinished():
print('Finished')
plot.saveDataTxt('results_ptt3.txt',vars=('compDens','sz'))
O.pause()
totalMass = 0
for b in O.bodies:
totalMass+=b.state.mass
def addPlotData():
plot.addData(
unbalanced=unbalancedForce(),
i=O.iter,
sz=abs(triax.stress[2]),
compDens=totalMass/O.cell.volume
)
# ez=abs(triax.strain[2]),
print("The total mass is", totalMass)
print("The number of particles is",len(O.bodies))
#print("The volume is", O.cell.volume)
# define what to plot
plot.plots = {
'i': ('unbalanced'),
'compDens': ('sz'),
}
# show the plot
plot.plot()
O.run()
--
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
