Motion integration

Motion integration consists in computing accelerations from forces/torques acting on nodes carrying mass and updating their positions/oritentations. Velocities are update based on accelerations, and positions are advanced.

Note

Let us repeat: It must be carefully observed for which time-point is any value valid; this is show by using time-indices \(\prev{x}\equiv x^{(t-\Dt)}\), \(\pprev{x}\equiv x^{(t-\Dt/2)}\), \(\curr{x}\equiv x^{(t)}\), \(\nnext{x} \equiv x^{(t+\Dt/2)}\), \(\next{x} \equiv x^{(t+\Dt)}\). Because of the leap-frog integration scheme, even derivatives of position (position, acceleration, forces, torques) are known at full-steps while odd derivatives (velocity) are known at mid-steps.

Motion integration is implemented in woo.dem.Leapfrog.

The summary force \(\vec{F}_{\Sigma}\) and torque \(\vec{T}_{\Sigma}\) include contributions from all relevant contacts, fields (such as gravity) and internal forces (e.g. from Membrane element).

Todo

Mention how are forces from contacts summed up.

Position

Integrating position consists in using current acceleration \(\curr{\vec{a}}\equiv\curr{\ddot{\vec{u}}}\) on a node to update its position from the current value of position \(\curr{\vec{u}}\) to its value at the next timestep \(\next{\vec{u}}\); we know the node’s velocity \(\pprev{\vec{v}}\equiv\pprev{\dot{\vec{u}}}\). Computation of acceleration, knowing the current force \(\vec{F}_{\Sigma}\) acting on the node in question and its mass \(m\), is simply

\[\curr{\vec{a}}=\frac{\curr{\vec{F}_{\Sigma}}}{m}\]

Using the 2-nd order finite difference approximation with step \(\Dt\), we obtain

\[\curr{\vec{a}}\cong\frac{\prev{\vec{u}}-2\curr{\vec{u}}+\next{\vec{u}}}{\Dt^2}\]

from which we express

\[\begin{split}\next{\vec{u}}&=2\curr{\vec{u}}-\prev{\vec{u}}+\curr{\vec{a}}\Dt^2 =\\ &=\curr{\vec{u}}+\Dt\underbrace{\left(\frac{\curr{\vec{u}}-\prev{\vec{u}}}{\Dt}+\curr{\vec{a}}\Dt\right)}_{(\dagger)}.\end{split}\]

Typically, \(\prev{\vec{u}}\) is already not known (only \(\curr{\vec{u}}\) is); we notice, however, that

\[\pprev{\vec{v}}\simeq\frac{\curr{\vec{u}}-\prev{\vec{u}}}{\Dt},\]

i.e. the mean velocity during the previous step, which is known. Plugging this approximate into the \((\dagger)\) term, we also notice that mean velocity during the current step can be approximated as

\[\nnext{\vec{v}}\simeq\pprev{\vec{v}}+\curr{\vec{a}}\Dt,\]

which is \((\dagger)\); we arrive finally at

\[\next{\vec{u}}=\curr{\vec{u}}+\Dt\left(\pprev{\vec{v}}+\curr{\vec{a}}\Dt\right).\]

The algorithm can then be written down by first computing current mean velocity \(\nnext{\vec{v}}\) which we need to store for the next step (just as we use its old value \(\pprev{\vec{v}}\) now), then computing the position for the next time step \(\next{\vec{u}}\):

(28)\[\begin{split}\nnext{\vec{v}}=\pprev{\vec{v}}+\curr{\vec{a}}\Dt \\ \next{\vec{u}}=\curr{\vec{u}}+\nnext{\vec{v}}\Dt.\end{split}\]

Positions are known at times \(i\Delta t\) (if \(\Delta t\) is constant) while velocities are known at \(i\Delta t+\frac{\Delta t}{2}\). The fact that they interleave (jump over each other) in such way gave rise to the colloquial name “leap-frog”.

Note

