Equations of motion

The aim of this topic is to explain how the six-degrees of freedom equations of motion with respect to rotating Earth are derived

The flat-Earth aircraft equations of motion are used for short-distance flight only and for observation of the aircraft responses to control inputs. Here the aircraft equations of motion for spherical, rotating Earth are derived. These equations are widely used for aircraft simulation.

Coordinate systems definitions

Derivation

Nomenclature

$[\mathbf{C}]_{b/a}$- Direction cosine matrix transformation from coordinate system $a$ to $b$ $()^{a}$- Vector in coordinate system $a$ $^b\dot{()}^{a}$- Derivative of vector in coordinate system $a$with respect to coordinate system $b$ $()_{b/a}$- $b$ with respect to $a$ $\mu$- latitude $l$- longitude $h$- altitude $Re$- Earth radius

In order to derive the aircraft equations for spherical, rotating Earth four coordinate systems are defined. The Earth-centred inertial (ECI) coordinate system has its origin at the centre of the Earth with its $x_{b}$-axis pointing towards the vernal equinox and its $z_{b}$-axis pointing towards the Earth's north pole. The Earth-fixed centred coordinate system (EFEC) has its origin at the centre of the Earth with its $x_{b}$-axis pointing towards the prime meridian and $z_{b}$-axis towards the Earth's north pole. The north-east-down (NED) coordinate system has its origin at a particular longitude, latitude and altitude with respect to the EFEC coordinate system, its $x_{b}$-axis is pointing north and its $z_{b}$-axis is pointing toward the centre of the Earth. The last coordinate system is the body-fixed coordinate system. Its origin coincides with the NED coordinate system origin.

The direction cosine matrices between the EFEC, NED and body-fixed coordinate systems must be defined. The direction cosine matrix from the EFEC to NED coordinate system in terms of longitude ($\mu$) and latitude ($l$) is:

$[\mathbf{C}]_{f/e}= \begin{bmatrix} -C_{l}S_{\mu} &-S_{l}S_{\mu} &C_{\mu}\\ -S_{l} &C_{l} &0\\ -C_{l}C_{\mu} &-C_{\mu}S_{l} &-S_{\mu}\\ \end{bmatrix}\text{,}$

and the direction cosine matrix from the NED to the body-fixed coordinate system in terms of Euler angles ($\phi$,$\theta$,$\psi$) is:

$[\mathbf{C}]_{b/f}= \begin{bmatrix} C_{\theta}C{\psi} &S_{\phi}S_{\theta}C_{\psi}-C_{\phi}S_{\psi} &C_{\phi}S_{\theta}C_{\psi}+S_{\phi}S_{\psi}\\ C_{\theta}S{\psi} &S_{\phi}S_{\theta}S_{\psi}+C_{\phi}C_{\psi} &C_{\phi}S_{\theta}S_{\psi}-S_{\phi}C_{\psi}\\ -S_{\theta} &S_{\phi}C_{\theta} &C_{\phi}C_{\theta}\end{bmatrix}^T\text{,}$

then the direction cosine matrix from the ECEF to the body-fixed coordinate system is:

$[\mathbf{C}]_{b/e}=[\mathbf{C}]_{b/f}[\mathbf{C}]_{f/e}\text{.}$

A position vector $\mathbf{r}$ in the ECEF coordinate system gives the position of the NED and the body-fixed coordinate systems (the NED and body-fixed coordinate systems coincide, as mentioned earlier). The time derivative of $\mathbf{r}$with respect to the ECI coordinate system is given by:

$^i\dot{\mathbf{r}}^{e}=^e\dot{\mathbf{r}}^{e}+\mathbf{\omega}_{e/i}^{e}\times\mathbf{r}^{e}\text{,}$

where $\mathbf{\omega}_{e/i}$is the earth rotational velocity. Differentiating again with respect to time gives:

$^i\ddot{\mathbf{r}}^{e}=^e\dot{\mathbf{r}}^{e}+^e\dot{\mathbf{\omega}}_{e/i}^{e}\times\mathbf{r}^{e}+\mathbf{\omega}_{e/i}^{e}\times ^e\dot{\mathbf{r}}^{e}+\mathbf{\omega}_{e/i}^{e}\times ^e\dot{\mathbf{r}}^{e}+\mathbf{\omega}_{e/i}^{e}\times\mathbf{\omega}_{e/i}^{e}\times\mathbf{r}^{e}\text{.}$

