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. _______________________________________________ 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