The initial (often zero) nodal velocity DemData.vel, before motion is integrated for the first time, is usually understood to be \(\curr{\vec{v}}\) rather than \(\pprev{\vec{v}}\). In that case, (28) is modified to read \(\nnext{\vec{v}}=\curr{\vec{v}}+\curr{\vec{a}}\Dt/2\). The DemData.velMidstep bit is used to know whether the node has already been encountered by the integrator or not. This adjustment can be turned off by unsetting woo.dem.Leapfrog.adjVel0.

Orientation

Spherical nodes

The basic integration procedure applies to nodes which have inertia tensor such that \(\tens{I}_{11}=\tens{I}_{22}=\tens{I}_{33}\) (this tensor is always diagonal, since local node coordinates coincide with principal axes, and is stored as a 3-vector in woo.dem.DemData.inertia; the sphericity can be queried via isAspherical).

\begin{align*} \curr{\dot{\vec{\omega}}}&=\frac{\curr{\vec{T}_{\Sigma}}}{\tens{I}_{11}} \\ \nnext{\vec{\omega}}&=\pprev{\vec{\omega}}+\curr{\dot{\vec{\omega}}}\Dt \\ \next{\vec{q}}&=\mathrm{Quaternion}(\nnext{\vec{\omega}}\Dt)\curr{\vec{q}} \end{align*}

Note

As with linear velocity, the first motion integration modifies the equation to read \(\nnext{\vec{\omega}}=\curr{\vec{\omega}}+\curr{\dot{\vec{\omega}}}\Dt/2\), under the conditions explained above.

Aspherical nodes

Aspherical nodes have different moment of inertia along each principal axis. Their positions are integrated in the same ways as with spherical nodes.

Integrating rotation is considerably more complicated as the local reference frame is not inertial. Rotation of rigid body in the local frame, where inertia tensor \(\mat{I}\) is diagonal, is described in the continuous form by Euler’s equations (\(i\in\{1,2,3\}\), and \(i\), \(j\), \(k\) being sequential):

\[\vec{T}_i=\mat{I}_{ii}\dot{\vec{\omega}}_i+(\mat{I}_{kk}-\mat{I}_{jj})\vec{\omega}_j\vec{\omega}_k.\]

Due to the presence of the current values of both \(\vec{\omega}\) and \(\dot{\vec{\omega}}\), they cannot be solved using the standard leapfrog algorithm.

The algorithm presented here is described by [AT89] (pg. 84–89) and was designed by Fincham for molecular dynamics problems; it is based on extending the leapfrog algorithm by mid-step/on-step estimators of quantities known at on-step/mid-step points in the basic formulation. Although it has received criticism and more precise algorithms are known ([Ome99], [NB06], [JWC08]), this one is currently implemented in Woo for its relative simplicity.

Each node has its local coordinate system aligned with the principal axes of inertia; we use \(\tilde{\bullet}\) to denote vectors in local coordinates. The orientation of the local system is given by the current node orientation \(\curr{q}\) as a quaternion; this quaternion can be expressed as the (current) rotation matrix \(\mat{A}\). Therefore, every vector \(\vec{a}\) is transformed as \(\tilde{\vec{a}}=q\vec{a}q^{*}=\mat{A}\vec{a}\). Since \(\mat{A}\) is a rotation (orthogonal) matrix, the inverse rotation \(\mat{A}^{-1}=\mat{A}^{T}\).

For the node in question, we know

• \(\curr{\tilde{\mat{I}}}\) (constant) inertia tensor diagonal (non-diagonal items are zero, since local coordinates are principal),

• \(\curr{\vec{T}}\) external torque,

• \(\curr{q}\) current orientation ( and its equivalent rotation matrix \(\mat{A}\)),

• \(\pprev{\vec{\omega}}\) mid-step angular velocity,

• \(\pprev{\vec{L}}\) mid-step angular momentum; this is an auxiliary variable that must be tracked additionally for use in this algorithm. It will be zero in the initial step.