Simplifying by gathering terms gives:

$^i\ddot{\mathbf{r}}^{e}=^e\ddot{\mathbf{r}}^{e}+^e\dot{\mathbf{\omega}}_{e/i}^{e}\times\mathbf{r}^{e}+2\mathbf{\omega}_{e/i}^{e}\times ^e\dot{\mathbf{r}}^{e}+\mathbf{\omega}_{e/i}^{e}\times\mathbf{\omega}_{e/i}^{e}\times\mathbf{r}^{e}\text{.}$

Using Newton's Second Law (often stated as $F=ma$) gives:

$^i\ddot{\mathbf{r}}^{e}=\mathbf{a}^{e}=\frac{\mathbf{F}^{e}}{m}+\mathbf{G}^{e}\text{,}$

where $\mathbf{F}^{e}$is the force vector and $\mathbf{G}^{e}$is the gravitational vector in ECEF coordinates.

Some of the derivatives can be rewritten as:

$^e\ddot{\mathbf{r}}^{e}=^e\dot{\mathbf{V}}^{e}\text{,}\\ ^e\dot{\mathbf{r}}^{e}=\mathbf{V}^{e}\text{.}$

If the above definitions are substituted in the acceleration equation the following equation is obtained:

$\frac{\mathbf{F}^{e}}{m}+\mathbf{G}^{e}=^e\dot{\mathbf{V}}^{e}+^e\dot{\mathbf{\omega}}_{e/i}^{e}\times\mathbf{r}^{e}+2\mathbf{\omega}_{e/i}^{e}\times \mathbf{V}^{e}+\mathbf{\omega}_{e/i}^{e}\times\mathbf{\omega}_{e/i}^{e}\times\mathbf{r}^{e}\text{.}$

Assuming the Earth's rotational velocity is constant:

$^e\dot{\mathbf{\omega}}_{e/i}^{e}\times\mathbf{r}^{e}=0\text{,}$

and multiplying by the direction cosine matrix $[\mathbf{C}]_{b/e}$ gives the following equation:

$\frac{\mathbf{F}^{b}}{m}+\mathbf{G}^{b}=^e\dot{\mathbf{V}}^{b}+2[\mathbf{C}]_{b/e}\mathbf{\omega}_{e/i}^{e}\times \mathbf{V}^{b}+[\mathbf{C}]_{b/e}\mathbf{\omega}_{e/i}^{e}\times\mathbf{\omega}_{e/i}^{e}\times\mathbf{r}^{e}\text{.}$

Rearranging for $^e\dot{\mathbf{V}}^{b}$ gives:

$^e\dot{\mathbf{V}}^{b}=\frac{\mathbf{F}^{b}}{m}-2[\mathbf{C}]_{b/e}\mathbf{\omega}_{e/i}^{e}\times \mathbf{V}^{b}+\mathbf{G}^{b}-[\mathbf{C}]_{b/e}\mathbf{\omega}_{e/i}^{e}\times\mathbf{\omega}_{e/i}^{e}\times\mathbf{r}^{e}\text{,}$

where $\mathbf{r}^{e}$is given by:

$\mathbf{r}^{e}=\begin{Bmatrix}(Re+h)cos(\mu)cos(l)\\(Re+h)cos(\mu)sin(l)\\(Re+h)sin(\mu)\end{Bmatrix}\text{.}$

This equation gives the derivative of the velocity vector in the body-fixed coordinate system with respect to the EFEC coordinate system. However, since the body-fixed coordinate system rotates with respect to the EFEC coordinate system the following equation is used:

$^e\dot{\mathbf{V}}^{b}=^b\dot{\mathbf{V}}^{b}+\mathbf{\omega}_{b/e}^{b}\times\mathbf{V}^{b}\text{.}$

Substitution of the above equation in the previous gives:

