# energy issues

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 .
#########################################################################
damping = 0.4
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
mass = b.state.mass
print(positionz)
print(mass)
print("#######")
print("############################################")
for b in O.bodies:
if isinstance(b.shape,Sphere):
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='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 ##############################

"""
# 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 #####################################
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 .
###################
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 , 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? (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!

## Question information

Language:
English Edit question
Status:
For:
Assignee:
No assignee Edit question
Last query:
 Revision history for this message Jan Stránský (honzik) said on 2020-11-11: #1

Hello,

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

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

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

it is the same value. See  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

 Revision history for this message Jérôme Duriez (jduriez) said on 2020-11-12: #2