Energy conservation of clump

Hello,

please provide a MWE [2]

> dW=sijdeij
> Etotal=Estrain+Efriction+Edamp-(Estrain_ini+Efriction_ini+Edamp_ini)

please provide how you compute / get individual components (should be clear form the MWE)

> So I check the script [1]

it is not script, it is source code. Just a terminology note.

> And i have another question

in general, please open a new question for a new question ([2], point 5).
Here it is the same topic, just have it in mind.

> If my opinion is right, the calculation of Kn,Ks is wrong for the clump, the calculation of the normal stress and shear stress in the interaction are also wrong.

It depends on the definition of "wrong".
The stiffnesses are computed for sphere-sphere interaction. Regardless if it is clump or not.
It is implementation / design decision. It can be done differently, but it is done this way.
You can adjust stiffness of individual spheres to match your requirements (e.g. based on clump size).
But this problem should have no effect on energy conservation.

Cheers
Jan

[2] https://www.yade-dem.org/wiki/Howtoask

Hi,

I endorse all the comments from Jan. I can add one more. It is very likely that there are some interactions between clump members. Those interactions do not affect the simulation. However, they may contribute to the potential energy if you are taking them into account.

Cheers,
Karol

Hi Jan,
Thank you of your response.

> please provide how you compute / get individual components
I compute the input work in the script with these lines:

dw=(0.5*(O.cell.preVelGrad+O.cell.velGrad)*utils.getStress())*O.dt*(O.cell.hSize).determinant()
W=W+dw

And the energy calculated from microscalle:

Law=Law2_ScGeom_FrictPhys_CundallStrack(traceEnergy=True)
Estrain=Law.elasticEnergy()
Efriction=Law.plasticDissipation()
Edamp=O.energy['nonviscDamp']
Ekin=utils.kineticEnergy()

> it is not script, it is source code. Just a terminology note.
Thank you for your notification.

> please open a new question for a new question ([2], point 5).
Ok, I'll open a new question later. Thank you.

Cheers,
Xu

Hi Karol,
Thank you for your comment.

> It is very likely that there are some interactions between clump members. Those interactions do not affect the simulation. However, they may contribute to the potential energy if you are taking them into account.

Sorry, I dont know how to distingiosh this part of interaction between clump members. Could you give a little more description of on this?

Cheers,
Xu

Hi,

If you provided MWE, I could show what I mean based on your example.

The simplest approach is to create a single clump of two overlapping spheres (no other bodies, no gravity). Run one iteration of the simulation. What is the energy state in such a simulation?

Cheers,
Karol

> It is very likely that there are some interactions between clump members.

if the clump is created before interactions are created, it should not happen [3]
But it depends on actual use case, so please provide a MWE.

Cheers
Jan

[3] https://gitlab.com/yade-dev/trunk/-/blob/master/pkg/common/Collider.cpp#L32

Hi Karol,
Thank you for your reply.

This is the script i use:

 ######Consolidation#######
from yade import pack,plot,qt,export,ymport,utils
import numpy as np
IsoSigma = -1.e5
O.periodic=True
O.trackEnergy=True
idSand=O.materials.append(FrictMat(young=2e8,poisson=.8,frictionAngle=0.09,density=2650000,label='spheres'))
sp1=pack.SpherePack()
sp1.makeCloud(maxCorner=(0.02, 0.02, 0.02), psdSizes=[0.00017, 0.000191, 0.0002285, 0.00026, 0.000292, 0.000325, 0.00035], psdCumm=[0.0, 0.1, 0.3, 0.5, 0.6, 0.9, 1], periodic=True,num=5000,seed=1)
sp1.toSimulation(color=(0,0,1),material=idSand)
v=utils.getSpheresVolume()
print utils.getSpheresVolume()
fout = file('Packing_volume.txt','w')
fout.write(str(v))
fout.close()

#### show how to use makeClumpTemplate():
relRadList2 = [1,1]
relPosList2 = [[0.4,0,0],[-0.4,0,0]]

templates= []
templates.append(clumpTemplate(relRadii=relRadList2,relPositions=relPosList2))
O.bodies.replaceByClumps(templates,[1.],discretization=10)

law=Law2_ScGeom_FrictPhys_CundallStrack(traceEnergy=True)
O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb()]),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom()],
     [Ip2_FrictMat_FrictMat_FrictPhys()],
     [law]),
# PyRunner(command='fabric()',iterPeriod=10000),
        GlobalStiffnessTimeStepper(),
 NewtonIntegrator(damping=0.20),
 PeriTriaxController( goal=(IsoSigma,IsoSigma,IsoSigma), # Vector6 of prescribed final values
    stressMask=7,
    dynCell=True,
    maxStrainRate=(5.e+0,5.e+0,5.e+0),
    maxUnbalanced=0.0001,
    relStressTol=1.e-3,
    doneHook='Finished()',
    label='p3d'
 ),
