Differential Equations for Rigid Aircraft Flight Dynamics
(Equations of Motion)
On this page we shall derive the equations of motion needed to fully describe the motion of an aircraft in response to forces and moments and describe in broad terms how to solve them numerically on a computer. In doing so, we shall draw on some of the results we have obtained already in the article “Aircraft Attitude and Angles”, which discusses kinematics, i.e. the geometry of the motion (without any physics, i.e. without forces and moments).
The reader shall note that “dynamics” and “kinetics” are used interchangeably, with the former being more popular in recent literature. Both terms refer to the physics of motion, which has to do with forces, moments, mass, inertia tensor, etc. In contrast, the term “kinematics” refers to the geometry of motion (without any physics, and answer questions such as given the body angular rates of the aircraft, how do the values of the Euler angles change with time.
Derivation of Equations for Rigid Aircraft Flight Dynamics
Starting Point
We will need four equations to describe aircraft motion completely:
- two dynamic (kinetic) equations which capture the physics, and
- two kinematic equations, which capture the geometry of the problem.
This article is primarily about the first two, since the last two were already derived in the article “Aircraft Attitude and Euler Angles”.
The four equations are:
- Translational Dynamic Equation:
\begin{equation}
\mathbf{F} = \frac{d\mathbf{P}}{dt},
\end{equation}
where \(\mathbf{F}\) represents the sum of all total external forces (from all sources, including gravity), \(\mathbf{P}:=m\mathbf{V}\) is the translational momentum of the aircraft (with \(m\) being the mass of the aircraft – in general time dependent because of fuel burn, though we shall ignore this in subsequent calculations -, and \(\mathbf{V}\) being the velocity vector of the aircraft). - Rotational Dynamic Equation:
\begin{equation}
\mathbf{M} = \frac{d\mathbf{H}}{dt},
\end{equation}
where \(M\) is the total external moment vector, and \(\mathbf{H}:=([I]\underline{\omega})\) is the angular momentum of the aircraft (with \([I]\) being the inertia tensor – a 3×3 matrix when written in components with respect to a basis, which plays the role of “rotational mass” and captures the 3-dimensional mass distribution of the aircraft -, and \(\underline{\omega}\) being the angular momentum vector. Note that unlike in translational motion, where velocity \(\mathbf{V}\) and translational momentum \(\mathbf{P}\) are always alined, because \(m\) is a scalar (a single number, not a vector or matrix), in translational motion in general the angular velocity vector \(\underline{\omega}\) and the angular momentum vector \(\mathbf{H}\) are not alined, because the inertia tensor \([I]\) (being a matrix and not a scalar) mixes the components of \(\underline{\omega}\) into \(\mathbf{H}\). - Translational Kinematic Equation:
\begin{equation}
\frac{{}^Ed\mathbf{r}_E}{dt} = \mathbf{V}_E = [EB] \mathbf{V}_B.
\end{equation}
This equation translates the velocity vector in the body frame, \(\mathbf{V}_B\) to the world frame, \(\mathbf{V}_B\), and then expresses it as the first time derivative of position \(\mathbf{r}_E\). \([EB]\) is the directional cosine matrix (DCM) which takes vector components from the body frame into the world frame (see our article “Aircraft Attitude and Euler Angles”). - Rotational Kinematic Equation:
\begin{equation}
\frac{{}^{E}\underline{\Theta}}{dt} = [B]\underline{\omega}_B
\end{equation}
This equation relates the body angular rates \(\underline{\omega}_B:=(p, q, r)\) to the Euler rates (first time derivative of the Euler angles) \(\underline{\dot\Theta}=(\dot\phi, \dot\theta, \dot\psi)\). \([B]\) is the body-to-Euler-rate transformation matrix, which we have introduced and explained in the article “Aircraft Attitudes and Euler Angles”; it reflects the fact that the body rates are not the same as the Euler rates, because the body rates are expressed in the body frame, while the Euler angles have been constructed in a mixture of three frames.
With the four equations above we can compute the aircraft position and attitude at any later time, if we know for the initial time:
- Aircraft Position
- Velocity Vector of the Aircraft
- Aircraft Attitude (Euler Angles)
- Angular Velocity of the Aircraft (Body Rates)
This amounts to knowing the 12-dimensional state vector \((\mathbf{r}_E, \mathbf{V}_B, \underline{\Theta}, \underline{\omega}_B)\) or, written in components, \((x, y, z; u, v, w; \phi, \theta, \psi; p, q, r)\).
Origin of the Equations
The above equations are written in a very general form which does not let us calculate anything yet. The first two equations express the conservation of translational and angular momentum in the absence of any forces and moments, as well as how external forces and moments change translational and angular momentum.
The third and fourth equations (kinematic equations) relate that velocities in the body frame can be expressed in the inertial frame, and that body rates can be converted to Euler rates with a linear transformation. Luckily, we have already introduced the matrices [EB] and [B] that accomplish this feat in article “Aircraft Attitude and Euler Angles”, so these equations have already been covered there, and the only thing that remains for us to do it to derive the dynamic equations (the first two equations).
Intermezzo: Time Derivative of Vector-Valued Functions in a Rotating Frame
Before we can proceed, we must learn how to take the time derivative in a rotating reference frame. Let us have two reference frames A and B, where frame B is rotating with respect to the first with angular velocity vector \(\underline{\omega}_{B/A}\). Then we obtain the following rules how to do time derivatives of functions in the rotating frame B, denoted by \(\frac{{}^B d}{dt}\).
For a scalar function \(f\), i.e. a function from the space of real numbers \(R\) into the space of real numbers \(R\),
\begin{eqnarray}
f: R &\rightarrow& R\\
t &\mapsto& f(t)
\end{eqnarray}
we have:
\begin{equation}
\frac{{}^Adf}{dt} = \frac{{}^Bdf}{dt}.
\end{equation}
In other words, the time derivative of a scalar function \(f\) is the same if taken in frame A or frame B. An example of this is a temperature gauge, which will show the same number in both the inertial and the body frame of an airplane.
For vector-valued functions, i.e. a function \(\mathbf{f}\) from \(R\) into a 3-dimensional vector space \( V \) with
\begin{eqnarray}
\mathbf{f}: R &\rightarrow& V\\
t &\mapsto& \mathbf{f}(t)=f_1(t){\mathbf{v}}_1(t)+ f_2(t){\mathbf{v}}_2(t)+ f_3(t){\mathbf{v}}_3(t)
\end{eqnarray}
the relationship is more complicated. This is because not only are the scalar coefficients \(f_i(t)\), \(i=1,2,3\), time dependent, but the linearly independent basis vectors \({\mathbf{v}}_i(t)\) spanning \(V\) are time dependent as well, because the basis is rotating. The product rule of calculus applies, and the equation obtains an additional term, coming from the time dependence of the basis vectors.
For vector valued functions \(\mathbf{f}\), we have (we give this result here without derivation):
\begin{equation}\label{eq:df}
\frac{{}^A d\mathbf{f}}{dt} = \frac{{}^B d\mathbf{f}}{dt} + (\underline{\omega}_{B/A}\times \mathbf{f}),
\end{equation}
where the symbol \(\times\) in the last term denotes the cross product (or vector product) of two vectors, and \(\underline{\omega}_{B/A}\) denotes the angular velocity with which reference frame B is rotating with respect to reference frame A.
In our application for the derivation of the equations of motion, we shall use the inertial/world/Earth frame E for frame A above, and the body (aircraft) frame B for frame B above. The angular velocity \(\underline{\omega}_{B/A}\)between the two then simply is the angular velocity \(\underline{\omega}_{B}\) of the aircraft.
Let us give two examples to illustrate the equation above: 1) a needle fixed in the aircraft cockpit pointing at the right wing. The first term on the right will be zero, because the needle does not move in the body frame of the aircraft. However, it does move in the inertial frame as the aircraft moves. This is captured by the second term on the right, as well as the non-zero term on the left. 2) A 3D compass needle always pointing north. In this example, the term on the left-hand side is zero, because the needle does not rotate in the inertial frame of the world. However, it does rotate in the aircraft frame, as the aircraft moves, which is captured by the first term on the right. The second term on the right exactly cancels the first term in this particular case, such that the equation holds and gives the proper zero answer in the inertial frame.
This time derivative expression for rotating frames will allow us to write the dynamic equations in the body frame. In doing so, they will acquire an addition term containing the cross product with the angular velocity. A huge benefit will be however, that the inertial tensor will become time independent in the body frame (neglecting aeroelasticity).
Because the cross product is mathematically cumbersome to deal with, the way it mixes the vector components, we shall employ a mathematical trick and replace it by matrix multiplication, using the cross-product-equivalent matrix
\begin{equation}[\underline{\tilde \omega}]
:=\begin{pmatrix}
0 & -\omega_3 & \omega_2\\
\omega_3 & 0 & -\omega_1\\
-\omega_2 & \omega_1 & 0
\end{pmatrix}.
\end{equation}
With this matrix, we can turn the cross product into a standard matrix multiplication, using the mathematical equality:
\begin{equation}
\underline\omega \times \mathbf{f} = [\underline{\tilde \omega}] \mathbf{f}.
\end{equation}
That way we can apply all the standard rules of matrix manipulations to this term as well, if needed.
Derivation of the Translational Dynamic Equation
In Vectorial Form
With the above, the translational dynamic equation of motion in the world frame
\begin{equation}
\mathbf{F}_E = \frac{{}^E d(m\mathbf{V}_E)}{dt}
\end{equation}
can be written in vectorial form in the body frame as (we here assume mass is independent of time)
\begin{eqnarray}
\mathbf{F}_B &=& m\frac{{}^B d\mathbf{V}_B}{dt} + m\underline{\omega}_B\times \mathbf{V}_B = \\ &=& m\frac{{}^B d\mathbf{V}_B}{dt} + m[\underline{\tilde \omega}_B]\mathbf{V}_B
\end{eqnarray}
or, solved for acceleration in the body frame, as:
\begin{eqnarray}
\frac{{}^B \mathbf{V}_B}{dt} &=& \frac{\mathbf{F}_B}{m} + \underline{ \omega}_B \times\mathbf{V}_B = \\ &=& \frac{\mathbf{F}_B}{m} + [\underline{\tilde \omega}_B]\mathbf{V}_B.
\end{eqnarray}
In Scalar Form (Component Form)
In scalar form (component form) the above equation is actually a system of three coupled differential equations:
\begin{eqnarray}
F_{x_B} &=& m (\dot u + qw -rv)\\
F_{y_B} &=& m(\dot v + ru -pw)\\
F_{z_B} &=& m(\dot w + pv – qu),
\end{eqnarray}
where we have used \(\mathbf{F}_B=(F_{x_B}, F_{y_B}, F_{z_B})^T\) for the force components in the body frame and \(\mathbf{V}_B=(u, v, w)^T\) for the velocity vector in the body frame.
The above system of equations can be solved for the acceleration in the body frame:
\begin{eqnarray}
\dot u &=& \frac{F_{x_B}}{m} – qw + rv\\
\dot v &=& \frac{F_{y_B}}{m} – ru + pw\\
\dot w &=& \frac{F_{z_B}}{m} – pv + qu.
\end{eqnarray}
It is in this form that they will be needed to find a solution by numerical integration.
For the total force vector in the body frame, \(\mathbf{F}_B\), we have \(\mathbf{F}_B = \mathbf{F}_{B_\mathrm{aero}}+\mathbf{F}_{B_\mathrm{thrust}}+ m[BE]\mathbf{g}_E\), with gravitational acceleration \(\mathbf{g}_E=(0, 0, g)^T\) in the inertial frame. With the expression for [BE] derived in the article “Aircraft Attitude and Euler Angles”, we obtain \(\mathbf{g}_B = (-g \sin \theta, g \sin \phi \cos \theta, g \cos \phi \cos \theta)^T\), and therefore:
\begin{eqnarray}
\dot u & = & \frac{\mathbf{F}_{x_{B_\mathrm{aero}}}+\mathbf{F}_{x_{B_\mathrm{thrust}}}}{m} – g \sin\theta – qw + rv \\
\dot v & = & \frac{\mathbf{F}_{y_{B_\mathrm{aero}}}+\mathbf{F}_{y_{B_\mathrm{thrust}}}}{m} + g \cos\phi\cos\theta – ru + pw \\
\dot w & = & \frac{\mathbf{F}_{z_{B_\mathrm{aero}}}+\mathbf{F}_{z_{B_\mathrm{thrust}}}}{m} + g \cos\phi\cos\theta – pv + qu
\end{eqnarray}
\(\underline{\omega}_B=(p, q, r)) for the angular velocity.
Derivation of the Rotational Dynamic Equation
In Vectorial Form
The derivation of the rotational equation of motion follows similar lines. Using the above formula for time derivatives of vector valued functions, we get:
\begin{equation}
\frac{{}^E d([I]\underline{\omega})}{dt} = \frac{{}^B d([I]\underline{\omega})}{dt} + \underline{\omega}_{B/E} \times ([I]\underline{\omega}) = \underbrace{\frac{{}^B d[I]}{dt}}_{=0}\underline{\omega} + [I]\frac{{}^B d\underline{\omega}}{dt} + {\underline{\omega}_{B/E}} \times ([I]\underline{\omega}),
\end{equation}
where we have indicated in the last step that the time derivative of the inertia tensor \([I]\) will be zero, because the body frame moves with the aircraft, and the mass distribution therefore does not move in the body frame (if we ignore aeroelasticity, i.e. the deformation of the aircraft structure due to aerodynamic forces).
The above expression transform the right hand side of the rotational dynamic equation into the body frame. We thus obtain:
\begin{equation}
\mathbf{M}_B = [I_B] \frac{{}^B d\underline{\omega}_B}{dt} + \underline{\omega}_B \times [I_B] \underline{\omega}_B
= [I_B] \frac{{}^B d\underline{\omega}_B}{dt} + [\underline{\tilde\omega}_B] [I_B] \underline{\omega}_B
\end{equation}
This is Euler’s equation for rigid body dynamics.
Solved for angular acceleration the equation becomes:
\begin{equation}
\frac{d\underline{\omega}_B}{dt} = [I_B]^{-1} ( \mathbf{M}_B -[\underline{\tilde\omega}_B] [I_B] \underline{\omega}_B )
\end{equation}
Scalar Form (Component Form)
Rigid Body Equations of Motion for a Left-Right Symmetric Airplane
Notice the multiplicative inverse of the inertial moment tensor, \([I_B]^{-1}\), on the right-hand side in the last equation above. When writing this equation in scalar form (in components), this would make the system of equations very complicated to write down. While it can be done, we shall instead make a common simplifying assumption before giving the equations in component form: we shall assume that we are dealing with a left-right symmetric aircraft (most airplanes are, though there are notable exceptions). The following entries of the inertia tensor matrix vanish, \(I_{xy}=I_{yx}=I_{yz}=I_{zy}=0\), for aircraft with left-right mirror symmetry.
The system of equations then takes the following form:
\begin{eqnarray}
\mathcal{L} &=& \dot p I_{xx} + qr (I_{zz} – I_{yy}) – (\dot r + pq) I_{xz}\\
\mathcal{M} &=& \dot q I_{yy} – pr (I_{zz} – I_{xx}) + (p^2 + r^2) I_{xz}\\
\mathcal{N} &=& \dot r I_{zz} + pq (I_{yy}-I_{xx}) + (qr-\dot p) I_{xz},
\end{eqnarray}
where we have used \(\mathbf{M}_B=(\mathcal{L}, \mathcal{M}, \mathcal{N})^T\) for the components of the total moment vector in the body frame.
Solving for the angular accelerations, we obtain:
\begin{eqnarray}
\dot p &=& \frac{I_{zz}\mathcal{L} + I_{xz}\mathcal{N} – \{ I_{xz}(I_{yy}-I_{xx}-I_{zz})p + [I_{xz}^2 + I_{zz}(I_{zz}-I_{yy})]r\}q}{I_{xx}I_{zz}-I_{xz}^2}\\
\dot q &=& \frac{\mathcal{M} – (I_{xx}-I_{zz})pr – I_{xz})p^2-r^2)}{I_{yy}}\\
\dot r &=& \frac{I_{xz}\mathcal{L}+I_{xx}\mathcal{N}+\{I_{xz}(I_{yy}-I_{xx}-I_{zz})r+[I_{xz}^2+I_{xx}(I_{xx}-I_{yy})]p\}q}{I_{xx}I_{zz}-I_{xz}^2}
\end{eqnarray}
This is the form in which we need the equations in order to find a numerical solution.
Simplified Form of Euler's Equations
The \(I_{xz}\) component is typically not zero. But it is usually smaller than the other three, and if we neglect it, the three equations above (Eqs. \eqref{eq:pdot}, \eqref{eq:qdot}, and \eqref{eq:rdot}) take on the even simpler form:
\begin{eqnarray}
I_{xx} \dot p & = & \mathcal{L} + (I_{zz} – I_{yy})qr \label{eq:pdot simple}\\
I_{yy} \dot q & = & \mathcal{M} + (I_{zz} – I_{xx})pr \label{eq:qdot simple}\\
I_{zz} \dot r & = & \mathcal{N} + (I_{xx} – I_{yy})qp \label{eq:rdot simple}
\end{eqnarray}
Setting \(I_{xz}\) to zero eliminated a lot of the terms. In this form the equations are often referred to as the “simplified” form of Euler’s equations. (Note that here \(I_{xx}\), \(I_{yy}\), and \(I_{zz}\) are scalar quantities and can simply be taken to the other side by regular division (no matrix inversion necessary).
In this simplified form, the interpretation of these three equations also becomes more tractable: \(I_{zz}>I_{xx}\) and \(I_{zz}>I_{yy}\) for all fixed-wing aircraft, because both wings and fuselage contribute to \(I_{zz}\), but primarily only the wing contributes to \(I_{xx}\) and only the fuselage to \(I_{zz}\). From this we can infer from the first two equations, respectively, that pitching up and yawing to the right induces a rolling moment to the right (Eq. \eqref{eq:pdot simple}, while rolling to the right and yawing to the right induces a pitch up moment (Eq. \eqref{eq:qdot simple}). The latter can be used, for instance, for a roll-coupling entry into an inverted spin in aircraft which do not have sufficient down elevator authority to induce an inverted spin otherwise. This third equation, Eq. \eqref{eq:rdot simple} shows us that pitching up and simultaneously rolling to the right can induce either a yaw moment to the left or to the right: it depends on whether \(I_{xx}\gt I_{yy}\) or \(I_{xx}\lt I_{yy}\), which depends on the airplane make and model. If the yaw moment is to the left, as is the case in fuselage-heavy airplanes, it can be used as an anti-spin input (in an upright spin both yaw and roll are in the same direction): the pilot does so by providing in-spin aileron (towards the spin axis); while this increases the roll rate, it also slows down the yaw rate, and is in some aircraft the only way to perform a recovery. This is why the spin recovery procedure in some business jets and military aircraft (e.g. the McDonnell Douglas F-4 Phantom) calls for in-spin aileron to slow down the yaw rate.
Numerical Solution of Equations for Rigid Aircraft Flight Dynamics
This section describes the general numerical integration procedure of the full set of equations of motion, starting from an initial state vector \((x(t_0), y(t_0), z(t_0); u(t_0), v(t_0), w(t_0); \phi(t_0), \theta(t_0), \psi(t_0); p(t_0), q(t_0), r(t_0))\) at some initial time \(t_0\). First we advance the rotational quantities numerically by a small, finite, time step \(\Delta t\), then the translational quantities.
Advancing Rotational Quantities
- 1a) Compute Moments \(\mathbf{M}_B=(\mathcal{L}(t), \mathcal{M}(t), \mathcal{N}(t))\): in the body frame, using angle of attack \(\alpha\) and sideslip angle \(\beta\), as well as magnitude of the velocity vector, \(V = |V|\). The latter two angles can be obtained from the velocity vector in the body frame, \(\mathbf{V}_B = (u, v, w)\), with the equations \(\alpha=\arctan(u/w)\) and \(\beta=\arcsin(v/V)\). Gravity exerts no moment around the center of gravity (CG), so the moments do not depend on attitude with respect to the ground, only with respect to the relative wind. Thrust creates a moment if the thrust line does not pass through the CG of the aircraft (typically a pitch up or pitch down moment, if the engine(s) are mounted below or above the CG, respectively).
- 1b) Update Angular Velocity \(\underline{\omega}_B=(p, q, r)\): Compute \(\dot p(t), \dot q(t), \dot r(t)\), using the rotational dynamic equation, at time \(t\) from the current values of \(p(t), q(t), r(t)\) and \(\mathcal{L}(t), \mathcal{M}(t), \mathcal{N}(t)\). Advance body rates \(p(t), q(t), r(t)\) by time step \(\Delta t\):
\begin{eqnarray} p(t+\Delta t) &=& p(t) + \dot p(t)\Delta t\\ q(t+\Delta t) &=& q(t) + \dot q(t)\Delta t\\ r(t+\Delta t) &=& r(t) + \dot r(t)\Delta t
\end{eqnarray} - 2) Update Attitude (i.e. Euler Angles \((\phi, \theta, \psi)\)): Compute \(\dot \phi(t), \dot \theta(t), \dot \psi(t)\) from \(\dot p(t), \dot q(t), \dot r(t)\) using the differential rotational kinematic equation for the Euler angles (matrix \([B]\) is discussed in the article “Aircraft Attitude and Euler Angles”) – for this we will need \(\phi(t), \theta(t), \psi(t)\). Advance Euler angles \(\phi(t), \theta(t), \psi(t)\) by time step \(\Delta t\):
\begin{eqnarray}
\phi(t+\Delta t) &=& \phi(t) + \dot \phi(t)\Delta t \\
\theta(t+\Delta t) &=& \theta(t) + \dot \theta(t)\Delta t \\
\psi(t+\Delta t) &=& \psi(t) + \dot \psi(t)\Delta t
\end{eqnarray}
Advancing Translational Quantities
- 3a) Compute Forces \((F_{x_B}(t), F_{y_B}(t), F_{z_B}(t))\): in the body frame. Gravity will need to be brought in from the inertial frame, using the DCM \([BE]\). This will give a dependence on \(\phi(t)\) and \(\theta(t)\), because the direction of the weight force vector depends on attitude. Aerodynamic forces (lift \(L\), drag \(D\), sideforce \(Y\)) will have to be brought in from the wind frame by the DCM \([BW]\), as they depend on the angle of attack \(\alpha\) and sideslip angle \(\beta\).
- 3b) Update Translational Velocity \((u, v, w)\): Compute \(\dot u(t), \dot v(t), \dot w(t)\), using the translational dynamic equation at time \(t\) from the current values of \(u(t), v(t), w(t), p(t), q(t), r(t)\), and \(F_{x_B}(t), F_{y_B}(t), F_{z_B}(t)\). Advance velocity vector components \(u(t), v(t), w(t)\) (in the body frame) by time step \(\Delta t\):
\begin{eqnarray} u(t+\Delta t) &=& u(t)+\dot u(t)\Delta t\\ v(t+\Delta t) &=& v(t)+\dot v(t) \Delta t\\ w(t+\Delta t) &=& w(t)+\dot w(t)\Delta t
\end{eqnarray} - 4) Update Position \((x, y, z)\): in the inertial (world) frame. From the velocity components in the body frame, \(u(t), v(t), w(t)\), compute the velocity components in the world frame \(V_x(t), V_y(t), V_z(t)\) using the translational kinematic equation. This is a basis transformation accomplished by the DCM \([EB]\), discussed the article “Aircraft Attitude and Euler Angles”. For the DCM we will need \(\phi(t), \theta(t), \psi(t)\). Advance position vector components \(x(t), y(t), z(t)\) (in the world frame) by time step \(\Delta t\):
\begin{eqnarray} x(t+\Delta t) &=& x(t)+V_x(t)\Delta t \\
y(t+\Delta t) &=& y(t)+V_y(t)\Delta t \\
z(t+\Delta t) &=& z(t)+V_z(t)\Delta t
\end{eqnarray} - 5) Restart the whole process for the next time step \(\Delta t\), i.e. for time \(t+\Delta t\) (define the new t to be the old \(t+\Delta t\)).