Final Examples¶

Aircraft Modes from Stability Derivatives¶

Consider the following aircraft data for the Pipe Cherokee r PA-28-180 flying at 50m/s.


\(m\) = 1090kg	\(S\) = 15	\(I_{yy}\) = 1700
\(I_{xx}\) = 3100	\(I_{zz}\) = 1400	\(I_{xz}\) = 0
\(U_e\) = 50	\(\rho\)=1.06	\(\bar{c}\) = 1.6
\(C_{L_0}\) = 0.543	\(C_{D_0}\) = 0.0615	\(b\) = 9.11

\(X_u\) = -0.06728	\(Z_u\) = -0.396	\(M_u\) = 0.0
\(X_w\) = 0.02323	\(Z_w\) = -1.729	\(M_w\) = -0.2772
\(X_q\) = 0.0	\(Z_q\) = -1.6804	\(M_q\) = -2.207
\(M_{\dot{w}}\) = -0.0197	\(Z_{\delta_e}\) = -17.01	\(M_{\delta_e}\) = -44.71

\(Y_v\) = -0.1444	\(L_v\) = -0.1166	\(N_v\) = 0.174
	\(L_p\) = -2.283	\(N_p\) = -1.732
	\(L_r\) = 1.053	\(N_r\) = -1.029
\(Y_{\delta_r}\) = 2.113	\(L_{\delta_r}\) = 0.6133	\(N_{\delta_r}\) = -6.583
	\(L_{\delta_a}\) = 3.101	\(N_{\delta_a}\) = 0.0

The following example has been completed using Python, and I’m a better coder than when I first wrote the examples in MATLAB nearly three years ago, so my codes look more like programs than a list of calculator instructions.

I know that this can make things seem less accessible to some, but it’s smarter to keep data in a class like below.

# The means that I'm using the compute this is to store the Cherokee stability derivatives in a class,
# and then make a function to produce the stability deritatives. This is the easiest way to make readable code
# that's easily extensible to other aircraft without copying and pasting, or having all the derivatives live
# in the general namespace

import numpy as np
import sympy as sp
from IPython.display import display, Math, Latex, Markdown

# All data for Piper Cherokee PA-28-180
class Cherokee():
    h = 1200
    V = 50
    S = 15
    m = 1090
    Ixx = 1300 
    Izz = 1400
    Ixz = 0
    Iyy = 1700 
    CL0 = 0.543
    CD0 = 0.0615
    b = 9.11
    cbar = 1.3
    U0 = 50
    theta_0=0

    Ue = V

    g = 9.80665

    Xu = -0.06728
    Zu = -0.396
    Mu = 0.0
    Xw = 0.02323
    Zw = -1.729
    Mw = -0.2772
    Xq = 0.0
    Zq = -1.6804
    Mq = -2.207
    Mdw = -0.0197
    Zde = -17.01
    Mde = -44.71

    Yv = -0.1444
    Lv = -0.1166
    Nv = 0.174
    Lp = -2.283
    Np = -1.732
    Lr = 1.053
    Nr = -1.029
    Ydr = 2.113
    Ldr = 0.6133
    Ndr = -6.583
    Lda = 3.101
    Nda = 0.0
    
# Make A and B matrices for longitudinal motion
def MakeLongitudinal(ac):    
    # Determine the M* derivatives (assuming Theta_e = 0)
    Mu_star = ac.Mu + ac.Mdw*ac.Zu;
    Mw_star = ac.Mw + ac.Mdw*ac.Zw;
    Mq_star = ac.Mq + ac.Mdw*ac.Ue;
    Mth_star = 0;
    Mde_star = ac.Mde + ac.Mdw*ac.Zde;
    
    A = np.array([[ac.Xu, ac.Xw, 0, -ac.g],
                [ac.Zu, ac.Zw, ac.U0, 0],
                [Mu_star, Mw_star, Mq_star, Mth_star],
                [0, 0, 1, 0]])
    
    # Make the control/input matrix, for good measure
    B = np.array([[0], [ac.Zde], [Mde_star], [0]])
    return A, B

