{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "You guys have the supplementary notes in PDF form for the aircraft modes of motion, and the last bit is best taught by way of an example.\n", "\n", "# Final Examples\n", "\n", "## Aircraft Modes from Stability Derivatives\n", "\n", "Consider the following aircraft data for the Pipe Cherokee PA-28-180 flying at 50m/s.\n", "\n", "| | | |\n", "|:------------------------|:------------------------|:------------------------|\n", "| $m$ = 1090kg | $S$ = 15 | $I_{yy}$ = 1700 |\n", "| $I_{xx}$ = 3100 | $I_{zz}$ = 1400 | $I_{xz}$ = 0 |\n", "| $U_e$ = 50 | $\\rho$=1.06 | $\\bar{c}$ = 1.6 |\n", "| $C_{L_0}$ = 0.543 | $C_{D_0}$ = 0.0615 | $b$ = 9.11 |\n", "| | | |\n", "| $X_u$ = -0.06728 | $Z_u$ = -0.396 | $M_u$ = 0.0 |\n", "| $X_w$ = 0.02323 | $Z_w$ = -1.729 | $M_w$ = -0.2772 |\n", "| $X_q$ = 0.0 | $Z_q$ = -1.6804 | $M_q$ = -2.207 |\n", "| $M_{\\dot{w}}$ = -0.0197 | $Z_{\\delta_e}$ = -17.01 | $M_{\\delta_e}$ = -44.71 |\n", "| | | |\n", "| $Y_v$ = -0.1444 | $L_v$ = -0.1166 | $N_v$ = 0.174 |\n", "| | $L_p$ = -2.283 | $N_p$ = -1.732 |\n", "| | $L_r$ = 1.053 | $N_r$ = -1.029 |\n", "| $Y_{\\delta_r}$ = 2.113 | $L_{\\delta_r}$ = 0.6133 | $N_{\\delta_r}$ = -6.583 |\n", "| | $L_{\\delta_a}$ = 3.101 | $N_{\\delta_a}$ = 0.0 |\n", "\n", "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](https://aircraftflightmechanics.com/otherfiles/PiperCherokee.m), so my codes look more like programs than a list of calculator instructions.\n", "\n", "I know that this can make things seem _less_ accessible to some, but it's smarter to keep data in a class like below." ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "tags": [ "hide-input" ] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The longitudinal A and B matrices for this aircraft are\n" ] }, { "data": { "text/latex": [ "$\\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]$" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/latex": [ "$\\displaystyle [B_{lon}] = \\left[\\begin{matrix}0\\\\-17.01\\\\-44.375\\\\0\\end{matrix}\\right]$" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "The lateral/directional A and B matrices for this aircraft are\n" ] }, { "data": { "text/latex": [ "$\\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]$" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/latex": [ "$\\displaystyle [B_{lat}] = \\left[\\begin{matrix}2.113 & 0\\\\0.6133 & 3.101\\\\-6.583 & 0\\\\0 & 0\\\\0 & 0\\end{matrix}\\right]$" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# The means that I'm using the compute this is to store the Cherokee stability derivatives in a class,\n", "# and then make a function to produce the stability deritatives. This is the easiest way to make readable code\n", "# that's easily extensible to other aircraft without copying and pasting, or having all the derivatives live\n", "# in the general namespace\n", "\n", "import numpy as np\n", "import sympy as sp\n", "from IPython.display import display, Math, Latex, Markdown\n", "import plotly.graph_objects as go\n", "\n", "\n", "# All data for Piper Cherokee PA-28-180\n", "class Cherokee():\n", " h = 1200\n", " V = 50\n", " S = 15\n", " m = 1090\n", " Ixx = 1300 \n", " Izz = 1400\n", " Ixz = 0\n", " Iyy = 1700 \n", " CL0 = 0.543\n", " CD0 = 0.0615\n", " b = 9.11\n", " cbar = 1.3\n", " U0 = 50\n", " theta_0=0\n", "\n", " Ue = V\n", "\n", " g = 9.80665\n", "\n", " Xu = -0.06728\n", " Zu = -0.396\n", " Mu = 0.0\n", " Xw = 0.02323\n", " Zw = -1.729\n", " Mw = -0.2772\n", " Xq = 0.0\n", " Zq = -1.6804\n", " Mq = -2.207\n", " Mdw = -0.0197\n", " Zde = -17.01\n", " Mde = -44.71\n", "\n", " Yv = -0.1444\n", " Lv = -0.1166\n", " Nv = 0.174\n", " Lp = -2.283\n", " Np = -1.732\n", " Lr = 1.053\n", " Nr = -1.029\n", " Ydr = 2.113\n", " Ldr = 0.6133\n", " Ndr = -6.583\n", " Lda = 3.101\n", " Nda = 0.0\n", " \n", "# Make A and B matrices for longitudinal motion\n", "def MakeLongitudinal(ac): \n", " # Determine the M* derivatives (assuming Theta_e = 0)\n", " Mu_star = ac.Mu + ac.Mdw*ac.Zu;\n", " Mw_star = ac.Mw + ac.Mdw*ac.Zw;\n", " Mq_star = ac.Mq + ac.Mdw*ac.Ue;\n", " Mth_star = 0;\n", " Mde_star = ac.Mde + ac.Mdw*ac.Zde;\n", " \n", " if not hasattr(ac, 'theta0'):\n", " ac.theta0 = 0\n", " \n", " A = np.array([[ac.Xu, ac.Xw, 0, -ac.g*np.cos(ac.theta0)],\n", " [ac.Zu, ac.Zw, ac.U0, -ac.g*np.sin(ac.theta0)],\n", " [Mu_star, Mw_star, Mq_star, Mth_star],\n", " [0, 0, 1, 0]])\n", " \n", " # Make the control/input matrix, for good measure\n", " B = np.array([[0], [ac.Zde], [Mde_star], [0]])\n", " return A, B\n", "\n", "def MakeLateral(ac):\n", " # Making the starred terms\n", " Imess = ac.Ixx * ac.Izz / (ac.Ixx * ac.Izz - ac.Ixz**2)\n", " I2 = ac.Ixz / ac.Ixx\n", "\n", " # Lstarred terms\n", " Lvstar = Imess * (ac.Lv + I2 * ac.Nv)\n", " Lpstar = Imess * (ac.Lp + I2 * ac.Np)\n", " Lrstar = Imess * (ac.Lr + I2 * ac.Nr)\n", " Ldrstar = Imess * (ac.Ldr + I2 * ac.Ndr)\n", " Ldastar = Imess * (ac.Lda + I2 * ac.Nda)\n", "\n", " # Nstarred terms\n", " I2 = ac.Ixz / ac.Izz\n", " Nvstar = Imess * (ac.Nv + I2 * ac.Lv)\n", " Npstar = Imess * (ac.Np + I2 * ac.Lp)\n", " Nrstar = Imess * (ac.Nr + I2 * ac.Lr)\n", " Ndrstar = Imess * (ac.Ndr + I2 * ac.Ldr)\n", " Ndastar = Imess * (ac.Nda + I2 * ac.Lda)\n", " \n", " # Make the lateral matrix\n", " A = np.matrix([[ac.Yv, 0, -ac.U0, ac.g*np.cos(ac.theta_0), 0],\n", " [Lvstar, Lpstar, Lrstar, 0, 0],\n", " [Nvstar, Npstar, Nrstar, 0, 0],\n", " [0, 1, np.tan(ac.theta_0), 0, 0],\n", " [0, 0, 1/np.cos(ac.theta_0), 0, 0]])\n", "\n", " # And the control matrix\n", " if not hasattr(ac, 'Yda'):\n", " ac.Yda = 0\n", " \n", " B = np.matrix([[ac.Ydr, ac.Yda], [Ldrstar, Ldastar], [Ndrstar, Ndastar], [0, 0], [0, 0]])\n", " \n", " return A, B\n", "\n", "# Make the two sets of matrices\n", "Alon, Blon = MakeLongitudinal(Cherokee)\n", "Alat, Blat = MakeLateral(Cherokee)\n", "\n", "print(\"The longitudinal A and B matrices for this aircraft are\")\n", "display(Math('[A_{lon}] = ' + sp.latex(sp.Matrix(Alon).evalf(5))))\n", "display(Math('[B_{lon}] = ' + sp.latex(sp.Matrix(Blon).evalf(5))))\n", "\n", "print(\"The lateral/directional A and B matrices for this aircraft are\")\n", "display(Math('[A_{lat}] = ' + sp.latex(sp.Matrix(Alat).evalf(5))))\n", "display(Math('[B_{lat}] = ' + sp.latex(sp.Matrix(Blat).evalf(5))))" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "tags": [ "hide-input" ] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The eigenvalues associated with longtudinal motion are the complex pair -2.4663±3.4056j, and -0.0279±0.2452j\n" ] } ], "source": [ "# Get the eigenvalues\n", "eigs_lon, eigv_lon = np.linalg.eig(Alon)\n", "\n", "# We know these will always be a complex pair, so can separate. You could get the second one\n", "# by fixing the index, but I like to do boolean indexing\n", "eig1 = eigs_lon[0]\n", "eig2 = eigs_lon[eigs_lon.real != eig1.real][0]\n", "\n", "# A quick function to print the single eigenvalue as a complex pair\n", "def makecomplexpair(eigin, roundto=4):\n", " return str(eigin.round(roundto)).replace('+', '±').strip('(').strip(')')\n", "\n", "print(f\"The eigenvalues associated with longtudinal motion are the complex pair {makecomplexpair(eig1)}, and {makecomplexpair(eig2)}\")" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "tags": [ "remove-input" ] }, "outputs": [ { "data": { "application/papermill.record/text/plain": "'-2.4663±3.4056j'" }, "metadata": { "scrapbook": { "mime_prefix": "application/papermill.record/", "name": "eig1" } }, "output_type": "display_data" }, { "data": { "application/papermill.record/text/plain": "'-0.0279±0.2452j'" }, "metadata": { "scrapbook": { "mime_prefix": "application/papermill.record/", "name": "eig2" } }, "output_type": "display_data" } ], "source": [ "from myst_nb import glue\n", "glue(\"eig1\", f\"{makecomplexpair(eig1)}\", display=False)\n", "glue(\"eig2\", f\"{makecomplexpair(eig2)}\", display=False)\n", " \n" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "tags": [ "hide-input" ] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The eigenvalues associated with lateral-direction motion are:\n", "The complex pair -0.3468±3.3718j\n", "and 3 other real eigenvalues: 0. and -2.7823 and 0.0194\n" ] } ], "source": [ "# do the same for the lateral-directional modes\n", "# Get the eigenvalues\n", "eigs_lat, eigv_lat = np.linalg.eig(Alat)\n", "\n", "\n", "# Find the dutch roll - the only part with a non-zero imaginary component\n", "eig_DR = eigs_lat[eigs_lat.imag > 0][0]\n", "other_lat_eigs = eigs_lat[eigs_lat.imag == 0]\n", "\n", "print(f\"The eigenvalues associated with lateral-direction motion are:\\nThe complex pair {makecomplexpair(eig_DR)}\")\n", "# This next line is complicated but elegant because I like to enjoy how wonderfully easy it is to manipulate strings in Python by contrast to MATLAB\n", "print(f\"and {len(other_lat_eigs)} other real eigenvalues: {' '.join(str(other_lat_eigs.real.round(4)).split())[2:].replace(' ', ' and ').strip('[').strip(']')}\")\n" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "tags": [ "remove-input" ] }, "outputs": [ { "data": { "application/papermill.record/text/plain": "'-0.3468±3.3718j'" }, "metadata": { "scrapbook": { "mime_prefix": "application/papermill.record/", "name": "eig3" } }, "output_type": "display_data" }, { "data": { "application/papermill.record/text/plain": "'0.0'" }, "metadata": { "scrapbook": { "mime_prefix": "application/papermill.record/", "name": "eig4" } }, "output_type": "display_data" }, { "data": { "application/papermill.record/text/plain": "'-2.7823'" }, "metadata": { "scrapbook": { "mime_prefix": "application/papermill.record/", "name": "eig5" } }, "output_type": "display_data" }, { "data": { "application/papermill.record/text/plain": "'0.0194'" }, "metadata": { "scrapbook": { "mime_prefix": "application/papermill.record/", "name": "eig6" } }, "output_type": "display_data" } ], "source": [ "from myst_nb import glue\n", "glue(\"eig3\", f\"{makecomplexpair(eig_DR)}\", display=False)\n", "glue(\"eig4\", f\"{str(other_lat_eigs[0].real.round(4))}\", display=False)\n", "glue(\"eig5\", f\"{str(other_lat_eigs[1].real.round(4))}\", display=False)\n", "glue(\"eig6\", f\"{str(other_lat_eigs[2].real.round(4))}\", display=False)\n", " \n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Longitudinal Modes\n", "\n", "The two eigenvalues are {glue:text}`eig1` and {glue:text}`eig2`, 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:" ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "tags": [ "hide-input" ] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "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\n", "\n", "From inspection of the imaginary part of the eigenvalues, -0.0279±0.2452j is clearly the Phugoid mode, with a period of 25.63s\n", "\n" ] } ], "source": [ "# A bit of logic to determine which mode is which\n", "if eig1.imag > eig2.imag:\n", " SP = eig1\n", " PH = eig2\n", "else:\n", " SP = eig2\n", " PH = eig1\n", "\n", "\n", " \n", "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\")\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\")\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Since both modes are stable, the time to _half_ amplitude can be found from the real part of the eigenvalue" ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "tags": [ "hide-input" ] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The time to half amplitude of the Short Period mode is 0.28s\n", "The time to half amplitude of the Phugoid mode is 24.75s\n" ] } ], "source": [ "def timetohalf(eig):\n", " return -0.69/eig.real\n", "\n", "for e, n in zip([SP, PH], ['Short Period', \"Phugoid\"]):\n", " print(f\"The time to half amplitude of the {n} mode is {timetohalf(e):1.2f}s\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The undamped natural frequency, $\\omega_n$, is given by the absolute value of the eigenvalue,\n", "\n", "$$\\omega_n=\\sqrt{\\Re\\left({\\lambda}\\right)^2 + \\Im\\left({\\lambda}\\right)^2}$$\n", "\n", "and hence the damping ratio, $\\zeta$ can be determine either from comparison with the standard form characteristic equation and\n", "\n", "$\\zeta = -\\frac{\\Re(\\lambda)}{\\omega_n}$\n", "\n", "or from the definition of the damped natural frequency\n", "\n", "$\\zeta = \\sqrt{1-\\frac{\\omega_d^2}{\\omega_n^2}}$" ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "tags": [ "hide-input" ] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The undamped natural frequency of the Short Period mode is 4.2048rad/s\n", "which gives a damping ratio, from method 1 of 0.5865\n", "and a damping ratio, from method 2 of 0.5865\n", "\n", "The undamped natural frequency of the Phugoid mode is 0.2468rad/s\n", "which gives a damping ratio, from method 1 of 0.1130\n", "and a damping ratio, from method 2 of 0.1130\n", "\n" ] } ], "source": [ "for e, n in zip([SP, PH], ['Short Period', \"Phugoid\"]):\n", " print(f\"The undamped natural frequency of the {n} mode is {np.abs(e):1.4f}rad/s\")\n", " \n", " zeta = -e.real/np.abs(e)\n", " print(f\"which gives a damping ratio, from method 1 of {zeta:1.4f}\")\n", " \n", " zeta = np.sqrt(1-e.imag**2/np.abs(e)**2)\n", " print(f\"and a damping ratio, from method 2 of {zeta:1.4f}\\n\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "the two methodologies obviously are, and were always going to be equivalent, but it's nice to show redundancy in methods." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Lateral Directional Modes\n", "\n", "From the A matrix for lateral directional motion the eigenvalues \n", "\n", "{glue:text}`eig3`, \n", "{glue:text}`eig4`, \n", "{glue:text}`eig5`, \n", "{glue:text}`eig6` \n", "\n", "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:" ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "tags": [ "hide-input" ] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "From inspection of the imaginary part of the eigenvalue, -0.3468±3.3718j, the Dutch Roll mode, has a period of 1.863s\n", "and a time to half amplitude of 1.99s.\n" ] } ], "source": [ "print(f\"From inspection of the imaginary part of the eigenvalue, {makecomplexpair(eig_DR)}, the Dutch Roll mode, has a period of {2*np.pi/eig_DR.imag:1.4}s\\n\\\n", "and a time to half amplitude of {-0.69/eig_DR.real:1.4}s.\")\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "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; {glue:text}`eig5` and {glue:text}`eig6`. \n", "\n", "We know that one denotes the **spiral** mode and one denotes the **roll mode**. The absolute magnitude of the damping indicates that {glue:text}`eig5` is the **roll mode**, whilst it can also be seen that {glue:text}`eig6` is _unstable_ so _cannot be the roll mode_.\n", "\n", "In summary:" ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "tags": [ "hide-input" ] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "For this aircraft the spiral mode is unstable with a time to double amplitude of 35.56s\n", "For this aircraft the roll mode is stable (which is *has* to be) with a time to half amplitude of 0.25s\n" ] } ], "source": [ "\n", "print(f\"For this aircraft the spiral mode is unstable with a time to double amplitude of {0.69/other_lat_eigs[2].real:1.2f}s\")\n", "\n", "print(f\"For this aircraft the roll mode is stable (which is *has* to be) with a time to half amplitude of {0-0.69/other_lat_eigs[1].real:1.2f}s\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Extracting Aircraft Modes from Data\n", "\n", "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](https://www.aircraftflightmechanics.com/Data/RollDisturbance.csv).\n", "\n", "The data is the free-response of the aircraft, just as seen in the examples following the _impulsive_ forcing of the aircraft." ] }, { "cell_type": "code", "execution_count": 12, "metadata": { "tags": [ "remove-input" ] }, "outputs": [], "source": [ "\n", "# Make up some data - first, spiral mode\n", "Phi_dist = 1.23\n", "Phi_final = 10\n", "t_final = 80\n", "\n", "# get the eigenvalue\n", "lam_spiral_input = np.log(Phi_final/Phi_dist)/t_final\n", "\n", "\n", "# Simulate the data\n", "t = np.linspace(0, t_final, 1000)\n", "Phi_spiral = Phi_dist * np.exp(lam_spiral_input * t)\n", "\n", "# # Just plot the roll attitude\n", "# fig = go.Figure()\n", "# fig.update_layout(title=\"Roll Attitude\", title_x=0.5)\n", "# fig.add_trace(go.Scatter(x=t, y=Phi_spiral.real))\n", "\n", "# Add in a Dutch Roll Mode\n", "a = -0.15 # Real part\n", "w = 1.24 # Frequency\n", "Phi_dist_dr = 2.57\n", "\n", "Phi_dr = Phi_dist_dr*np.exp(a*t + w*1j*(t+np.random.random(1)[0]*0))\n", "\n", "Phi_total = Phi_spiral + Phi_dr\n", "\n", "\n", "\n", "# # Make a data file of these \n", "data = np.array([t.transpose(), Phi_total.real.transpose()]).transpose()\n", "\n", "np.savetxt('RollDisturbance.csv', data, delimiter=',', header=\"Time, Phi\");\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Plotting these data vs time is the first step, and we can see:" ] }, { "cell_type": "code", "execution_count": 13, "metadata": { "tags": [ "remove-input" ] }, "outputs": [ { "data": { "text/html": [ "