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

Hello,

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

Cheers,

Robert

[1]https://yade-dem.org/doc/yade.wrapper.html#yade.wrapper.Cell.velGrad

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.periodic=True
O.materials.append(FrictMat(density=2000,young=100e6,poisson=0.2,frictionAngle=radians(30),label='spheres'))
O.cell.hSize=Matrix3(0.03,0,0, 0,0.03,0, 0,0,0.03)
sp=pack.SpherePack()
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])

O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom(), Ig2_Box_Sphere_ScGeom()],
  [Ip2_FrictMat_FrictMat_MindlinPhys(betan=0.3,betas=0.3)],
  [Law2_ScGeom_MindlinPhys_Mindlin()] ),
 NewtonIntegrator(gravity=(0,0,-9.81),damping=0.0),
 PyRunner(command='check()',iterPeriod=1000),
]
O.dt=0.3*utils.PWaveTimeStep()
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():
 kineE=0
 for i in O.bodies:
  m=i.state.mass
  v2=numpy.dot(i.state.vel, i.state.vel)
  I=sqrt(numpy.dot(i.state.inertia,i.state.inertia))
  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,
Yuxuan

Hello,

both values are computed differently.
E=0.5*m*v^2+0.5*I*w^2
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).

Cheers
Jan

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

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,
Yuxuan

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,
Yuxuan

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]

Cheers
Jan

[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

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!
Yuxuan