def MakeLateral(ac):
    # Making the starred terms
    Imess = ac.Ixx * ac.Izz / (ac.Ixx * ac.Izz - ac.Ixz**2)
    I2 = ac.Ixz / ac.Ixx

    # Lstarred terms
    Lvstar = Imess * (ac.Lv + I2 * ac.Nv)
    Lpstar = Imess * (ac.Lp + I2 * ac.Np)
    Lrstar = Imess * (ac.Lr + I2 * ac.Nr)
    Ldrstar = Imess * (ac.Ldr + I2 * ac.Ndr)
    Ldastar = Imess * (ac.Lda + I2 * ac.Nda)

    # Nstarred terms
    I2 = ac.Ixz / ac.Izz
    Nvstar = Imess * (ac.Nv + I2 * ac.Lv)
    Npstar = Imess * (ac.Np + I2 * ac.Lp)
    Nrstar = Imess * (ac.Nr + I2 * ac.Lr)
    Ndrstar = Imess * (ac.Ndr + I2 * ac.Ldr)
    Ndastar = Imess * (ac.Nda + I2 * ac.Lda)
    
    # Make the lateral matrix
    A = np.matrix([[ac.Yv, 0, -ac.U0, ac.g*np.cos(ac.theta_0), 0],
                   [Lvstar, Lpstar, Lrstar, 0, 0],
                   [Nvstar, Npstar, Nrstar, 0, 0],
                   [0, 1, np.tan(ac.theta_0), 0, 0],
                   [0, 0, 1/np.cos(ac.theta_0), 0, 0]])

    # And the control matrix
    if not hasattr(ac, 'Yda'):
        ac.Yda = 0
        
    B = np.matrix([[ac.Ydr, ac.Yda], [Ldrstar, Ldastar], [Ndrstar, Ndastar], [0, 0], [0, 0]])
    
    return A, B

# Make the two sets of matrices
Alon, Blon = MakeLongitudinal(Cherokee)
Alat, Blat = MakeLateral(Cherokee)

print("The longitudinal A and B matrices for this aircraft are")
display(Math('[A_{lon}] = ' + sp.latex(sp.Matrix(Alon).evalf(5))))
display(Math('[B_{lon}] = ' + sp.latex(sp.Matrix(Blon).evalf(5))))

print("The lateral/directional A and B matrices for this aircraft are")
display(Math('[A_{lat}] = ' + sp.latex(sp.Matrix(Alat).evalf(5))))
display(Math('[B_{lat}] = ' + sp.latex(sp.Matrix(Blat).evalf(5))))

The longitudinal A and B matrices for this aircraft are

\[\begin{split}\displaystyle [A_{lon}] = \left[\begin{matrix}-0.06728 & 0.02323 & 0 & -9.8067\\-0.396 & -1.729 & 50.0 & 0\\0.0078012 & -0.24314 & -3.192 & 0\\0 & 0 & 1.0 & 0\end{matrix}\right]\end{split}\]

\[\begin{split}\displaystyle [B_{lon}] = \left[\begin{matrix}0\\-17.01\\-44.375\\0\end{matrix}\right]\end{split}\]

The lateral/directional A and B matrices for this aircraft are

\[\begin{split}\displaystyle [A_{lat}] = \left[\begin{matrix}-0.1444 & 0 & -50.0 & 9.8067 & 0\\-0.1166 & -2.283 & 1.053 & 0 & 0\\0.174 & -1.732 & -1.029 & 0 & 0\\0 & 1.0 & 0 & 0 & 0\\0 & 0 & 1.0 & 0 & 0\end{matrix}\right]\end{split}\]

\[\begin{split}\displaystyle [B_{lat}] = \left[\begin{matrix}2.113 & 0\\0.6133 & 3.101\\-6.583 & 0\\0 & 0\\0 & 0\end{matrix}\right]\end{split}\]

# Get the eigenvalues
eigs_lon, eigv_lon = np.linalg.eig(Alon)

# We know these will always be a complex pair, so can separate. You could get the second one
# by fixing the index, but I like to do boolean indexing
eig1 = eigs_lon[0]
eig2 = eigs_lon[eigs_lon.real != eig1.real][0]

# A quick function to print the single eigenvalue as a complex pair
def makecomplexpair(eigin, roundto=4):
    return str(eigin.round(roundto)).replace('+', '±').strip('(').strip(')')

print(f"The eigenvalues associated with longtudinal motion are the complex pair {makecomplexpair(eig1)}, and {makecomplexpair(eig2)}")