$^b\dot{\mathbf{V}}^{b}=\frac{\mathbf{F}^{b}}{m}-\mathbf{\omega}_{b/e}^{b}\times\mathbf{V}^{b}-2[\mathbf{C}]_{b/e}\mathbf{\omega}_{e/i}^{e}\times \mathbf{V}^{b}+\mathbf{G}^{b}-[\mathbf{C}]_{b/e}\mathbf{\omega}_{e/i}^{e}\times\mathbf{\omega}_{e/i}^{e}\times\mathbf{r}^{e} \text{.}$

The above equation gives the derivative of the velocity vector in the body-fixed coordinate system with respect to the body-fixed coordinate system. The flat-Earth aircraft equations of motion can be identified. However, there are additional terms for the Coriolis and Centripetal acceleration due to the Earth's rotational velocity. Since angular velocities are additive, the angular velocity of the body-fixed coordinate system relative to the ECI the coordinate system in body coordinates is given by:

$\mathbf{\omega}_{b/i}^{b}= \mathbf{\omega}_{b/e}^{b}+ [\mathbf{C}]_{b/e}\mathbf{\omega}_{e/i}^{e}\text{.}$

Rearranging the above equation for $\mathbf{\omega}_{b/e}^{b}$and substituting in the previous equation gives:

$^b\dot{\mathbf{V}}^{b}=\frac{\mathbf{F}^{b}}{m}-(\mathbf{\omega}_{b/i}^{b}-[\mathbf{C}]_{b/e}\mathbf{\omega}_{e/i}^{e})\times\mathbf{V}^{b}-2[\mathbf{C}]_{b/e}\mathbf{\omega}_{e/i}^{e}\times \mathbf{V}^{b}+\mathbf{G}^{b}-[\mathbf{C}]_{b/e}\mathbf{\omega}_{e/i}^{e}\times\mathbf{\omega}_{e/i}^{e}\times\mathbf{r}^{e}\text{.}$

Further simplifying the above equation gives:

$^b\dot{\mathbf{V}}^{b}=\frac{\mathbf{F}^{b}}{m}-\mathbf{\omega}_{b/i}^{b}\times\mathbf{V}^{b}-[\mathbf{C}]_{b/e}\mathbf{\omega}_{e/i}^{e}\times \mathbf{V}^{b}+\mathbf{G}^{b}-[\mathbf{C}]_{b/e}\mathbf{\omega}_{e/i}^{e}\times\mathbf{\omega}_{e/i}^{e}\times\mathbf{r}^{e}\text{.}$

The above are the translational aircraft equations of motion. The rotational aircraft equations of motion are given by:

$\dot{\mathbf{\omega}}_{b/i}^{b}=[\mathbf{I}]^{-1}(\mathbf{M}-\mathbf{\omega}_{b/i}^{b}\times([\mathbf{I}]\mathbf{\omega}_{b/i}^{b}))\text{,}$

and are similar to the flat-Earth rotational equations of motion with the difference that $\mathbf{\omega}_{b/i}^{b}$ is the aircraft angular velocity vector with respect to the ECI coordinate system. The components of this angular velocity vector are not directly equal to $p$, $q$, and $r$ which is the case for the flat-Earth equations of motion. Here they are given by $\mathbf{\omega}_{b/f}^{b}$. In order to determine $p$, $q$ and $r$, the rotation of the NED with respect to the EFEC coordinate system must be taken into account. Why is rotation present between the NED and EFEC coordinate systems? The NED coordinate system is always tangent to the Earth. If the aircraft starts its journey at $l=0$ deg and ends its journey at $l=90$ deg, the NED coordinate system would have rotated in order to stay tangent to the Earth. So the angular velocity of the body-fixed to the ECEF coordinate system can be represented as:

$\mathbf{\omega}_{b/e}^{b}=\mathbf{\omega}_{b/f}^{b}+[\mathbf{C}]_{b/f}\mathbf{\omega}_{f/e}^{f}\text{,}$

where $\mathbf{\omega}_{b/f}^{b}$is the angular velocity of the body-fixed with respect to the NED coordinate system and $\mathbf{\omega}_{f/e}^{f}$is the angular velocity of the NED with respect to the ECEF coordinate system.

The angular velocity of the NED coordinate system with respect to the ECEF coordinate system is given by:

