Final Examples

Final Examples#

Aircraft Modes from Stability Derivatives#

Consider the following aircraft data for the Pipe Cherokee 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.

Show code cell source Hide code cell source

# 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
import plotly.graph_objects as go


# 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;
    
    if not hasattr(ac, 'theta0'):
        ac.theta0 = 0
    
    A = np.array([[ac.Xu, ac.Xw, 0, -ac.g*np.cos(ac.theta0)],
                [ac.Zu, ac.Zw, ac.U0, -ac.g*np.sin(ac.theta0)],
                [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}\]

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

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.2452j, 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:

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

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}}\)

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.

Lateral Directional Modes#

From the A matrix for lateral directional motion the eigenvalues

-0.3468±3.3718j, 0.0,
-2.7823, 0.0194

are found. The first value is the only complex pair, and therefore has to be the Dutch Roll mode (the only oscillatory lateral/directional mode). We can see that is stable (negative real part) and the following can be yielded using the same as before:

From inspection of the imaginary part of the eigenvalue, -0.3468±3.3718j, the Dutch Roll mode, has a period of 1.863s
and a time to half amplitude of 1.99s.

For the remaining eigenvalues, one simply denotes the aircraft’s ability to yaw - the zero value shows that the aircraft has directional freedom and no real concept of heading stability (you could at most say the aircraft possesses neutral heading stability). The other two eigenvalues are real-only; -2.7823 and 0.0194.

We know that one denotes the spiral mode and one denotes the roll mode. The absolute magnitude of the damping indicates that -2.7823 is the roll mode, whilst it can also be seen that 0.0194 is unstable so cannot be the roll mode.

In summary:

For this aircraft the spiral mode is unstable with a time to double amplitude of 35.56s
For this aircraft the roll mode is stable (which is *has* to be) with a time to half amplitude of 0.25s

Extracting Aircraft Modes from Data#

Say you’ve been given some data for an aircraft subjected to a roll disturbance in a datafile - the information you have is roll attitude vs time in a text file.

The data is the free-response of the aircraft, just as seen in the examples following the impulsive forcing of the aircraft.

Plotting these data vs time is the first step, and we can see:

Note that the file itself doesn’t actually tell you the units of the roll attitude, but for determination of the modal qualities, it doesn’t acutally matter. The y axis values are only used to get the damping ratios.

You can determine that this file contains information about the lateral-directional modes since we’re looking at roll attitude. Furthermore, you can see:

Spiral mode and Dutch roll mode are superposed on the data.
The Dutch roll mode is stable (decreasing amplitude).
The Spiral mode is unstable (increasing amplitude).
The Roll mode is not visible.

This is the same as the example went over in class except the spiral mode is unstable. The theory is the same, though. We need to remove the effect of the spiral mode on the roll mode. You can do this by hand, or using computational methods.

The first step is to graphically separate the two modes. We will look at what the starting value of the spiral mode is. Visually, this is the intercept if we were to imagine the curve without the Dutch roll mode.

The spiral roll attitude is of the form

\[\phi_{sp}=\phi_0 e^{\lambda_{sp}t}\]

where \(\phi_{0}\) is intercept on the y-axis if the Dutch roll were not present.

We can see that if we extend the curve to the y-axis, it looks like it would intercept around 1.2 degrees, and has the two values given above. We can try this number and see if it works well by looking at the final value on the curve. That is, if we hvae \(\phi_0=1.2\), then use the fact that we know \(\phi(80s)=10\) and solve for \(\lambda_{sp}\)

Below shows two methods - one is guessing, which you can perform by trial and error and see which is a good fit through different data points, and the other is a means to fit the curve. That’s just there for advanced folks.

With a guess of 1.2 as the intercept, the eigenvalue for the spiral mode is 0.0265

Using the intial guess as the starting point for a least-squares fit gives 1.23as the intercept,
 the eigenvalue for the spiral mode is 0.02619

With the eigenvalue for the spiral mode now known, it can be removed from the data to look at the Dutch Roll on its own.

If you were doing this by hand, you could look at the peaks of the DR on its own, but it’s relatively easy to do this on the whole array. This will be completed with the least-squares fit and the guessed values from above.

# Using guessed values for the spiral mode intial value and eigenvalue, it can be subtracted from the Dutch Roll
just_dutch_roll = Phi_total - phi_0_guess * np.exp(t * lam_sp_guess)

