TorqueEngine, angular acceleration and moment of inertia of Boxes

Hi Behzad,

The inertia of a box body is computed in [1], assuming you're using the box() function of utils module. How this inertia is used to compute angular acceleration according to torque is in NewtonIntegrator.cpp

I am not familiar with this part of the code for non-spherical bodies, but it is possible you get "wrong" results with default parameters in your simulation. The default computations seem to assume spherical bodies. To get consistent results with aspherical bodies like your box, I guess you have to use aspherical flags. Defined a the body scale [2], and at NewtonIntegrator engine scale [3]

A side note : do you need dynamic (whose movements depend on sum of forces / torques, and are not user-controlled) boxes ?

Jerome

[1] https://github.com/yade/trunk/blob/dfbb8dde95d305ed95add159a77d8bcb0b3de935/pkg/dem/Shop_01.cpp#L260
[2] https://yade-dem.org/doc/yade.wrapper.html#yade.wrapper.Body.aspherical
[3] https://yade-dem.org/doc/yade.wrapper.html#yade.wrapper.NewtonIntegrator.exactAsphericalRot

Hi Jerome,

Thanks for the answer.

A side note : do you need dynamic (whose movements depend on sum of forces / torques, and are not user-controlled) boxes ?

- I sandwich a sphere pack between two plates (boxes). Then, I apply a torque on the upper plate and I measure it's rotational velocity. So, I think it must be dynamic.

I applied a torque on the box without any spheres in contact and the acceleration/velocity was wrong.

I give a try to aspherical flag to see if it can solve the problem.

cheers,
Behzad

It doesn't work even with exactAsphericalRot=True!

The equations in [1] are fine.

we've got:

body->state->inertia=Vector3r(mass*(4*extents[1]*extents[1]+4*extents[2]*extents[2])/12.,mass*(4*extents[0]*extents[0]+4*extents[2]*extents[2])/12.,mass*(4*extents[0]*extents[0]+4*extents[1]*extents[1])/12.);

It's the correct theoretical equation. However, I don't know why the answer of this equation is not correct even when "exactAsphericalRot=True"

[1] https://github.com/yade/trunk/blob/dfbb8dde95d305ed95add159a77d8bcb0b3de935/pkg/dem/Shop_01.cpp#L260

here's the script:

O.reset()
from yade import utils, plot
from yade import pack, qt

#=============================================Engines===============================================================

O.engines=[
ForceResetter(),
TorqueEngine(ids=[0],moment=(0,0,1e5)),
GlobalStiffnessTimeStepper(),
InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb(),Bo1_Facet_Aabb()]),
InteractionLoop(
[Ig2_Sphere_Sphere_ScGeom6D(),Ig2_Box_Sphere_ScGeom6D()],
[Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(setCohesionNow=True, setCohesionOnNewContacts=True)],
[Law2_ScGeom6D_CohFrictPhys_CohesionMoment(),Law2_ScGeom_FrictPhys_CundallStrack()]
),
NewtonIntegrator(damping=0.7,gravity=[0,0,0],exactAsphericalRot=True)
]

id_Mat1=O.materials.append(CohFrictMat(normalCohesion= 1e15, shearCohesion= 1e13, isCohesive= True, young=1e5,
density=2000, poisson=0.3, frictionAngle= 0.4))
CohMat=O.materials[id_Mat1]

upper_plate = O.bodies.append(box((0,0,0),(1,1,0.1),fixed=True,material=CohMat))

O.bodies.append(sphere((0,0,1),0.05,fixed=True, material=CohMat))

#=========================================Boundary Conditions=======================================================

O.bodies[upper_plate].aspherical=True

O.engines=O.engines+[PyRunner(iterPeriod=1,command='addPlotData()',label='plotDataCollector')]

plot.plots={'t':('eps',),' t':('rot'),}

def addPlotData():
       yade.plot.addData(t=O.time,pos=O.bodies[upper_plate].state.pos[2],eps=(O.bodies[upper_plate].state.angVel[2])/O.time,
   rot=O.bodies[upper_plate].state.rot()[2])

plot.plot()

>"NewtonIntegrator(damping=0.7,..."

Erm... do you know what newton damping is?
Did you put 0.7 because it is the default in PFC??

Yes, this value came from PFC!

I changed it to zero. Now the ratio of the applied torque to the obtained angular acceleration is correct.
However, the problem is the calculation of inertia is not correct.

we have:

upper_plate = O.bodies.append(box((0,0,0),(1,1,0.1),fixed=True,material=CohMat))

]: O.bodies[upper_plate].state.mass
 -> [2]: 1600.0

So from

body->state->inertia=Vector3r(mass*(4*extents[1]*extents[1]+4*extents[2]*extents[2])/12.,mass*(4*extents[0]*extents[0]+4*extents[2]*extents[2])/12.,mass*(4*extents[0]*extents[0]+4*extents[1]*extents[1])/12.);

we must have inertia as (538.66,538.66,1066.66)

But we get:

 O.bodies[upper_plate].state.inertia
 -> [3]: Vector3(8080,8080,16000)

Where does it come from?

Just a sidenote, if you use box, the mass and inertia is computed in
py/utils.py file [1,2,3], not in c++ as you refer to.. in utils.py, the
inertia is computed as you cited, but with density instead of mass.. so
evidently a bug :-) feel free to report it
Jan

PS: thanks for notifying

b1 = b=box((0,0,0),(1,1,1))
b2 = b=box((0,0,0),(2,2,2))
print b1.state.inertia, b2.state.inertia

Vector3(8000,8000,8000) Vector3(32000,32000,32000) # but should be 32x more
(8x more mass x 4x more geometry inertia), not 4x more...

[1]
http://bazaar.launchpad.net/~yade-pkg/yade/git-trunk/view/head:/py/utils.py#L214
[2]
http://bazaar.launchpad.net/~yade-pkg/yade/git-trunk/view/head:/py/utils.py#L116
[3]
http://bazaar.launchpad.net/~yade-pkg/yade/git-trunk/view/head:/py/utils.py#L128

2015-03-02 19:36 GMT+01:00 behzad <email address hidden>:

> Question #262922 on Yade changed:
> https://answers.launchpad.net/yade/+question/262922
>
> behzad posted a new comment:
> Yes, this value came from PFC!
>
> I changed it to zero. Now the ratio of the applied torque to the obtained
> angular acceleration is correct.
> However, the problem is the calculation of inertia is not correct.
>
> we have:
>
> upper_plate =
> O.bodies.append(box((0,0,0),(1,1,0.1),fixed=True,material=CohMat))
>
> ]: O.bodies[upper_plate].state.mass
> -> [2]: 1600.0
>
>
> So from
>
>
> body->state->inertia=Vector3r(mass*(4*extents[1]*extents[1]+4*extents[2]*extents[2])/12.,mass*(4*extents[0]*extents[0]+4*extents[2]*extents[2])/12.,mass*(4*extents[0]*extents[0]+4*extents[1]*extents[1])/12.);
>
> we must have inertia as (538.66,538.66,1066.66)
>
> But we get:
>
> O.bodies[upper_plate].state.inertia
> -> [3]: Vector3(8080,8080,16000)
>
> Where does it come from?
>
> --
> You received this question notification because you are a member of
> yade-users, which is an answer contact for Yade.
>
> _______________________________________________
> Mailing list: https://launchpad.net/~yade-users
> Post to : <email address hidden>
> Unsubscribe : https://launchpad.net/~yade-users
> More help : https://help.launchpad.net/ListHelp
>