yade-users team mailing list archive

[Bruno Chareyre <question703169@xxxxxxxxxxxxxxxxxxxxx>]

Question #703169 on Yade changed:
https://answers.launchpad.net/yade/+question/703169

    Status: Open => Answered

Bruno Chareyre proposed the following answer:
It will be fixed with https://gitlab.com/yade-dev/trunk/-/merge_requests/891. 
I confirmed the fix with a modified version of #10.

Note that this example with walls has another issue (which does not
affect periodic BCs): the volume computed by default in getStress is
overestimated because it includes enlarged bounding boxes. When  passing
tha actual volume as argument, like below, the error falls to 1% with
the new version (still high, but that's because equilibrium is only
approximately satisfied).

Bruno

################################################
from yade import pack, plot

useClumps = True

readParamsFromTable(rMean=.075, rRelFuzz=.3, maxLoad=1e5)

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 = [0.8,0.8,0.8,0.8]
 relPosList2 = [[0.6,0,0],[-0.6,0,0],[0,0.7,0],[0,0.25,0.7]]
 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.2),
        PyRunner(command='checkUnbalanced()', iterPeriod=20, label='checker'),
]
O.dt = .5 * PWaveTimeStep()

def checkUnbalanced():
 # 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, -.3)
 # start plotting the data now, it was not interesting before
 O.engines = O.engines + [PyRunner(command='addPlotData()', iterPeriod=20)]
 # 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 = (0,0,0)
 elif plate.state.vel[2]==0 and unbalancedForce() < 0.01:
  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(plate.state.pos[2])[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()

################################################

-- 
You received this question notification because your team yade-users is
an answer contact for Yade.