10.0 - Law of Motion

The mechanical laws of motion, along with thermal body attribute updates, are undertaken at cycle point 10.0. The mechanical law of motion will be presented in this section. See PFC Thermal Formulation for details of the thermal formulation.

The motion of a single, rigid particle is determined by the resultant force and moment vectors acting upon it, and can be described in terms of the translational motion of a point in the particle and the rotational motion of the particle. The motion of the center of mass is described in terms of its position \(\mathbf x\), velocity \(\dot{\mathbf x}\) , and acceleration \(\ddot {\mathbf x}\). The rotational motion of the particle is described in terms of its angular velocity \(\pmb \omega\) and angular acceleration \(\dot {\pmb \omega}\).

The equations of motion can be expressed as two vector equations: one relates the resultant force to the translational motion and the other relates the resultant moment to the rotational motion.

Translational Motion

The equation for translational motion can be written in the vector form

(1)\[\mathbf F=m(\ddot{\mathbf x} - \mathbf g)\qquad\hbox{(translational motion)}\]

where \(\mathbf F\) is the resultant force, or the sum of all externally applied forces acting on the particle; \(m\) is the mass of the particle; and \(\mathbf g\) is the body force acceleration vector (e.g., gravitational loading).

The translational equations of motion are solved for balls and clumps via the second order Velocity Verlet algorithm [Verlet1967]. This method of integration offers second order accuracy and is often used in molecular dynamics simulations because, for conservative systems, the energy oscillates around a constant value corresponding to the exact energy of the system.

Suppose that the previous cycle solved Equation (1) at time \(t\) and that the timestep resolved for the current cycle is \(\Delta t\). The 1/2 step velocity, \(\dot{\mathbf{x}}^{(t+\Delta t/2)}\), is calculated as

(2)\[\dot{\mathbf{x}}^{(t+\Delta t/2)} = \dot{\mathbf{x}}^{(t)} + \frac{1}{2} \left(\frac{\mathbf{F}^{(t)}}{m}+\mathbf{g}\right) \Delta t\]

The position at time \(t+\Delta t\) is updated using the 1/2 step velocity

(3)\[\mathbf{x}^{(t+\Delta t)} = \mathbf{x}^{(t)} + \dot{\mathbf{x}}^{(t+\Delta t/2)} \Delta t\]

During the force displacement cycle point, the forces are updated for the current cycle, leading to the updated acceleration \(\ddot{\mathbf{x}}^{(t+\Delta t)}\). The velocity is subsequently updated

(4)\[\dot{\mathbf{x}}^{(t+\Delta t)} = \dot{\mathbf{x}}^{(t+\Delta t/2)} + \frac{1}{2} \left(\frac{\mathbf{F}^{(t+\Delta t)}}{m}+\mathbf{g}\right) \Delta t\]

In PFC, the final velocity update of Equation (4) occurs as the initial step in the timestep determination cycle point or during the finalization stage when finishing a series of cycles as discussed here. Thus, if one queries the velocity of a ball or clump after cycle point 10.0, they are actually querying the 1/2 step velocity \(\dot{\mathbf{x}}^{(t+\Delta t/2)}\).

Rotational Motion

The fundamental equation of rotational motion for a rigid body is

(5)\[\mathbf L = \mathbf I \pmb \omega \qquad\hbox{(rotational motion)}\]

where \(\mathbf L\) is the angular momentum, \(\mathbf I\) is the inertia tensor and \(\pmb \omega\) is the angular velocity. Euler’s equation can be obtained from Equation (5) by taking the time derivative

(6)\[\mathbf M = \dot {\mathbf L} = \mathbf I \dot{\pmb\omega} + \pmb{\omega} \times \mathbf L\]

where \(\mathbf M\) is the resultant moment acting on the rigid body. This relation refers to a local coordinate system that is attached to the body at its center of mass.

Rotational Motion for Balls and 2D Clumps

If the local coordinate system is oriented such that it lies along the principal axes of inertia of the body, then Equation (6) reduces to Euler’s equations of motion:

(7)\[\begin{split}M_1 &= I_1\ \dot\omega_1 + (I_3-I_2)\omega_3\omega_2 \\ M_2 &= I_2\ \dot\omega_2 + (I_1-I_3)\omega_1\omega_3 \\ M_3 &= I_3\ \dot\omega_3 + (I_2-I_1)\omega_2\omega_1 \\\end{split}\]

where \(I_1\), \(I_2\), and \(I_3\) are the principal moments of inertia of the body; \(\dot \omega_1\), \(\dot \omega_2\), and \(\dot \omega_3\) are the angular accelerations about the principal axes; \(\omega_1\), \(\omega_2\), and \(\omega_3\) are the angular velocities about the principal axes; and \(M_1\), \(M_2\), and \(M_3\) are the components of the resultant moment as referred to in the principal axis system.

For a {disk-shaped in 2D; spherical in 3D} body with radius \(R\), whose mass is distributed uniformly throughout its volume, the center of mass coincides with its centroid. For a disk-shaped body whose axis remains in the out-of-plane direction, \(\omega_1 = \omega_2 \equiv 0\), yielding