# PyRunner(command='plotAddData()',iterPeriod=100),
]

def Finished():
 O.save('post_iso.xml')
 print 'Finished'
 print (O.cell.size[0]*O.cell.size[1]*O.cell.size[2]-v)/v
        print getStress()
        print unbalancedForce()
 O.pause()

O.run()

########Undrained triaxial test#########
from yade import pack,plot,qt,export
from math import fabs
import numpy as np
O.load('post_iso.xml')
O.trackEnergy=True
setContactFriction(0.5)
filename1='Packing_volume.txt'
a=np.loadtxt(filename1)
Ep_ini=O.energy['elastPotential']
Edamp_ini=O.energy['nonviscDamp']
Ekin_ini=utils.kineticEnergy()
work=0

law=Law2_ScGeom_FrictPhys_CundallStrack(traceEnergy=True)
O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb()]),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom()],
     [Ip2_FrictMat_FrictMat_FrictPhys()],
     [law]),
 NewtonIntegrator(damping=0.2),
        GlobalStiffnessTimeStepper(),
 PyRunner(command='inputwork()',iterPeriod=1),
 PyRunner(command='plotAddData()',iterPeriod=10),
]

def inputwork():
 global work
 dwork=(0.5*(O.cell.prevVelGrad+O.cell.velGrad)*utils.getStress()).trace()*O.dt*(O.cell.hSize).determinant()
 work=work+dwork

def plotAddData():
 global Ep_ini
 global Edamp_ini
 global Ekin_ini
        global work
 plot.addData(
  iter=O.iter,iter_=O.iter,
  sxx=utils.getStress()[0,0],
  syy=utils.getStress()[1,1],
  szz=utils.getStress()[2,2],
  exx=O.cell.getSmallStrain()[0,0],
  eyy=O.cell.getSmallStrain()[1,1],
  ezz=O.cell.getSmallStrain()[2,2],
  exx1=O.cell.size[0],
  eyy1=O.cell.size[1],
  ezz1=O.cell.size[2],
  ex1=O.cell.hSize[0,0],
  ey1=O.cell.hSize[1,1],
  ez1=O.cell.hSize[2,2],
  exy1=O.cell.hSize[0,1],
  eyx1=O.cell.hSize[1,0],
  Z=avgNumInteractions(),
  Zm=avgNumInteractions(skipFree=True),
                voidratio=(O.cell.size[0]*O.cell.size[1]*O.cell.size[2]-a)/a,
# poros=voxelPorosity(500,(0,0,0),O.cell.size),
  unbalanced=utils.unbalancedForce(),
  t=O.time,
  inwork=work
  gWork=O.energy['gravWork'],
# Ep=O.energy['elastPotential'],
                Ep=law.elasticEnergy()-Ep_ini,
  Edamp=O.energy['nonviscDamp']-Edamp_ini,
                Ediss=law.plasticDissipation(),
# Ediss=O.energy['plastDissip'],
  Ekin=utils.kineticEnergy(),
  Etot=law.elasticEnergy()-Ep_ini+O.energy['nonviscDamp']-Edamp_ini+law.plasticDissipation()+utils.kineticEnergy()-Ekin_ini,
# Etot=O.energy.total(),**O.energy

 )
 plot.saveDataTxt('macroFile',vars=('t','exx','eyy','ezz','sxx','syy','szz','Z','Zm','voidratio'))
 plot.saveDataTxt('energyFile',vars=('t','unbalanced','gWork','Edamp','Ekin','Ep','Ediss','Etot','inwork'))

O.cell.trsf=Matrix3.Identity
O.cell.velGrad=Matrix3(-0.1,0,0,0,0.05,0,0,0,0.05)
O.run()

Hi,

Sorry, but I cannot run this. There are some syntax errors and the script requires loading external files. Please try to prepare MWE [2]:

"Do your best to send it as a MWE (Minimal Working Example™). The MWE should actually be the smallest script possible showing the same problem you faced probably in a much longer script, initially. It must be:

- minimal = short. It is difficult to find mistakes in long pieces of code and it is usually much more annoying. Any post-processing command should for instance be banned.
- working = before sending, test it to be sure it is working (it has no syntax errors etc.)."

Cheers,
Karol

[2] https://www.yade-dem.org/wiki/Howtoask

Hi Karol,

Sorry, the format of the code changes when i copy and paste it. I have already corrected it shown as below. There are two scripts, first is the consolidation and the second is the undrained test. After you run through the first script, a file of consolidated sample will be get, which can be directly loaded by the second script.