Our goal is to compute new values of the latter three, i.e. \(\nnext{\vec{L}}\), \(\next{q}\), \(\nnext{\vec{\omega}}\). We first estimate current angular momentum and compute current local angular velocity:

\begin{align*} \curr{\vec{L}}&=\pprev{\vec{L}}+\curr{\vec{T}}\frac{\Dt}{2}, &\curr{\tilde{\vec{L}}}&=\mat{A}\curr{\vec{L}}, \\ \nnext{\vec{L}}&=\pprev{\vec{L}}+\curr{\vec{T}}\Dt, &\nnext{\tilde{\vec{L}}}&=\mat{A}\nnext{\vec{L}}, \\ \tilde{\curr{\vec{\omega}}}&=\curr{\tilde{\mat{I}}}{}^{-1}\curr{\tilde{\vec{L}}}, \\ \nnext{\tilde{\vec{\omega}}}&=\curr{\tilde{\mat{I}}}{}^{-1}\nnext{\tilde{\vec{L}}}. \\ \end{align*}

Then we compute \(\curr{\dot{q}}\) (see Quaternion differentiation), using \(\curr{q}\) and \(\curr{\tilde{\vec{\omega}}}\):

\begin{align*} \begin{pmatrix}\curr{\dot{q}}_w \\ \curr{\dot{q}}_x \\ \curr{\dot{q}}_y \\ \curr{\dot{q}}_z\end{pmatrix}&= \def\cq{\curr{q}} \frac{1}{2}\begin{pmatrix} \cq_w & -\cq_x & -\cq_y & -\cq_z \\ \cq_x & \cq_w & -\cq_z & \cq_y \\ \cq_y & \cq_z & \cq_w & -\cq_x \\ \cq_z & -\cq_y & \cq_x & \cq_w \end{pmatrix} \begin{pmatrix} 0 \\ \curr{\tilde{\vec{\omega}}}_x \\ \curr{\tilde{\vec{\omega}}}_y \\ \curr{\tilde{\vec{\omega}}}_z \end{pmatrix}, \\ \nnext{q}&=\curr{q}+\curr{\dot{q}}\frac{\Dt}{2}.\\ \end{align*}

We evaluate \(\nnext{\dot{q}}\) from \(\nnext{q}\) and \(\nnext{\tilde{\vec{\omega}}}\) in the same way as in (30) but shifted by \(\Dt/2\) ahead. Then we can finally compute the desired values

\begin{align*} \next{q}&=\curr{q}+\nnext{\dot{q}}\Dt, \\ \nnext{\vec{\omega}}&=\mat{A}^{-1}\nnext{\tilde{\vec{\omega}}}. \end{align*}

Motion in uniformly deforming space

In some simulations, nodes can be considered as moving within uniformly deforming medium with velocity gradient tensor \(\tens{L}\), being stationary at origin. This functionality is provided by the woo.core.Cell class (Scene.cell), along with periodic boundaries. The equations written here only valid when Cell.homoDeform has the value Cell.homoGradV2, which is the default.

Note

It is (currently) not possible to simulate deforming space without periodicity. As workaround, set the periodic space big enough to encompass the whole simulation.

Note

The node’s linear and angular velocities comprise both space and fluctuation velocity. It is (almost) only in the integrator that the two must be distinguished.

Spin is skew-symmetric (rotational) part of velocity gradient \(\tens{W}=\frac{1}{2}(\tens{L}-\tens{L}^T)\); its Hodge dual vector (noted \(*\bullet\)) is position-independent medium angular velocity. Local medium velocities (in cartesian coordinates) read

\[\begin{split}\vec{v}_L&=\tens{L}\vec{x} \\ \vec{\omega}_L&=*\tens{W}\!/2\end{split}\]

where \(\vec{\omega}_L\) is written in cartesian coordinates component-wise as

\[\omega_{Lk}=\frac{1}{2}\epsilon_{ijk}W_{ij}\]

where \(\epsilon_{ijk}\) is the Levi-Civita symbol (also known as “permutation tensor”).

