In this tutorial, we are going to derive the 6 DOF (degree of freedom) equations of motion. Let's get started!
Newton's 2nd Law of motion says the force on an object is proportional to the mass times the acceleration of the object, or
$$ \Sigma \vec{F} = m \vec{a} $$
However, force is also equal to the rate of change of momentum.
[ 1 ]
$$ \Sigma \vec{F} = \frac{d \vec{p}}{dt} = \frac{d}{dt} (m \vec{v}) $$
This can then be split into its x, y, and z components.
$$ F_x = \frac{d}{dt} (m u) $$
$$ F_y = \frac{d}{dt} (m v) $$
$$ F_z = \frac{d}{dt} (m w) $$
Here, u, v, and w are the linear velocities in the x, y, and z directions, respectively.
Fig. 1 - Aircraft split into units of $ \delta m $.
Now let's say we split up the aircraft into infinitesimally small units of mass, $ \delta m $, like in Figure 1, then we can write Eq. (1) as
[ 2 ]
$$ \delta \vec{F} = \delta m \frac{d\vec{v}}{dt} $$
And summing all of the forces would mean $ \Sigma \delta \vec{F} = \vec{F} $. Now looking at a specific point on the plane, $ \delta m $, we could write its velocity as the sum of the plane center of gravity velocity ($ \vec{v_c} $ ) and the change in distance over time between the plane's center of gravity and the point ($ \frac{d\vec{r}}{dt} $). This would give a point velocity of
[ 3 ]
$$ \vec{v} = \vec{v_c} + \frac{d\vec{r}}{dt} $$
on the plane. Eq. (3) can then be substituted into Eq. (2).
$$ \Sigma \delta \vec{F} = \frac{d}{dt} \Sigma (\vec{v_c} + \frac{d\vec{r}}{dt} \delta m) $$
$$ = \frac{d}{dt} \Sigma \vec{v_c} \delta m + \frac{d}{dt} \Sigma \frac{d\vec{r}}{dt} \delta m $$
Now, due to the mass being constant, we can write $ \vec{F} $ as
$$ \vec{F} = m \frac{d\vec{v_c}}{dt} + \frac{d^2}{dt^2} \Sigma \vec{r} \delta m $$
However, since $ \vec{r} $ is being measured from the center of gravity of the plane, $ \Sigma \vec{r} \delta m = 0 $. This leaves us with Newton's 2nd Law reduced to a point mass format.
[ 4 ]
$$ \vec{F} = m \frac{d\vec{v_c}}{dt} $$
Now, remember back to the tutorial Rotating Reference Frames and Chasle's Theorem when we talked about how to account for vector rotation as well as translation in the time rate of change of a vector. We ended up with a theorem called Chasle's Theorem which was defined for a vector $ \vec{A} $ as
[ 5 ]
$$ \frac{d}{dt} \vec{A} = (\frac{d}{dt} \vec{A})_{xyz} + \vec{\omega} \times \vec{A} $$
We can use Eq. (5) to substitute into Eq. (4) to get
[ 6 ]
$$ \vec{F} = m (\frac{d\vec{v_c}}{dt} |_b + \vec{\omega} \times \vec{v_c}) $$
where the subscript b is for the body frame. I should mention that $ \vec{v_c} = \begin{Bmatrix} u \\ v \\ w \end{Bmatrix} $ and $ \vec{\omega} = \begin{Bmatrix} p \\ q \\ r \end{Bmatrix} $. This leaves us with
[ 7 ]
$$ \vec{F} = m \begin{pmatrix} \begin{Bmatrix} \dot{u} \\ \dot{v} \\ \dot{w} \end{Bmatrix} + \begin{Bmatrix} p & q & r \end{Bmatrix} \times \begin{Bmatrix} u \\ v \\ w \end{Bmatrix} \end{pmatrix} $$
The cross product can be computed and then the x, y, and z components of the force can be written out
[ 8 ]
$$ F_x = m(\dot{u} + qw - vr) $$
$$ F_y = m(\dot{v} - pw + ur) $$
$$ F_z = m(\dot{w} + pv - uq) $$
What this means is the acceleration of an aircraft on all axes depends on the force applied along that axis as well as the linear and rotational velocities of the other 2 axes.
Now let's start to look at equations of motion that come from rotation. We can use the same idea with splitting the plane up into smaller mass units with calculating the moment from the rate of change of angular momentum.
[ 9 ]
$$ \delta \vec{M} = \frac{d}{dt} \delta \vec{H} = \frac{d}{dt} (\vec{r} \times \vec{v}) \delta m $$
We can also use Chasle's Theorem like before to obtain the full expression we seek.
$$ \vec{H} = \Sigma \delta \vec{H} = \Sigma (\vec{r} \times \vec{v_c}) \delta m + \Sigma [\vec{r} \times (\vec{\omega} \times \vec{r})] \delta m $$
However, since $ \vec{r} $ is from the center of gravity, $ \Sigma \vec{r} m = 0 $. So,
[ 10 ]
$$ \vec{H} = \Sigma [\vec{r} \times (\vec{\omega} \times \vec{r}] \delta m $$
where $ \vec{r} = \begin{Bmatrix} x \\ y \\ z \end{Bmatrix} $ and $ \vec{\omega} = \begin{Bmatrix} p \\ q \\ r \end{Bmatrix} $. Substituting these into Eq. (10), we can get the components of angular momentum as
[ 11 ]
$$ H_x = p \Sigma (y^2 + z^2) \delta m - q \Sigma xy \delta m - r \Sigma xz \delta m $$
$$ H_y = q \Sigma (x^2 +z^2) \delta m + p \Sigma yx \delta m + r \Sigma yz \delta m $$
$$ H_z = r \Sigma (x^2 + y^2) \delta m - q \Sigma yz \delta m - p \Sigma zx \delta m $$
Now notice that $ I_x = \Sigma (y^2 + z^2) \delta m $, where $ I_x $ is the moment of inertia about the x axis, and $ I_{xy} = \Sigma xy \delta m $. There are 6 different equations for the moments of inertia of a 3D object and they can be substituted into the components of $ \vec{H} $.
[ 12 ]
$$ H_x = p I_x - q I_{xy} - r I_{xz} $$
$$ H_y = -p I_{xy} + q I_y - r I_{yz} $$
$$ H_z = -p I_{xz} - q I_{yz} + r I_z $$
$$ \Rightarrow \vec{H} = I \vec{\omega} $$
If we go back to Eq. (9), we can expand $ \vec{M} = \frac{d \vec{H}}{dt} $ using Chasle's Theorem.
[ 13 ]
$$ \Sigma \vec{M} = \frac{d\vec{H}}{dt} |_b + \vec{\omega} \times \vec{H} $$
$$ \Rightarrow \Sigma \vec{M} = I \dot{\vec{\omega}} + \vec{\omega} \times \vec{H} $$
This leaves us with the components of the moment of inertia, L, M, and N.
[ 14 ]
$$ L = I_x \dot{p} - I_{xz} \dot{r} - I_{xz} pq + (I_z - I_y) qr $$
$$ M = I_y \dot{q} + (I_x - I_z) pr + I_{xz} (p^2 - r^2) $$
$$ L = I_z \dot{r} - I_{xz} \dot{p} + I_{xz} rq + (I_y - I_x) pq $$
We can then rearrange Eq. (13) and solve for $ \dot{\vec{\omega}} $.
[ 15 ]
$$ \dot{\vec{\omega}} = I^{-1} (\Sigma \vec{M} - \vec{\omega} \times I \vec{\omega}) $$
Finally, we can write the equations for $ \dot{\vec{\omega}} $.
[ 16 ]
$$ \dot{p} = (I_zM_x + I_{xz}M_z + I_1pq + I_2qr)/\Delta $$
$$ \dot{q} = (M_y + (I_z - I_x)pr - I_{xz} (p^2 - r^2))/\Delta $$
$$ \dot{r} = (I_xM_z + I_{xz}M_x - I_1qr + I_3pq)/\Delta $$
where $ I_1 = I_{xz}(I_x - I_y + I_z) $, $ I_2 = I_yI_z - I_z^2 + I_{xz}^2 $, $ I_3 = I_{xz}^2 - I_xI_y + I_x^2 $, and $ \Delta = I_xI_z - I_{xz}^2 $.
This was a lot to digest, but to summarize, here are the 6 DOF equations of motion.
Below are the translational equations of motion.
$$ F_x = m(\dot{u} + qw - vr) $$
$$ F_y = m(\dot{v} - pw + ur) $$
$$ F_z = m(\dot{w} + pv - uq) $$
$$ \vec{F} = m(\dot{\vec{v_b}} + \vec{\omega} \times \vec{v_c}) $$
Below are the rotational equations of motion.
$$ L = I_x \dot{p} - I_{xz} \dot{r} - I_{xz} pq + (I_z - I_y) qr $$
$$ M = I_y \dot{q} + (I_x - I_z) pr + I_{xz} (p^2 - r^2) $$
$$ L = I_z \dot{r} - I_{xz} \dot{p} + I_{xz} rq + (I_y - I_x) pq $$
$$ \Sigma \vec{M} = \dot{\vec{H_b}} + \vec{\omega} \times \vec{H} $$