$\mathbf{\omega}_{f/e}^{f}=\begin{Bmatrix}\dot{l}cos(\mu)\\-\dot{\mu}\\-\dot{l}sin(\mu)\end{Bmatrix}\text{,}$

where

$\dot{\mu}=\frac{1}{Re+h}V_{f_{x}}\text{,}\\ \dot{l}=\frac{1}{Re+h}cos(\mu)V_{f_{y}}\text{,}\\ \dot{h}=-V_{f_{z}}\text{,}\\ \mathbf{V_{f}}=[\mathbf{C}]_{b/f}^T\mathbf{V_{b}}\text{.}$

Then the angular velocity of the body-fixed relative to the NED coordinate system is given by:

$\mathbf{\omega}_{b/f}^{b}=\mathbf{\omega}_{b/e}^{b}-[\mathbf{C}]_{b/f}\mathbf{\omega}_{f/e}^{f}\text{,}$

which can be further simplified if using the previously derived definition for $\mathbf{\omega}_{b/e}^{b}$:

$\mathbf{\omega}_{b/f}^{b}=\mathbf{\omega}_{b/i}^{b}-[\mathbf{C}]_{b/f}\mathbf{\omega}_{f/e}^{f}-[\mathbf{C}]_{b/e}\mathbf{\omega}_{e/i}^{e}\text{.}$

Non-linear state-space form

$\mathbf{V^{b}}= \begin{Bmatrix} u\\ v\\ w \end{Bmatrix}, \mathbf{\omega}_{b/f}= \begin{Bmatrix} p\\ q\\ r \end{Bmatrix}\\\\ ^b\dot{\mathbf{V}}^{b}=\begin{Bmatrix} \dot{u}_{x}\\ \dot{v}_{y}\\ \dot{w}_{z}\\ \end{Bmatrix}=\frac{\mathbf{F}^{b}}{m}-\mathbf{\omega}_{b/i}^{b}\times\mathbf{V}^{b}-[\mathbf{C}]_{b/e}\mathbf{\omega}_{e/i}^{e}\times \mathbf{V}^{b}+\mathbf{G}^{b}-[\mathbf{C}]_{b/e}\mathbf{\omega}_{e/i}^{e}\times\mathbf{\omega}_{e/i}^{e}\times\mathbf{r}^{e}\\\dot{\mathbf{\omega}}_{b/i}^{b}=\begin{Bmatrix} \dot{\omega}_{x}\\ \dot{\omega}_{y}\\ \dot{\omega}_{z}\\ \end{Bmatrix}=[\mathbf{I}]^{-1}(\mathbf{M}-\mathbf{\omega}_{b/i}^{b}\times([\mathbf{I}]\mathbf{\omega}_{b/i}^{b}))\\ \mathbf{\omega}_{b/f}^{b}=\mathbf{\omega}_{b/i}^{b}-[\mathbf{C}]_{b/f}\mathbf{\omega}_{f/e}^{f}-[\mathbf{C}]_{b/e}\mathbf{\omega}_{e/i}^{e}$

$\dot{\mathbf{x}}= \begin{Bmatrix} \dot{x}_{1} \\ \dot{x}_{2} \\ \dot{x}_{3} \\ \dot{x}_{4} \\ \dot{x}_{5} \\ \dot{x}_{6} \\ \dot{x}_{7} \\ \dot{x}_{8} \\ \dot{x}_{9} \\ \dot{x}_{10} \\ \dot{x}_{11} \\ \vdots \\ \dot{x}_{18} \\ \dot{x}_{19} \\ \dot{x}_{20} \\ \vdots \\ \dot{x}_{27} \\ \end{Bmatrix} =\begin{Bmatrix} \dot{u}_{x}\\ \dot{v}_{y}\\ \dot{w}_{z}\\ \dot{\omega}_{x}\\ \dot{\omega}_{y}\\ \dot{\omega}_{z}\\ \dot{\mu}\\ \dot{l}\\ \dot{h}\\ \dot{C11}_{f/e}\\ \dot{C12}_{f/e}\\ \vdots \\ \dot{C33}_{f/e}\\ \dot{C11}_{b/f}\\ \dot{C12}_{b/f}\\ \vdots \\ \dot{C33}_{b/f}\\ \end{Bmatrix}$