Linear velocity

Velocity \(\vec{v}\) of each particle is sum of fluctuation velocity \(\vec{v}_f\) and local medium velocity \(\vec{v}_L\). \(\vec{v}\) evolves not only by virtue of acceleration, but also of \(\tens{L}\), defined at mid-step time-points. We add this term to (28)

(31)\[\nnext{\vec{v}}=\pprev{\vec{v}}+\Dt\curr{\vec{a}}+\Dt\curr{\vec{\dot v}_L}\]

with

\begin{align*} \Dt\curr{\vec{\dot v}_L}&=\Dt\partial_t(\curr{\tens{L}}\curr{\vec{x}})=\Dt(\curr{\tens{\dot L}}\curr{\vec{x}}+\curr{L}\curr{v})\approx\\ &\approx \Dt \left[\frac{\nnext{\tens{L}}-\pprev{\tens{L}}}{\Dt}\curr{\vec{x}}+\frac{\pprev{\tens{L}}+\nnext{\tens{L}}}{2}\curr{\vec{v}}\right]. \end{align*}

This equation can be rearranged as

\[\Dt\curr{\vec{\dot v}_L}=-\pprev{\tens{L}}(\underbrace{\curr{\vec{x}}-\curr{\vec{v}}\frac{\Dt}{2}}_{\approx\pprev{\vec{x}}})+\nnext{\tens{L}}(\underbrace{\curr{\vec{x}}+\curr{\vec{v}}\frac{\Dt}{2}}_{\approx\nnext{\vec{x}}})\approx-\pprev{\tens{L}}\pprev{\vec{x}}+\nnext{\tens{L}}\nnext{\vec{x}}\]

showing that the correction \(\Dt\curr{\vec{\dot v}_L}\) corresponds to subtracting the previous field velocity and adding the current field velocity.

Going back to (32), we write the unknown on-step velocity as \(\curr{\vec{v}}\approx (\pprev{\vec{v}}+\nnext{\vec{v}})/2\) and substitute \(\Dt\curr{\vec{\dot v}_L}\) into (31) obtaining

\[\nnext{\vec{v}}=\left(\mat{1}-\frac{\nnext{\tens{L}}+\pprev{\tens{L}}}{4}\Dt\right)^{-1}\left[(\nnext{\tens{L}}-\pprev{\tens{L}})\curr{\vec{x}}+\left(\mat{1}+\frac{\nnext{\tens{L}}+\pprev{\tens{L}}}{4}\Dt\right)\pprev{\vec{v}}+\curr{\vec{a}}\Dt\right].\]

The position-independent terms are stored in ImLL4hInv, LmL, IpLL4h and are updated at each timestep. Both \(\pprev{\tens{L}}\) and \(\nnext{\tens{L}}\) must be know in this equations; they are stored in woo.core.Cell.gradV and woo.core.Cell.nextGradV respectively.

Angular velocity

As \(\vec{\omega}_L\) is only a function of time (not of space, unlike \(\vec{v}_L\)), (29) simply becomes

\[\nnext{\vec{\omega}}=\pprev{\vec{\omega}}+\Dt\curr{\dot{\vec{\omega}}} \underbrace{- \pprev{\vec{\omega}_L}+\nnext{\vec{\omega}_L}}_{(\dagger)}\]

where the \((\dagger)\) term is stored as deltaSpinVec.

Warning

Uniformly deforming medium integration is not implemented for aspherical particles, their rotation will be integrated pretending \(\tens{W}\) is zero.

Periodic space

Periodic space is defined as parallelepiped located at origin, expressed as matrix \(\mat{H}\) of which columns are its side vectors. Periodic space volume is given by \(|\det\mat{H}|\).

\(\mat{H}\) follows motion of the space which can be written down as motion of a point \(p\) with zero fluctuation velocity and acceleration; linearizing its motion around mid-step, we obtain its next position as