> - minimal = short. It is difficult to find mistakes in long pieces of code and it is usually much more annoying.
Sorry for that, I have done my best to shorten the code.

This is the corrected script
######Consolidation#######
from yade import pack,plot,qt,export,ymport,utils
import numpy as np
IsoSigma = -1.e5
O.periodic=True
O.trackEnergy=True
idSand=O.materials.append(FrictMat(young=1e8,poisson=.8,frictionAngle=0.2,density=26500000,label='spheres'))
sp1=pack.SpherePack()
sp1.makeCloud(maxCorner=(0.1, 0.1, 0.1), psdSizes=[0.0017, 0.00191, 0.002285, 0.0026, 0.00292, 0.00325, 0.0035], psdCumm=[0.0, 0.1, 0.3, 0.5, 0.6, 0.9, 1], periodic=True,num=2000,seed=1)
sp1.toSimulation(color=(0,0,1),material=idSand)
v=utils.getSpheresVolume()
print utils.getSpheresVolume()
fout = file('Packing_volume.txt','w')
fout.write(str(v))
fout.close()

relRadList2 = [1,1]
relPosList2 = [[0.4,0,0],[-0.4,0,0]]

templates= []
templates.append(clumpTemplate(relRadii=relRadList2,relPositions=relPosList2))
O.bodies.replaceByClumps(templates,[1.],discretization=10)

law=Law2_ScGeom_FrictPhys_CundallStrack(traceEnergy=True)
O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb()]),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom()],
 [Ip2_FrictMat_FrictMat_FrictPhys()],
 [law]),
 GlobalStiffnessTimeStepper(),
 NewtonIntegrator(damping=0.20),
 PeriTriaxController( goal=(IsoSigma,IsoSigma,IsoSigma), # Vector6 of prescribed final values
    stressMask=7,
    dynCell=True,
    maxStrainRate=(5.e+0,5.e+0,5.e+0),
    maxUnbalanced=0.001,
    relStressTol=1.e-3,
    doneHook='Finished()',
    label='p3d'
 ),
]

def Finished():
 O.save('post_iso.xml')
 print 'Finished'
 print (O.cell.size[0]*O.cell.size[1]*O.cell.size[2]-v)/v
 print getStress()
 print unbalancedForce()
 O.pause()

O.run()

########Undrained triaxial test#########
from yade import pack,plot,qt,export
from math import fabs
import numpy as np
O.load('post_iso.xml')
O.trackEnergy=True
setContactFriction(0.5)
filename1='Packing_volume.txt'
a=np.loadtxt(filename1)
Ep_ini=O.energy['elastPotential']
Edamp_ini=O.energy['nonviscDamp']
Ekin_ini=utils.kineticEnergy()
work=0

law=Law2_ScGeom_FrictPhys_CundallStrack(traceEnergy=True)
O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb()]),
 InteractionLoop(
 [Ig2_Sphere_Sphere_ScGeom()],
 [Ip2_FrictMat_FrictMat_FrictPhys()],
 [law]),
 NewtonIntegrator(damping=0.2),
 GlobalStiffnessTimeStepper(),
 PyRunner(command='inputwork()',iterPeriod=1),
 PyRunner(command='plotAddData()',iterPeriod=10),
]

def inputwork():
 global work
 dwork=(0.5*(O.cell.prevVelGrad+O.cell.velGrad)*utils.getStress()).trace()*O.dt*(O.cell.hSize).determinant()
 work=work+dwork

def plotAddData():
 global Ep_ini
 global Edamp_ini
 global Ekin_ini
 global work
 plot.addData(
  iter=O.iter,iter_=O.iter,
  sxx=utils.getStress()[0,0],
  syy=utils.getStress()[1,1],
  szz=utils.getStress()[2,2],
  exx=O.cell.getSmallStrain()[0,0],
  eyy=O.cell.getSmallStrain()[1,1],
  ezz=O.cell.getSmallStrain()[2,2],
  Z=avgNumInteractions(),
  Zm=avgNumInteractions(skipFree=True),
  voidratio=(O.cell.size[0]*O.cell.size[1]*O.cell.size[2]-a)/a,
  unbalanced=utils.unbalancedForce(),
  t=O.time,
  inwork=work,
  gWork=O.energy['gravWork'],
# Ep=O.energy['elastPotential'],
  Ep=law.elasticEnergy()-Ep_ini,
  Edamp=O.energy['nonviscDamp']-Edamp_ini,
  Ediss=law.plasticDissipation(),
# Ediss=O.energy['plastDissip'],
  Ekin=utils.kineticEnergy(),
  Etot=law.elasticEnergy()-Ep_ini+O.energy['nonviscDamp']-Edamp_ini+law.plasticDissipation()+utils.kineticEnergy()-Ekin_ini,
# Etot=O.energy.total(),**O.energy

 )
 plot.saveDataTxt('macroFile',vars=('t','exx','eyy','ezz','sxx','syy','szz','Z','Zm','voidratio'))
 plot.saveDataTxt('energyFile',vars=('t','unbalanced','gWork','Edamp','Ekin','Ep','Ediss','Etot','inwork'))