The eigenvalues associated with longtudinal motion are the complex pair -2.4663±3.4056j, and -0.0279±0.2452j

# do the same for the lateral-directional modes
# Get the eigenvalues
eigs_lat, eigv_lat = np.linalg.eig(Alat)


# Find the dutch roll - the only part with a non-zero imaginary component
eig_DR = eigs_lat[eigs_lat.imag > 0][0]
other_lat_eigs = eigs_lat[eigs_lat.imag == 0]

print(f"The eigenvalues associated with lateral-direction motion are:\nThe complex pair {makecomplexpair(eig_DR)}")
# This next line is complicated only because I like to enjoy how wonderfully easy it is to manipulate strings in Python by contrast to MATLAB
print(f"and {len(other_lat_eigs)} other real eigenvalues:  {' '.join(str(other_lat_eigs.real.round(4)).split())[2:].replace(' ', ' and ').strip('[').strip(']')}")

The eigenvalues associated with lateral-direction motion are:
The complex pair -0.3468±3.3718j
and 3 other real eigenvalues:  0. and -2.7823 and 0.0194

Longitudinal Modes¶

The two eigenvalues are -2.4663±3.4056j and -0.0279±0.2452jj, and we can easily discriminate between the two modes because of the respective size. Recall that the damped natural frequency is given by the imaginary part of the eigenvalue, so we can see:

# A bit of logic to determine which mode is which
if eig1.imag > eig2.imag:
    SP = eig1
    PH = eig2
else:
    SP = eig2
    PH = eig1


    
print(f"From inspection of the imaginary part of the eigenvalues, {makecomplexpair(SP)} is clearly the short period mode, with a period of {2*np.pi/SP.imag:1.4}s\n")
print(f"From inspection of the imaginary part of the eigenvalues, {makecomplexpair(PH)} is clearly the Phugoid mode, with a period of {2*np.pi/PH.imag:1.4}s\n")

From inspection of the imaginary part of the eigenvalues, -2.4663±3.4056j is clearly the short period mode, with a period of 1.845s

From inspection of the imaginary part of the eigenvalues, -0.0279±0.2452j is clearly the Phugoid mode, with a period of 25.63s

Since both modes are stable, the time to half amplitude can be found from the real part of the eigenvalue

def timetohalf(eig):
    return -0.69/eig.real

for e, n in zip([SP, PH], ['Short Period', "Phugoid"]):
    print(f"The time to half amplitude of the {n} mode is {timetohalf(e):1.2f}s")

The time to half amplitude of the Short Period mode is 0.28s
The time to half amplitude of the Phugoid mode is 24.75s

The undamped natural frequency, \(\omega_n\), is given by the absolute value of the eigenvalue,

\[\omega_n=\sqrt{\Re\left({\lambda}\right)^2 + \Im\left({\lambda}\right)^2}\]

and hence the damping ratio, \(\zeta\) can be determine either from comparison with the standard form characteristic equation and

\(\zeta = -\frac{\Re(\lambda)}{\omega_n}\)

or from the definition of the damped natural frequency

\(\zeta = \sqrt{1-\frac{\omega_d^2}{\omega_n^2}}\)

for e, n in zip([SP, PH], ['Short Period', "Phugoid"]):
    print(f"The undamped natural frequency of the {n} mode is {np.abs(e):1.4f}rad/s")
    
    zeta = -e.real/np.abs(e)
    print(f"which gives a damping ratio, from method 1 of {zeta:1.4f}")
    
    zeta = np.sqrt(1-e.imag**2/np.abs(e)**2)
    print(f"and a damping ratio, from method 2 of {zeta:1.4f}\n")

The undamped natural frequency of the Short Period mode is 4.2048rad/s
which gives a damping ratio, from method 1 of 0.5865
and a damping ratio, from method 2 of 0.5865

The undamped natural frequency of the Phugoid mode is 0.2468rad/s
which gives a damping ratio, from method 1 of 0.1130
and a damping ratio, from method 2 of 0.1130

the two methodologies obviously are, and were always going to be equivalent, but it’s nice to show redundancy in methods.

Aircraft Flight Mechanics by Harry Smith, PhD

Final Examples¶

Aircraft Modes from Stability Derivatives¶

Longitudinal Modes¶