energy issues

Asked by ytang on 2020-11-06

Hi all,

recently, I'm trying to study energy consumption in YADE, so I just use the gravity deposition process to check the energy dissipation process.

here is the code: the code is almost the same as the example in YADE [1].
#########################################################################
from yade import pack, plot
damping = 0.4
O.bodies.append(geom.facetBox((.5,.5,.5),(.5,.5,.5),wallMask=31))
sp=pack.SpherePack()
sp.makeCloud((0,0,0),(1,1,1),rMean=.23,rRelFuzz=0.0,seed =1)
sp.toSimulation()
print(len(O.bodies))

for b in O.bodies:
 if isinstance(b.shape,Sphere):
  positionz = b.state.pos[2]
  mass = b.state.mass
  print(positionz)
  print(mass)
  print("#######")
print("############################################")
for b in O.bodies:
 if isinstance(b.shape,Sphere):
  radius = b.shape.radius
  print(radius)
O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb()]),
 InteractionLoop(
  # handle sphere+sphere and facet+sphere collisions
  [Ig2_Sphere_Sphere_ScGeom(),Ig2_Facet_Sphere_ScGeom()],
  [Ip2_FrictMat_FrictMat_FrictPhys()],
  [Law2_ScGeom_FrictPhys_CundallStrack(traceEnergy = True, label = 'law')]
  #[Ip2_FrictMat_FrictMat_MindlinPhys( betan =0.2 , betas = 0.2, krot = 0.2 , eta = 0.5 )],
  #[Law2_ScGeom_MindlinPhys_Mindlin()]
 ),
 #NewtonIntegrator(gravity=(0,0,-9.81),kinSplit=True,damping=damping),
 NewtonIntegrator(gravity=(0,0,-9.81),damping=0.4),
 PyRunner(command='checkUnbalanced()',realPeriod=2),
 PyRunner(command='addPlotData()',iterPeriod=2000),
 #PyRunner(command='record_all_info()',iterPeriod=5000)
]
O.dt=.5*PWaveTimeStep()
O.trackEnergy=True
def checkUnbalanced():
 if unbalancedForce()<.001:
  O.pause()
  plot.saveDataTxt('bbb.txt.bz2')
######################################################################################
######################################################################################
###################### method 1, energy 1 ##############################

"""
def addPlotData():
 # each item is given a names, by which it can be the unsed in plot.plots
 # the **O.energy converts dictionary-like O.energy to plot.addData arguments
 elastic_part = law.elasticEnergy()
 plastic_part = law.plasticDissipation()
 E_kin_translation = 0
 E_kin_rotation = 0
 E_pot = 0
 E_plastic = 0
 E_nonviscoDamp = 0
 E_tracker = dict(list(O.energy.items()))

 E_kin_translation = E_tracker['kinTrans']
 E_kin_rotation = E_tracker['kinRot']
 E_pot = E_tracker['gravWork']
 if('plastDissip' in E_tracker):
  E_plastic = E_tracker['plastDissip']
 if(damping!=0):
  E_nonviscoDamp = E_tracker['nonviscDamp']
 plot.addData( elasticEnergy = elastic_part, plasticenergy = plastic_part,
  translationenergy = E_kin_translation, rotationenergy = E_kin_rotation,
  gravitywork = E_pot, plasticdissip = E_plastic, nonviscodamp = E_nonviscoDamp,
  i=O.iter,unbalanced=unbalancedForce())
     plot.saveDataTxt('energy.txt')
"""
#####################################################################################
################### method 2, energy 2 #####################################
def addPlotData():
 plot.addData(i=O.iter,unbalanced=unbalancedForce(),**O.energy)
plot.plots={'i':('unbalanced',None,O.energy.keys)}
plot.plot()
O.saveTmp()
################################################################################
################################################################################
I use two methods to get the data, as you can see the annotation from the code. I put the data in this folder [2].
###################
I have some basic questions.
(1) how do we calculate the gravity work?
In my simulation, there are four particles. As we all know, the particles should have the maximum gravity work at the height position (which means at the start position).
I get the z position of the particles and use this equation (m*g*h) to calculate the gravity work. The result is 1022.689 for these four particles. However, the recorded data shows that gravity work changes from -461.95 to -510.83. (as for -461.95, this data is recorded at timesteps = 2000. ) I find the source code on how to calculate the gravity work [3], but I do not fully understand it.

(2) I plot the curves on the energy changes with the timesteps. it seems the energy is not conserved. I do not know whether the energy in YADE is conserved. Can anyone tell me about this?

The reason why I say the energy is not conserved: if we add all the energy at every certain timestep, the result is not zero. for instance :
i elastPotential gravWork kinRot kinTrans nonviscDamp
72000 0.344990744 -510.833081 8.65006E-28 2.57632E-27 291.8306236
Here the total energy is not zero.