# Repeats the original figure
fig.show()

# Adds on the bit we've just calculated
fig.add_trace(go.Scatter(x=t, y=just_dutch_roll.real, name="Spiral Removed"))

# Find some peak values using Scipy
from scipy.signal import find_peaks

# Get the indices of the peaks for the total signal and the bit we've just removed
ipeaks_total = find_peaks(Phi_total)[0]
ipeaks_dr = find_peaks(just_dutch_roll)[0]
# Add on the first couple of peaks to both
for i in range(2):
    fig.add_trace(go.Scatter(x=[t[ipeaks_total[i]]], y = [Phi_total[ipeaks_total[i]].real], mode="markers+text", text=f"{t[ipeaks_total[i]]:1.2f}, {Phi_total.real[ipeaks_total[i]]:1.2f}", textposition='top center', showlegend=False))
    fig.add_trace(go.Scatter(x=[t[ipeaks_dr[i]]], y = [just_dutch_roll[ipeaks_dr[i]].real], mode="markers+text", text=f"{t[ipeaks_dr[i]]:1.2f}, {just_dutch_roll.real[ipeaks_dr[i]]:1.2f}", textposition='bottom center', showlegend=False))

fig.show()

The damped natural frequency is available from the red plot. The logic used to display them on the graph isn’t necessary for you guys to do. You can do this by inspection of the graph - but it’s easier to do it as an example this way.

The damped natural frequency is:

# Get the damped natural frequency from the curve containing *just* the spiral mode
wd_dr = 2*np.pi / (t[ipeaks_dr[1] - ipeaks_dr[0]])

# If you didn't remove the spiral_mode
wd_dr_wrong = 2*np.pi / (t[ipeaks_total[1] - ipeaks_total[0]])

print(f"The damped natural frequency of the Dutch Roll mode is {wd_dr:1.4}rad/s, but you would have got {wd_dr_wrong:1.4f}rad/s if you didn't remove the spiral mode" )

The damped natural frequency of the Dutch Roll mode is 1.245rad/s, but you would have got 1.2260rad/s if you didn't remove the spiral mode

The damping ratio is given by the exponential decay:

# Remake the figure, again. Tis time we'll just include the DR
fig = go.Figure()
fig.update_xaxes(title_text="Time/s")
fig.update_yaxes(title_text="Roll Attitude/deg")
fig.add_trace(go.Scatter(x=t, y=just_dutch_roll.real, name="Spiral Removed"))


# Get the exponential decay thus giving the eigenvalue
lam_dr_real = np.log(just_dutch_roll.real[ipeaks_dr[1]]/just_dutch_roll.real[ipeaks_dr[0]])/t[ipeaks_dr[1] - ipeaks_dr[0]]
lam_dr_real_wrong = np.log(Phi_total.real[ipeaks_total[1]]/Phi_total.real[ipeaks_total[0]])/t[ipeaks_total[1] - ipeaks_total[0]]

print(f"The real part of the eigenvalue of the Dutch Roll mode is {lam_dr_real:1.3f}")

# Make the eigenvalue
eig_dr = lam_dr_real + wd_dr * 1j
eig_dr_wrong = lam_dr_real_wrong + wd_dr_wrong * 1j

print(f"This makes the total eigenvalue of the Dutch Roll mode is {makecomplexpair(eig_dr)}\n\n")

print(f"If you didn't get rid of the spiral mode, you would have got {makecomplexpair(eig_dr_wrong)}")

# Put on the graph to check
exp_bracket_upper = just_dutch_roll.real[0] * np.exp(eig_dr.real * t)
fig.add_trace(go.Scatter(x=t, y=exp_bracket_upper.real, name="Exponential Decay"))
fig.show()

The real part of the eigenvalue of the Dutch Roll mode is -0.144
This makes the total eigenvalue of the Dutch Roll mode is -0.1443±1.2454j


If you didn't get rid of the spiral mode, you would have got -0.0364±1.226j

The exponential decay ‘envelope’ is pretty good, but it neither tends to actual zero, which it should, not does it really go through the actual peaks. This is because the spiral mode hasn’t been fully removed from the data above due to the guess not being quite correct.

Now this will be repeated with the fitted values in place of the guessed ones for the Spiral mode: