← Back to team overview

yade-users team mailing list archive

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 yade-users is
an answer contact for Yade.