yadeusers team mailing list archive

yadeusers team

Mailing list archive

Message #26077
Re: [Question #698392]: 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
Question #698392 on Yade changed:
https://answers.launchpad.net/yade/+question/698392
Status: Needs information => Open
Yuxuan Wen gave more information on the question:
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

You received this question notification because your team yadeusers is
an answer contact for Yade.