Question #703169 on Yade changed:
https://answers.launchpad.net/yade/+question/703169
Karol Brzezinski posted a new comment:
Hi Xu,
Thank you for your MWE. It looks like it is not what I thought at the
beginning. I don't know what is the reason.
My last clue is that 'utils.getStress()' may give slightly different results
for clumps and spheres. I modified oedometric test example [1] to check this.
But the results for different simulations are scattered.
################################################
from yade import pack, plot
useClumps = True
readParamsFromTable(rMean=.075, rRelFuzz=.3, maxLoad=1e6)
from yade.params.table import *
# create box with free top, and ceate loose packing inside the box
O.bodies.append(geom.facetBox((.5, .5, 1), (.5, .5, 1), wallMask=31))
sp = pack.SpherePack()
sp.makeCloud((0, 0, 0), (1, 1, 2), rMean=rMean, rRelFuzz=rRelFuzz)
sp.toSimulation()
# make material frictionless to avoid force acting on vert walls
O.materials[0].frictionAngle = 0
# I just make spheres smaller to avoid overlapping after replacing by clumps
for b in O.bodies:
if isinstance(b.shape, Sphere):
b.shape.radius*=0.8
if useClumps:
relRadList2 = [1,1]
relPosList2 = [[0.5,0,0],[-0.5,0,0]]
templates= []
templates.append(clumpTemplate(relRadii=relRadList2,relPositions=relPosList2))
O.bodies.replaceByClumps(templates,[1.],discretization=10)
O.engines = [
ForceResetter(),
InsertionSortCollider([Bo1_Sphere_Aabb(), Bo1_Facet_Aabb(),
Bo1_Wall_Aabb()]),
InteractionLoop(
[Ig2_Sphere_Sphere_ScGeom(), Ig2_Facet_Sphere_ScGeom(),
Ig2_Wall_Sphere_ScGeom()],
[Ip2_FrictMat_FrictMat_FrictPhys()],
[Law2_ScGeom_FrictPhys_CundallStrack()]
),
NewtonIntegrator(gravity=(0, 0, 0), damping=0.1),
PyRunner(command='checkUnbalanced()', virtPeriod=1, label='checker'),
]
O.dt = .5 * PWaveTimeStep()
def checkUnbalanced():
# at the very start, unbalanced force can be low as there is only few
contacts, but it does not mean the packing is stable
if O.iter < 5000:
return
# the rest will be run only if unbalanced is < .1 (stabilized packing)
if unbalancedForce() > .1:
return
# add plate at the position on the top of the packing
# the maximum finds the z-coordinate of the top of the topmost particle
O.bodies.append(wall(max([b.state.pos[2] + b.shape.radius for b in
O.bodies if isinstance(b.shape, Sphere)]), axis=2, sense=-1))
global plate # without this line, the plate variable would only exist
inside this function
plate = O.bodies[-1] # the last particles is the plate
# Wall objects are "fixed" by default, i.e. not subject to forces
# prescribing a velocity will therefore make it move at constant
velocity (downwards)
plate.state.vel = (0, 0, -.01)
# start plotting the data now, it was not interesting before
O.engines = O.engines + [PyRunner(command='addPlotData()',
iterPeriod=200)]
# next time, do not call this function anymore, but the next one
(unloadPlate) instead
checker.command = 'unloadPlate()'
def unloadPlate():
# if the force on plate exceeds maximum load, start unloading
if abs(O.forces.f(plate.id)[2]) > maxLoad:
plate.state.vel *= -1
O.pause()
print('Stress underestimated by the factor
{:.2f}'.format(max(plot.data['Fz'])/max(plot.data['szz'])))
def addPlotData():
if not isinstance(O.bodies[-1].shape, Wall):
plot.addData()
return
Fz = O.forces.f(plate.id)[2]
szz=-1*utils.getStress()[2,2]
plot.addData(Fz=Fz,szz=szz, w=plate.state.pos[2] -
plate.state.refPos[2], unbalanced=unbalancedForce(), i=O.iter)
# besides unbalanced force evolution, also plot the displacement-force diagram
plot.plots = {'w': ('Fz','szz')}
plot.plot()
O.run()
################################################
Best wishes,
Karol
[1] https://yade-dem.org/doc/tutorial-examples.html#oedometric-test
PS Jan, sorry for suggesting this 'interaction inside clump' issue, but
I can swear I encountered this some time ago.
--
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