(3) I compare the results from method one and method two at the final stage:
##################################################################
energy1 i elasticEnergy gravitywork rotationenergy translationenergy nonviscodamp plasticdissip plasticenergy
 62000 0.344990744 -510.833081 1.44064E-27 4.73862E-27 291.8306236 0 145.4170567

energy2 i elastPotential gravWork kinRot kinTrans nonviscDamp
 72000 0.344990744 -510.833081 8.65006E-28 2.57632E-27 291.8306236
###############################################################################
The results for these two methods are almost the same for some of the energies. However, there is no plastic dissipation in method 2 (although plastic dissipation in method one is zero).
In my opinion, plastDissip can be automatically recorded when we use this command (O.trackEnergy=True) even if the value is zero.
why the plastDissip does not being recorded here?
what is the relationship between the plastDissip and the plasticDissipation[4]? (it seems they are not the same from method one. one is zero, the other one is 145.41. )

I'm sorry that I have several questions.
thanks!

references: [1] https://yade-dem.org/doc/tutorial-examples.html
[2]https://www.dropbox.com/sh/fd13w0yednmiz8g/AABkLkUYBiInTSVZpKZmcSYya?dl=0
[3] https://gricad-gitlab.univ-grenoble-alpes.fr/cailletr/yade/-/blob/a583fdbea2c5b63f5d26b96655a6e05e8bccf375/pkg/dem/NewtonIntegrator.cpp
[4] https://yade-dem.org/doc/yade.wrapper.html?highlight=law2_scgeom_frictphys_cundallstrack#yade.wrapper.Law2_ScGeom_FrictPhys_CundallStrack

Question information

Language:
English Edit question
Status:
Answered
For:
Yade Edit question
Assignee:
No assignee Edit question
Last query:
2020-11-06
Last reply:
2020-11-11

This question was reopened

Jan Stránský (honzik) said : #1

Hello,

> (1) how do we calculate the gravity work?
> ... but I do not fully understand it.

see [1]
delta(GravWork) = - gravity.dot(vel) * mass * dt # per time step

this is equivalent to delta(m*g*h) = m*g*delta(h), delta(h) being vel*dt approximation

> As we all know, the particles should have the maximum gravity work at the height position ... I get the z position of the particles and use this equation (m*g*h) to calculate the gravity work...

Yade computes actual work done by gravity in the simulation.
The *potential energy* mgh uses some reference frame (the h term). This frame is arbitrary (consider falling on a flat roof, different h would be w.r.t the roof and w.r.t the ground) and hence the initial value is arbitrary.
If needed, consider the Yade's reference frame such that the initial potential energy is 0.
If you know h in advance, you can add it "manually".

> (2) ...it seems the energy is not conserved. I do not know whether the energy in YADE is conserved.

Potential energy m*g*h is not included.
You have some initial potential energy and after some run, it is lowered by gravity work, but should be included to energy balance if you are interested in the balance itself.

I suggest to split the problem and investigate simpler scenarios:
- free fall (without any interactions)
- fall on a ground (with prescribed initial velocity, no gravity)
- and after understanding both previous scenarios, go to the combined case

> (3)
> why the plastDissip does not being recorded here?

plastDissip is only recorded when "needed", when the pastic dissipation actually happen, not before (it is inside "if (shearForce.squaredNorm() > maxFs)" branch).
You have it in method 1, because you explicitly add E_plastic (default value is 0) to plod.addData.

> what is the relationship between the plastDissip and the plasticDissipation[4]?

it is the same value. See [3] for the difference.
For Law2(traceEnergy=True), it is saved only to law.plasticDissipation() and NOT to O.energies. That's why you get 0 for plasticdissip.
It works oppositely for Law2(traceEnergy=False).

cheers
Jan

[1] https://gitlab.com/yade-dev/trunk/-/blob/master/pkg/dem/NewtonIntegrator.cpp#L90
[2] https://gitlab.com/yade-dev/trunk/-/blob/master/pkg/dem/ElasticContactLaw.cpp#L101
[3] https://gitlab.com/yade-dev/trunk/-/blob/master/pkg/dem/ElasticContactLaw.cpp#L98

Jérôme Duriez (jduriez) said : #2

Regarding your initial ref [3] https://gricad-gitlab.univ-grenoble-alpes.fr/cailletr, I would rather advice to online browse YADE source code at https://gitlab.com/yade-dev/trunk

I'm unsure how closely bonded these two repositories are, at the moment.

Can you help with this problem?

Provide an answer of your own, or ask ytang for more information if necessary.

To post a message you must log in.