(8)\[M_3 = I \dot{\omega_3} = (\frac{1}{2}mR^2)\dot{\omega_3}\]

For a spherical body, any local-axis system attached to the center of mass is a principal-axis system, and the three principal moments of inertia are equal to one another. For a spherical body, Equation (7) can be simplified in the global-axis system as

(9)\[\mathbf M = \mathbf I \dot{\pmb{\omega}} = (\frac{2}{5}mR^2)\dot{\pmb{\omega}}\]

For balls (2D and 3D) and 2D clumps, the rotational equations of motion are solved via the Velocity Verlet algorithm described above. Vector notation is used without loss of generality where \(\omega_1 = \omega_2 \equiv 0\) in 2D and the inertia tensor is diagonal. Note that, for clumps in 2D, the polar moment of inertia may be any value that the user specifies.

First the 1/2 step angular velocity is calculated

(10)\[{\pmb{\omega}}^{(t+\Delta t/2)} = {\pmb{\omega}}^{(t)} + \frac{1}{2} \left( \frac{\mathbf{M}^{(t)}}{\mathbf I} \right) \Delta t\]

During the force displacement cycle point, the moments are updated for the current cycle, leading to the updated acceleration \(\dot{\pmb{\omega}}^{(t+\Delta t)}\). The angular velocity is subsequently updated

(11)\[{\pmb{\omega}}^{(t+\Delta t)} = {\pmb{\omega}}^{(t+\Delta t/2)} + \frac{1}{2} \left( \frac{\mathbf{M}^{(t+\Delta t)}}{\mathbf I} \right) \Delta t\]

Rotational Motion for 3D Clumps

3D clumps may have full inertia tensors in the global coordinate system and, as a result, the rotational equations of motion cannot be reduced to Euler’s equations (Equation (7)) in the global coordinate system. In PFC3D, both second- and fourth-order algorithms are available to integrate the clump rotational motion (see the clump order command). In both cases, the angular momentum \(\mathbf L\), the quaternion representing the orientation relative to the principal axis system \(q_c\) and the principal moments of inertia \(\mathbf J\) are stored. An introduction to quaternions is presented prior to a thorough presentation of these second- and fourth-order schemes.


Quaternions are four-dimensional vectors that are commonly used to represent 3D rotations and orientations [Shoemake1985]. In contrast to Euler angle representations, quaternions do not suffer from gimbal lock (i.e., loss of a degree of freedom under certain configurations). A quaternion is defined as

(12)\[q = q_0 + q_i i + q_j j + q_k k = q_0 + \mathbf q\]

where \(i\), \(j\), and \(k\) are the basis elements. \(q_0\) is termed the scalar part, and \(\mathbf q\) is the vector part. The products of basis elements are defined by the identities

(13)\[\begin{split}\begin{array}{c} i^2 = j^2 = k^2 = ijk = -1 \\ ij = k \\ jk = i \\ ki = j \end{array}\end{split}\]

The product of two quaternions, termed the Hamilton product, is determined by the distributive law and the identities above. In scalar/vector notation, the product is given by

(14)\[q p = p_0 q_0 - \mathbf p \cdot \mathbf q + p_0 \mathbf q + q_0 \mathbf p + \mathbf p \times \mathbf q\]

Due to the cross product, the product of two quaternions is not generally commutative (\(q p \neq p q\)). The product of the vector \(\mathbf v\) and a quaternion is performed by creating a pure vector quaternion with 0 scalar part (i.e., \(v = 0 + \mathbf v\) in scalar/vector notation) and performing the multiplication as defined in Equation (14). Addition and subtraction of quaternions is accomplished on a per-component basis.

The conjugate of a quaternion is defined as

(15)\[q^* = q_0 - \mathbf q\]

and the length of a quaternion is given by

(16)\[\|q\| = q_0^2 + q_i^2 + q_j^2 + q_k^2 = q q^* = q^* q\]

A unit quaternion (where \(\|q\| = 1\)) can be converted to a pure rotation matrix. To rotate the vector \(\mathbf v\) by a unit quaternion, one creates the pure vector quaternion \(v=0+\mathbf v\) (with scalar part 0) and forms the product