\[\next{\vec{p}}=\curr{\vec{p}}+\Dt\nnext{\tens{L}}\nnext{\vec{p}}\approx\curr{\vec{p}}+\Dt\nnext{\tens{L}}\frac{\curr{\vec{p}}+\next{\vec{p}}}{2}\]

from which \(\next{\vec{p}}\) is found as

(33)\[\next{\vec{p}}=\left(\mat{1}-\nnext{\tens{L}}\frac{\Dt}{2}\right)^{-1}\left(\mat{1}+\nnext{\tens{L}}\frac{\Dt}{2}\right)\curr{\vec{p}}.\]

The \(\mat{H}\) matrix is updated in the same way, i.e.

\[\next{\mat{H}}=\left(\mat{1}-\nnext{\tens{L}}\frac{\Dt}{2}\right)^{-1}\left(\mat{1}+\nnext{\tens{L}}\frac{\Dt}{2}\right)\curr{\mat{H}}.\]

The transformation matrix accumulates deformation from \(\tens{L}\) and is updated as

\[\next{\mat{T}}=\curr{\mat{T}}+\pprev{\tens{L}}\curr{\mat{T}}\Dt.\]

If \(\tens{L}\) is diagonal, then the (initially zero) \(\mat{T}\) is also diagonal and is a strain matrix. Diagonal contains the Hencky strain along respective axes.

Clumps

Todo

Write about clump motion integration; write about woo.dem.ClumpSpherePack and how its properties are computed. Perhaps move clumps to a separate file.

Timestep

Todo

Write on \(\Dt\), how Scene queries fields and engines on their timestep, about woo.core.Scene.dtSafety, about woo.dem.DynDt.

Numerical damping

Warning

This damping scheme is only suitable for quasi-static simulation where the dynamics is only auxiliary and not important as such. Using this damping scheme in dynamic simulations will lead to gross errors, such as free fall not happening correctly.

In simulations of quasi-static phenomena, it it desirable to dissipate kinetic energy of particles. Since some contact laws (such as Linear (Cundall) contact model) do not include velocity-based damping (such as one in [Addetta2001]), it is possible to use artificial numerical damping. The formulation is described in the manual of the PFC3D code, but the version used here is slightly adapted.

The basic idea is to decrease forces which increase the particle velocities and increase forces which decrease particle velocities by \((\Delta F)_d\), comparing the current acceleration sense and particle velocity sense. This is done component-by-component, which makes the damping scheme clearly non-physical, as it is not invariant with respect to coordinate system rotation; on the other hand, it is very easy to compute. Cundall proposed the form (we omit particle indices $i$ since it applies to all of them separately):

\[\frac{(\Delta \vec{F})_{dw}}{\vec{F}_w}=-\lambda_d\operatorname{sgn}(\vec{F}_w\pprev{\dot{\vec{u}}}_{w}),\quad w\in\{x,y,z\}\]

where \(\lambda_d\) is the damping coefficient. This formulation has several advantages [Hentz2003]:

• it acts on forces (accelerations), not constraining uniform motion;

• it is independent of eigenfrequencies of particles, they will be all damped equally;

• it needs only the dimensionless parameter \(\lambda_d\) which does not have to be scaled.

Woo uses the adapted form

(34)\[\frac{(\Delta\vec{F})_{dw}}{\vec{F}_w}=-\lambda_d\operatorname{sgn}\vec{F}_w\underbrace{\left(\pprev{\dot{u}}_w+\frac{\curr{\ddot{\vec{u}}}_w\Dt}{2}\right)}_{\simeq\curr{\dot{u}}_w},\]

where we replaced the previous mid-step velocity \(\pprev{\dot{u}}\) by its on-step estimate in parentheses. This is to avoid locked-in forces that appear if the velocity changes its sign due to force application at each step, i.e. when the particle in question oscillates around the position of equilibrium with \(2\Dt\) period.

Equation (34) is implemented in Leapfrog; the damping coefficient \(\lambda_d\) is Leapfrog.damping.

Tip

Report issues or inclarities to github.