O.cell.trsf=Matrix3.Identity
O.cell.velGrad=Matrix3(-0.1,0,0,0,0.05,0,0,0,0.05)
O.run()

###########################################
My problom is the difference of "inwork" and "Etot" in the second script when using Clump particles.
I really appreciate if you can confirm the issue or help me to solve the problem.
Thanks a lot,

Xu

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.

> PS Jan, sorry for suggesting this 'interaction inside clump' issue, but I can swear I encountered this some time ago.

at the time you proposed this reason, i.e. without a any code, it was relevant option IMO (I am not sure what Yade does with existing interactions if you clump already interacting particles and it is possible it just leaves it untouched?)

Cheers
Jan

Hi Karol and Jan,

Thank you for your MWE. It looks like 'utils.getStress()' may give slightly different results for clumps and spheres. However, I have used the clump which is made up of more than 2 spheres (3 or 4 or dozens). The more the number of spheres members in clump, the larger the difference between 'inwork' and 'Etot', when the volum of clump remains the same. And i think the slightly different results of 'utils.getStress()' may not cause such obvious error between input work and microscale energy. Because the value of total microscale energy is more than twice the input work value when the number of clump members reaches 40-50.

So that's why I said the cause of the difference may lie in line 20,21 in [1]. As the particle volume remains the same, the radius of sphere member in clump becomes smaller when the number of memeber increases. And the kn and ks calculated in [2] is smaller, hence, the microscale energy computed by [3-4] is larger than the real value. So i guess this is the reason why 'Etot' is larger than 'inwork' .

I think my opion may be wrong, but I can't find any other reason. So if you can give me some more reasons for this phenomenan, i would really appreciate it.

Best wishes,
Xu

[1] https://gitlab.com/yade-dev/trunk/-/blob/master/pkg/dem/FrictPhys.cpp#L20-L21
[2] https://gitlab.com/yade-dev/trunk/-/blob/master/pkg/dem/FrictPhys.cpp#L35-L37
[3] https://gitlab.com/yade-dev/trunk/-/blob/master/pkg/dem/ElasticContactLaw.cpp#L91-L94
[4] https://gitlab.com/yade-dev/trunk/-/blob/master/pkg/dem/ElasticContactLaw.cpp#L97-L100

Hi, that changing members size results in different stiffness is a correct observation.
It is not "wrong", in the sense that it's a default behavior and you can always counter this effect by adapting Young parameter of the material.
Anyway, that stiffness is different is not a reason for not conserving energy.
I noticed that you don't take kinetic energy into account, is it negligible?
Bruno

Hi Bruno,
Thank you for your reply.

> I noticed that you don't take kinetic energy into account, is it negligible?
Kinetic energy has already been considered in MWE(#9) above. Anyway, kinetic energy is negligible for not conserving energy.

The microscale energy calculated by:
Law=Law2_ScGeom_FrictPhys_CundallStrack(traceEnergy=True)
Estrain=Law.elasticEnergy()
Efriction=Law.plasticDissipation()
Edamp=O.energy['nonviscDamp']
Ekin=utils.kineticEnergy()
Etot=law.elasticEnergy()-Ep_ini+O.energy['nonviscDamp']-Edamp_ini+law.plasticDissipation()+utils.kineticEnergy()-Ekin_ini,

The input work done by external stress:
dw=(0.5*(O.cell.preVelGrad+O.cell.velGrad)*utils.getStress())*O.dt*(O.cell.hSize).determinant()
W=W+dw

The two parts being the same when using sphere particles . However, the microscale energy ('Etot') is greater than input work when using clump particles. Other details can be seen in MWE#9.

Cheers,
Xu

I found the problem, getStress does not handle clumps...

In the product f*branch, "branch" should point to one single point per body (in this case: per clump). Instead, the function implementation disregards the clump status and uses the center of the clump members.

I reported the issue here: https://gitlab.com/yade-dev/trunk/-/issues/289

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()

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

Hi Bruno,

I modified the "branch" for clump, and it works and solved my problem.
I would follow up this issue later.

Thank you,
Xu