(17)\[v' = qvq^*\]

The vector component of \(v'\) is the rotated vector. The sequential rotation \(q\) followed by the rotation \(p\) is given by the quaternion product \(p q\), meaning that any number of rotations can by composed into a single quaternion and applied at once, provided they are unit quaternions. Properties of unit quaternions are \(q\) and \(-q\) represent the same rotation, the inverse rotation of \(q\) is the conjugate \(q^*\), and the null rotation is \(q = 1\) where the vector part is \(\mathbf 0\).

Second-Order Solution

[Buss2000] develops algorithms of varying order to solve the rotational equations of motion defined above. Those algorithms are based on the Taylor-series expansion of the angular velocity to provide more accurate estimates of the average angular velocity during a timestep. PFC implements the so-called augmented second-order method. The angular velocity at time \(t\) is calculated by rearranging Equation (5)

(18)\[\pmb \omega^{(t)} = \mathbf{I}^{(t)-1} \mathbf L^{(t)}\]

Note that the inertia tensor (or inverse inertia tensor) in the global axis system can be obtained by rotating the principal moments of inertia by the clump orientation \(q_c\)

(19)\[\begin{split}\begin{array}{c} \mathbf{I} = q_c \mathbf J q_c^* \\ \mathbf{I}^{-1} = q_c \mathbf J^{(-1)} q_c^* \end{array}\end{split}\]

The angular acceleration at time \(t\) is determined by rearranging Equation (6) and using Equation (19)

(20)\[\dot{\pmb \omega}^{(t)} = q_c^{(t)} \mathbf J^{(-1)} q_c^{(t)*} \left(\mathbf M^{(t)}-\mathbf \omega^{(t)} \times \mathbf L^{(t)} \right)\]

and the average angular velocity over the timestep is estimated as

(21)\[\bar{\pmb{\omega}} = \pmb \omega^{(t)} + \dot{\pmb \omega}^{(t)} \frac{\Delta t}{2} + \left(\dot{\pmb{\omega}}^{(t)} \times \pmb \omega^{(t)} \right) \frac{\Delta t^2}{12}\]

The average angular velocity is used to update the clump orientation \(q_c^{(t + \Delta t)}\) by calculating the angular displacement, \(\theta=\bar{\pmb{\omega}}\Delta t\), converting this to a unit quaternion and rotating \(q_c^{(t)}\) by this unit quaternion

(22)\[q_c^{(t + \Delta t)} = \left(cos\left(\frac{\theta}{2}\right) + sin\left(\frac{\theta}{2}\right) \frac{\bar{\pmb{\omega}}}{\|\bar{\pmb{\omega}}\|}\right) q_c^{(t)}\]

Once the clump orientation has been updated, the new angular velocity \(\pmb \omega ^{(t + \Delta t)}\) is given by

(23)\[\pmb \omega ^{(t + \Delta t)} = q_c^{(t+\Delta t)} \mathbf J^{(-1)} q_c^{(t+\Delta t)*} \mathbf L^{(t)}\]

Fourth-Order Solution

One can show that the time derivative of the quaternion is given by

(24)\[\dot q_c = \frac{1}{2} \omega q_c\]

where \(\omega\) is the pure vector quaternion with vector component \(\pmb\omega\) (i.e., the angular velocity). Substituting Equations (18) and (19) into (24) gives:

(25)\[\dot q_c = \frac{1}{2} \left(q_c \mathbf J^{(-1)} q_c^* \mathbf L \right) q_c\]

The standard fourth-order Runge-Kutta method is applied to the equation above. The orientation at time \(t + \Delta t\) is given by

(26)\[q_c^{(t+\Delta t)} = q_c^{(t)} + \frac{\Delta t}{6} \left(k_1+2k_2+2k_3+k_4\right)\]

The quaternion \(k_1\) is calculated as

(27)\[k_1 = \frac{1}{2} \left(q_c^{(t)} \mathbf J^{(-1)} q_c^{(t)*} \mathbf L^{(t)} \right) q_c^{(t)}\]

An intermediate value of the orientation quaternion

(28)\[ {q_{c1}} = q_c^{(t)}+k_1 \frac{\Delta t}{2}\]

is used to calculate the \(k_2\) quaternion

(29)\[k_2 = \frac{1}{2} \left( {q_{c1}} \mathbf J^{(-1)} {q_{c1}^*} \mathbf L^{(t)} \right) {q_{c1}}\]

The next intermediate value of the orientation quaternion

(30)\[ {q_{c2}} = q_c^{(t)}+k_2 \frac{\Delta t}{2}\]

is used to calculate the \(k_3\) quaternion as

(31)\[k_3 = \frac{1}{2} \left( {q_{c2}} \mathbf J^{(-1)} {q_{c2}^*} \mathbf L^{(t)} \right) {q_{c2}}\]

The final intermediate value of the orientation quaternion

(32)\[ {q_{c3}} = q_c^{(t)}+k_3 \Delta t\]

is used to calculate the \(k_4\) quaternion as

(33)\[k_4 = \frac{1}{2} \left( {q_{c3}} \mathbf J^{(-1)} {q_{c3}^*} \mathbf L^{(t)} \right) {q_{c3}}\]

Once the clump orientation has been updated, Equation (23) is used to update the angular velocity to time \(t + \Delta t\). This formulation is similar to the one presented in [Johnson2008].



Johnson, S. M., J. R. Williams and B. K. Cook. “Quaternion-based Rigid Body Rotation Integration Algorithms for use in Particle Methods,” Int. J. Num. Meth. Eng., 74, 1303-1313 (2008).


Buss, S. R. “Accurate and Efficient Simulation of Rigid-Body Rotations,” J. Comp. Phys., 164, 377-406 (2000).


Shoemake, K. “Animating Rotation using Quaternions,” SIGGRAPH85, San Francisco, CA, USA (1985).


Verlet, L. “Computer ‘Experiments’ on Classical Fluids. I. Thermodynamical Properties of Lennard-Jones Molecules,” Phys. Rev., 159, 98-103 (1967).