# 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 [1].

#######

from yade import pack, plot

damping = 0.4

O.bodies.

sp=pack.

sp.makeCloud(

sp.toSimulation()

print(len(

for b in O.bodies:

if isinstance(

positionz = b.state.pos[2]

mass = b.state.mass

print(positionz)

print(mass)

print("#######")

print("

for b in O.bodies:

if isinstance(

radius = b.shape.radius

print(radius)

O.engines=[

ForceResetter(),

InsertionSortC

InteractionLoop(

# handle sphere+sphere and facet+sphere collisions

[Ig2_

[Ip2_

[Law2_

#[Ip2_

#[Law2_

),

#NewtonIntegra

NewtonIntegrat

PyRunner(

PyRunner(

#PyRunner(

]

O.dt=.5*

O.trackEnergy=True

def checkUnbalanced():

if unbalancedForce

O.pause()

plot.

#######

#######

#######

"""

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

E_kin_translation = 0

E_kin_rotation = 0

E_pot = 0

E_plastic = 0

E_nonviscoDamp = 0

E_tracker = dict(list(

E_kin_translation = E_tracker[

E_kin_rotation = E_tracker['kinRot']

E_pot = E_tracker[

if('plastDissip' in E_tracker):

E_plastic = E_tracker[

if(damping!=0):

E_nonviscoDamp = E_tracker[

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,

plot.

"""

#######

################### method 2, energy 2 #######

def addPlotData():

plot.addData(

plot.plots=

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=

why the plastDissip does not being recorded here?

what is the relationship between the plastDissip and the plasticDissipat

I'm sorry that I have several questions.

thanks!

references: [1] https:/

[2]https:/

[3] https:/

[4] https:/

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

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 plasticDissipat

it is the same value. See [3] for the difference.

For Law2(traceEnerg

It works oppositely for Law2(traceEnerg

cheers

Jan

[1] https:/

[2] https:/

[3] https:/

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

Regarding your initial ref [3] https:/

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.