I calculated the kinetic energy by E=0.5mv^2+0.5Iw^2 but it is not equal to kineticEnergy() when O.cell.velGrad is set

Asked by Yuxuan Wen

Dear YADE developers and users,

I am conducting a simulation where the energy results are important for the study. So I want to study how the kinetic energy is calculated in YADE. In the documentation it shows kinetic energy is calculated by E=0.5*m*v^2+0.5*I*w^2, which is definitely normal and accurate. So I have tried to calculated kinetic energy by myself, as the particle mass m, velocity v, inertia I and angular velocity w can be all obtain by using O.bodies[]state.mass, O.bodies[].state.vel, O.bodies[].state.inertia and O.bodies[].state.angVel. But when I have set the value of O.cell.velGrad, the E=0.5*m*v^2+0.5*I*w^2 calculated by myself is much larger than utils.kineticEnergy(). And when I didn't use O.cell.velGrad, the E=0.5*m*v^2+0.5*I*w^2 calculated by myself is the same with utils.kineticEnergy()

I am thinking why this happened? Why setting the value of O.cell.velGrad will make E=0.5*m*v^2+0.5*I*w^2 larger than the utils.kineticEnergy()?

My code for calculating the E=0.5*m*v^2+0.5*I*w^2 is:

for i in bodies_id:
        v2=numpy.dot(O.bodies[i].state.vel, O.bodies[i].state.vel)
 w2=numpy.dot(O.bodies[i].state.angVel, O.bodies[i].state.angVel)

Thank you very much!!

Question information

English Edit question
Yade Edit question
No assignee Edit question
Solved by:
Yuxuan Wen
Last query:
Last reply:
Revision history for this message
Robert Caulk (rcaulk) said :


Please confirm you have read [1], and are not inaccurately setting velGrad.




Revision history for this message
Yuxuan Wen (wenyuxuan) said :

Thank you Robert for your reply!

Yes, I set the velGrad in the format as shown in the documentation: O.cell.velGrad=Matrix3(-0.5,0,0,0,-0.5,0,0,0,-0.5). And then the kinetic energy calculated by E=0.5*m*v^2+0.5*I*w^2 is different with the kineticEnergy(). I made a simplified code of the simulation, as shown follows. When you run this code directly, the kinetic energy is different, and when you comment the line of O.cell.velGrad, the kinetic energy will be the same.

from yade import pack, plot, qt, export, os
O.cell.hSize=Matrix3(0.03,0,0, 0,0.03,0, 0,0,0.03)
sp.makeCloud(Vector3(0,0,0), Vector3(0.03,0.03,0.03), num=1000, rMean=0.001, rRelFuzz=0, distributeMass=False, periodic=True, seed=1)
O.bodies.append([sphere(center,rad,material='spheres') for center,rad in sp])

  [Ig2_Sphere_Sphere_ScGeom(), Ig2_Box_Sphere_ScGeom()],
  [Law2_ScGeom_MindlinPhys_Mindlin()] ),
O.cell.velGrad=Matrix3(-0.5,0,0, 0,-0.5,0, 0,0,-0.5) # if comment this line, kineE will be the same with kineticEnergy()

def check():
 for i in O.bodies:
  v2=numpy.dot(i.state.vel, i.state.vel)
  w2=numpy.dot(i.state.angVel, i.state.angVel)
  kineE=kineE+0.5*m*v2+0.5*I*w2 # compare with utils.kineticEnergy()
 print(kineE, kineticEnergy())

Thank you so much!

Kind Regards,

Revision history for this message
Jan Stránský (honzik) said (last edit ):


both values are computed differently.
of course still holds, but the values of "v" and "w" are different.
Have a look at source code [2]: conditions "if (scene->isPeriodic)", variables like "spin = scene->cell->getSpin()" and their usage, comments "Only take in account the fluctuation velocity" etc.
If still anything is not clear, please let us know.

If you think this discrepancy should be mentioned in the documentation, feel free to update it yourself or create an issue on gitlab pages of the project (not to be forgotten).


[2] https://gitlab.com/yade-dev/trunk/-/blob/master/pkg/dem/Shop_01.cpp#L193

Revision history for this message
Yuxuan Wen (wenyuxuan) said :

Hello Jan,

Thank you for you help and sending me the position of the source code! I am much appreciated!

But I am still a little confused with the definition of the fluctuation velocity. As you said, the E=0.5*m*v^2+0.5*I*w^2 still holds, and "v" and "w" are the fluctuation velocity . In the source code, the fluctuation velocity is defined as:
bodyFluctuationVel(state->pos - state->vel * scene->dt, state->vel, scene->cell->velGrad)
which is a function of 1)updated position, 2) v, and 3) velGrad.
But I don't know how to calculate the fluctuation velocity by these 3 terms. I searched the documentation by "bodyFluctuationVel" and didn't find anything matched. I wonder if it is possible that you could tell me the expression of the fluctuation velocity? Thank you!

Kind Regards,

Revision history for this message
Yuxuan Wen (wenyuxuan) said :

Hello Jan,

Thanks for your suggestion, now I have created an issue on gitlab pages of the project [3]. Because I am still confused with the definition of fluctuation velocity, so I didn't update it myself but created an issue on the gitlab pages.

[3] https://gitlab.com/yade-dev/trunk/-/issues/236

Kind Regards,

Revision history for this message
Jan Stránský (honzik) said :

Probably again only source code is relevant here [4].

The idea is this:
Let's have a particle at (0,0,0), i.e. one corner of periodic cell.
Properties of the particle should (must) be identical with its periodic images (in this case e.g. particles in other corners of the periodic cell). See [5] and pictures there for a reference.
Let's say the (0,0,0) particle has zero velocity.
Now consider the cell has some velGrad. E.g. it is stretching in x direction. Then he "image particles" in the "right face" of the cell has some non-zero positive velocity.
Now, from the periodicity point of view, it should not matter if you choose the (0,0,0) particle or any of its periodic images.
The (0,0,0) particle has zero velocity, the periodic images have non-zero velocity, which one should be chosen?? The "fluctuation part" of the velocity is chosen in Yade.

Similarly it is with angular velocity and cell spin [6]


[4] https://gitlab.com/yade-dev/trunk/-/blob/master/core/Cell.hpp#L121
[5] https://yade-dem.org/doc/formulation.html#periodic-boundary-conditions
[6] https://gitlab.com/yade-dev/trunk/-/blob/master/core/Cell.hpp#L216

Revision history for this message
Yuxuan Wen (wenyuxuan) said :

Thank you so much Jan!! Now I understand the definition of fluctuation velocity and why it is chosen in Yade!
Thanks again for you help!! I am much appreciated!