yade-users team mailing list archive

["Sergei D." <sj2001@xxxxxxxxx>]

> Perhaps I don't understand what you need, but I am a bit surprised to
> hear it needs prevAngMom.
> In the current scheme,  at time t, what you find in body->velocity is in
> fact velocity at time t-dt/2.
> 
> Introducing prevAngMom means you want a dependancy between
> AngMom(t+dt/2) and AngMom(t-3*dt/2).
> It can happen for high (4th) order time discretisation, but it is really
> what you want now?
> 


I want to use modified leap-frog algorithm (see Allen, M. P. &
Tildesley, D. J. Computer simulation of liquids Clarendon Press, 1987, p.89)

Let l is the angular moment,

l=Iw,

where I is the inertia tenzor and w is the angular velocity; all in the
global reference frame (r.f.). Let also (.)' will be a value of (.) in
the local (or body's) r.f.

We have:

l(t) = l(t-dt/2) + dt/2 T(t),

where dt is the time step; T is the torque (in global r.f.).
When,

l'(t) = Q(t)*l(t)*Q_(t) = I'w'(t)

where Q={q0,q1,q2,q3} is the quaternion defining the orientation of a
body; I'={I1,I2,I3} is the principal inertia.

So, from above we have local r.f. angular velocity for time t:

w'(t)_i = l'(t)_i/I'_i, i=1,2,3.

Quaternion Q and angular velocity w' of a body are related by eq.(*):

d(q0)/dt = 1/2 ( - q1 w'1 - q2 w'2 - q3 w'3 )
d(q1)/dt = 1/2 (   q0 w'1 - q3 w'2 + q2 w'3 )  Eq.(*)
d(q2)/dt = 1/2 (   q3 w'1 + q0 w'2 - q1 w'3 )
d(q3)/dt = 1/2 ( - q2 w'1 + q1 w'2 + q0 w'3 )

So,

w'(t), Q(t) -> (*) => dQ(t)/dt.

and we have

Q(t+dt/2) = Q(t) + dt/2 dQ(t)/dt,

Next,

l(t+dt/2) = l(t-dt/2) + dt T(t),

when

l(t+dt/2) => l'(t+dt/2) => w'(t+dt/2),

and we have

w'(t+dt/2), Q(t+dt/2) -> (*) => dQ(t+dt/2)/dt.

And after this, finally

Q(t+dt) = Q(t) + dt dQ(t+dt/2)/dt.

For this algorithm we need to keep l(t-dt/2) for next time step.

Any comments?

==
Best regards, Sergei D.