{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Worked Example\n", "\n", "Data for several aircraft are presented in {cite}Heffley:1972wb and, honestly, this is a good exercise in showing you what a tremendous amount of being an engineer is - _extracting data from old reports and doing stuff with it_.\n", "\n", "## Data for Boeing 747\n", "\n", "Data are given on Page 217 of the aforementioned report for the B-747 in power approach configuration.\n", "\n", "### Configuration Information\n", "\n", "$$h=SL$$\n", "$$V_\\infty=165\\text{KTAS}$$\n", "$$S=5500\\text{ft}^2$$\n", "$$b=195.68\\text{ft}$$\n", "$$\\bar{c}=27.31\\text{ft}$$\n", "$$\\theta_0=0$$\n", "\n", "_note that the trim theta is not actually given, but assumed to be zero_\n", "\n", "### Nondimensional Stability Derivatives\n", "\n", "| Longitudinal \t| Lateral-Directional \t|\n", "|------------------------------------------\t|---------------------------------------------\t|\n", "| $C_L=1.11$ \t| $C_{y_\\beta}=-0.96\\text{rad}^{-1}$ \t|\n", "| $C_D=0.102$ \t| $C_{\\ell_\\beta}=-0.221\\text{rad}^{-1}$ \t|\n", "| $C_{L_\\alpha}=5.70\\text{rad}^{-1}$ \t| $C_{n_\\beta}=0.150\\text{rad}^{-1}$ \t|\n", "| $C_{D_\\alpha}=0.66\\text{rad}^{-1}$ \t| $C_{\\ell_\\hat{p}}=-0.45$ \t|\n", "| $C_{m_\\alpha}=-1.26\\text{rad}^{-1}$ \t| $C_{n_\\hat{p}}=-0.121$ \t|\n", "| $C_{L_\\dot{\\alpha}}=-6.7\\text{rad}^{-1}$ \t| $C_{\\ell_\\hat{r}}=0.101$ \t|\n", "| $C_{m_\\dot{\\alpha}}=-3.2\\text{rad}^{-1}$ \t| $C_{n_\\hat{r}}=-0.30$ \t|\n", "| $C_{L_\\hat{q}}=5.40$ \t| $C_{\\ell_{\\delta_a}}=0.0461\\text{rad}^{-1}$ \t|\n", "| $C_{m_\\hat{q}}=-20.8$ \t| $C_{n_{\\delta_a}}=0.0064\\text{rad}^{-1}$ \t|\n", "| $C_{L_\\text{M}}=-0.81$ \t| $C_{y_{\\delta_r}}=0.175\\text{rad}^{-1}$ \t|\n", "| $C_{m_\\text{M}}=0.27$ \t| $C_{\\ell_{\\delta_r}}=0.007\\text{rad}^{-1}$ \t|\n", "| $C_{L_{\\delta_e}}=0.338\\text{rad}^{-1}$ \t| $C_{n_{\\delta_r}}=-0.109\\text{rad}^{-1}$ \t|\n", "| $C_{m_{\\delta_e}}=-1.34\\text{rad}^{-1}$ \t| \t|\n", "\n", "Note that in {cite}Heffley:1972wb, the derivatives with respect to angular rates are erroneously given units of /rad, when they clearly are non-dimensional rates. The derivatives with respect to aerodynamic angles correctly are given units of /rad." ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "scrolled": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Xu = -0.0212\n", "Xw = 0.0467\n", "Zu = -0.2098\n", "Zw = -0.6027\n", "Zdw = -0.0341\n", "Zq = -7.6596\n", "Mu = 0.0001\n", "Mw = -0.0019\n", "Mw = -0.0019\n", "Mdw = -0.0002\n", "Mq = -0.4373\n", "Zde = -9.7779\n", "Mde = -0.5746\n" ] } ], "source": [ "import numpy as np\n", "from IPython.display import display, Math, Latex, Markdown\n", "\n", "\n", "# Condition and geometry information\n", "rho = 0.002377 # slugs/cubic feet\n", "# rho = rho * 32.174 # to lb / cubic feet\n", "rho = 0.0765 # lb/ft^3\n", "U0 = 165 # Knots\n", "U0 = U0 * 1.68781 # to feet/second\n", "S = 5500 # square feet\n", "b = 195.68 # feet\n", "c = 27.31 # feet\n", "theta_0 = 0\n", "m = 564000 # lb\n", "\n", "\n", "Iyy = 32.3e6 # slug-ft^2\n", "Ixx = 14.3e6\n", "Izz = 45.3e6\n", "Ixz = -2.23e6\n", "\n", "SIunits = False\n", "\n", "if SIunits:\n", " # Conversion factors\n", " ft_to_metre = 0.3048\n", " lb_to_kg = 0.453592\n", " slug_to_kg = 14.5939\n", " \n", " MOIconvert = slug_to_kg*ft_to_metre**2\n", "\n", " \n", " # Constants\n", " a = 340 # sonic velocity in m/s\n", " g = 9.80665 # acceleration due to gravity\n", " rho = 1.225 # density in kg/m^3\n", " \n", " # Convert\n", " U0 = U0 * ft_to_metre\n", " S = S * ft_to_metre**2\n", " b = b * ft_to_metre\n", " c = c * ft_to_metre\n", " m = m * lb_to_kg\n", " Ixx = Ixx * MOIconvert\n", " Iyy = Iyy * MOIconvert\n", " Izz = Izz * MOIconvert\n", " Ixz = Ixz * MOIconvert \n", " \n", "else:\n", " a = 1125.33 # sonic velocity in ft/s\n", " g = 32.174\n", " \n", " # Convert mass moments of inertia to consistent units --> INTO lb-ft^2 FROM slug-ft^2\n", " Ixx = Ixx * g\n", " Iyy = Iyy * g\n", " Izz = Izz * g\n", " Ixz = Ixz * g\n", "\n", "# Store the derivatives in two dictionaries\n", "B747_lon_ders = {'C_L': 1.11, 'C_D': 0.102, 'C_L_a' : 5.7, 'C_D_a' : 0.66, 'C_m_a' : -1.26,\n", " 'C_L_da' : -6.7, 'C_m_da' : -3.2, 'C_L_hq' : 5.4, 'C_m_hq' : -20.8, 'C_L_M' : -0.81,\n", " 'C_m_M' : 0.27, 'C_L_de' : 0.338, 'C_m_de' : -1.34}\n", "\n", "B747_lat_ders = {'C_y_b' : -0.96, 'C_l_b' : -0.221, 'C_n_b' : 0.150, 'C_l_hp' : -0.45, 'C_n_hp' : -0.121,\n", " 'C_l_hr' : 0.101, 'C_n_hr': -0.30, 'C_l_da' : 0.0461, 'C_n_da' : 0.0064, 'C_y_dr' : 0.175,\n", " 'C_l_dr' : 0.007, 'C_n_dr' : -0.109}\n", "\n", "# Put the dictionaries into the local namespace\n", "# This might seem a bit convoluted but it enables us to store values of derivatives in the dicionary, above,\n", "# and then put them all into the global namespace.\n", "#\n", "# It's fairly easy to get (key: value) pairs into a dictionary from a text file or xls so this will be handy\n", "# later.\n", "locals().update(B747_lon_ders)\n", "\n", "# Note that without doing the dictionary --> local namespace, the we'd have to write:\n", "# Xu = q * S / m / U0 * (2 * B747_lon_ders[\"C_D\"] + M * B747_lon_ders[\"C_D_M\"])\n", "\n", "# Convert to dimensional form\n", "q = 0.5 * rho * U0**2\n", "M = U0 / a\n", "\n", "Xu = -q * S / m / U0 * (2 * C_D) # No C_D_M term so assumed zero\n", "Xw = q * S / m / U0 * (C_L - C_D_a)\n", "Xq = 0 # No CDq term given in the table so assumed zero\n", "Zu = -q * S / m / U0 * (2 * C_L + M * C_L_M)\n", "Zw = -q * S / m / U0 * (C_D + C_L_a)\n", "Zdw = q * S * c / m / 2 / U0**2 * C_L_da # This is a NEW term for us, but since it was given as C_L_da, must be included\n", "Zq = -q * S * c / 2 / m / U0 * C_L_hq\n", "Mu = q * S * c / Iyy / U0 * M * C_m_M\n", "Mw = q * S * c / Iyy / U0 * C_m_a\n", "Mdw = q * S * c**2 / 2 / Iyy / U0**2 * C_m_da\n", "Mq = q * S * c**2 / 2 / Iyy / U0 * C_m_hq\n", "Mq = q * S * c**2 / 2 / Iyy / U0 * C_m_hq\n", "Zde = -q * S / m * C_L_de\n", "Mde = q * S * c / Iyy * C_m_de\n", "\n", "print(f\"Xu = {Xu:1.4f}\")\n", "print(f\"Xw = {Xw:1.4f}\")\n", "print(f\"Zu = {Zu:1.4f}\")\n", "print(f\"Zw = {Zw:1.4f}\")\n", "print(f\"Zdw = {Zdw:1.4f}\")\n", "print(f\"Zq = {Zq:1.4f}\")\n", "print(f\"Mu = {Mu:1.4f}\")\n", "print(f\"Mw = {Mw:1.4f}\")\n", "print(f\"Mw = {Mw:1.4f}\")\n", "print(f\"Mdw = {Mdw:1.4f}\")\n", "print(f\"Mq = {Mq:1.4f}\")\n", "\n", "print(f\"Zde = {Zde:1.4f}\")\n", "print(f\"Mde = {Mde:1.4f}\")\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Some units have dimensions of $[\\text{T}^{-1}]$ and hence have the same value in either unit system. Expand the below to see this. Other derivatives do not." ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "tags": [ "hide-input" ] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Using US Customary Units:\n", "Dynamic pressure, q = 2966.5143 lbs/ft^2\n", "Wing area, S = 5500 ft^2\n", "Mass, m = 564000 lbs\n", "Trim forward speed, U0 = 278.5 ft/s\n", "Gives a value for Xu = -0.2098, and a value for Zq of -7.6596\n", "\n", "Using SI Units:\n", "Dynamic pressure, q = 4414.9750 Pa\n", "Wing area, S = 511 m^2\n", "Mass, m = 255826 kg\n", "Trim forward speed, U0 = 84.9 m/s\n", "Gives a value for Xu = -0.2098, and a value for Zq of -2.3348\n", "\n" ] } ], "source": [ "# Demonstration of the equality of some derivatives in both unit systems:\n", "SIunits = {\"Pressure\": \"Pa\", \"Mass\": \"kg\", \"Area\" : \"m^2\", \"Speed\" : \"m/s\"}\n", "USunits = {\"Pressure\": \"lbs/ft^2\", \"Mass\": \"lbs\", \"Area\" : \"ft^2\", \"Speed\" : \"ft/s\"}\n", "\n", "# Conversion factors\n", "ft_to_metre = 0.3048\n", "lb_to_kg = 0.453592\n", "slug_to_kg = 14.5939\n", "\n", "for SIUNIT in [False, True]:\n", " if SIUNIT:\n", " UNITS = SIunits\n", " else:\n", " UNITS = USunits\n", " \n", " # Convert to some temporary variables to not overwrite the original values\n", " if SIUNIT:\n", " u0 = U0 * ft_to_metre\n", " s = S * ft_to_metre**2\n", " c2 = c * ft_to_metre\n", " mass = m * lb_to_kg\n", " rho = 1.225\n", " qinf = 0.5 * 1.2255 * u0**2\n", " else:\n", " qinf = q\n", " u0 = U0\n", " s = S\n", " c2 = c\n", " mass = m\n", "\n", " if SIUNIT: print(\"Using SI Units:\") \n", " else: print(\"Using US Customary Units:\")\n", " print(f\"Dynamic pressure, q = {qinf:1.4f} {UNITS['Pressure']}\")\n", " print(f\"Wing area, S = {s:1.0f} {UNITS['Area']}\")\n", " print(f\"Mass, m = {mass:1.0f} {UNITS['Mass']}\")\n", " print(f\"Trim forward speed, U0 = {u0:1.1f} {UNITS['Speed']}\")\n", " \n", " testXu = -qinf * s / mass / u0 * (C_L * 2 + C_L_M * M) \n", " testZq = -qinf * s * c2 / 2 / mass / u0 * C_L_hq\n", "\n", " print(f\"Gives a value for Xu = {testXu:1.4f}, and a value for Zq of {testZq:1.4f}\\n\")\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now the system matrix can be constructed:\n", "\n", "\\begin{aligned}\n", " \\begin{bmatrix} \\dot{u}\\\\\\dot{w}\\\\\\dot{q}\\\\\\dot{\\theta}\\end{bmatrix} &= \\begin{bmatrix}\n", " X_u & X_w & 0 & -g\\cdot\\cos\\Theta_e\\\\\n", " Z_u & Z_w & U_0 & -g\\cdot\\sin\\Theta_e\\\\\n", " M_u^* & M_w^* & M_q^* & M_\\theta^*\\\\\n", " 0 & 0 & 1 & 0 \n", " \\end{bmatrix}\\begin{bmatrix}\n", " {u}\\\\{w}\\\\{q}\\\\{\\theta} \n", " \\end{bmatrix} + \\begin{bmatrix}\n", " 0\\\\Z_{\\delta_e}\\\\M_{\\delta_e}^*\\\\0 \n", " \\end{bmatrix}\\left[\\delta_e\\right]\\end{aligned}\n", " \n", "well...almost - you'll see two things that need to be performed:\n", "1. The starred terms need to be made from the products of the stability derivatives.\n", "2. We have a $Z_q$ term from conversion of the stability derivatives, but no corresponding place to put it.\n", "\n", "The first problem is easy. From Equations {eq}eq:conciseqterms" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "scrolled": true }, "outputs": [], "source": [ "Mustar = Mu + Mdw * Zu\n", "Mwstar = Mw + Mdw * Zw\n", "Mqstar = Mq + Mdw * Zq\n", "Mthetastar = -Mdw * g * np.sin(theta_0)\n", "Mdestar = Mde + Mdw * Zde" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we've got all the terms, we don't know where to put $Z_q$. The answer should be obvious - it goes where the $U_0$ term is - it doesn't _replace_ it, but is _added to it_. So:" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "scrolled": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The system matrix for longitudinal motion is \n" ] }, { "data": { "text/latex": [ "$\\displaystyle [A_{lon}] = \\left[\\begin{matrix}-0.021191 & 0.046745 & 0 & -32.174\\\\-0.20979 & -0.6027 & 270.83 & 0\\\\0.00015355 & -0.0017943 & -0.43544 & 0\\\\0 & 0 & 1.0 & 0\\end{matrix}\\right]$" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "The eigenvalues are: [-0.52829455+0.69158467j -0.52829455-0.69158467j -0.00137077+0.14113918j\n", " -0.00137077-0.14113918j]\n", "The coefficients of the CE are: [1. 1.05933064 0.7802033 0.02312592 0.01508872]\n" ] } ], "source": [ "import sympy as sp\n", "\n", "Alon = np.matrix([[Xu, Xw, 0, -g*np.cos(theta_0)],\n", " [Zu, Zw, U0 + Zq, -g*np.sin(theta_0)],\n", " [Mustar, Mwstar, Mqstar, Mthetastar],\n", " [0, 0, 1, 0]])\n", "\n", "print(\"The system matrix for longitudinal motion is \")\n", "display(Math('[A_{lon}] = ' + sp.latex(sp.Matrix(Alon).evalf(5))))\n", "\n", "# Get the eigenvalues - this is my own checking I got the matrix correct. It'll make sense later.\n", "eigs, _ = np.linalg.eig(Alon)\n", "\n", "print(\"The eigenvalues are:\", eigs)\n", "\n", "print(\"The coefficients of the CE are:\", np.poly(Alon))\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The same can be performed for the lateral-directional matrix:" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "scrolled": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Yv = -0.0997\n", "Yp = 0.0000\n", "Yr = 0.0000\n", "Lv = -0.0055\n", "Lp = -1.0971\n", "Lr = 0.2462\n", "Nv = 0.0012\n", "Np = -0.0931\n", "Nr = -0.2309\n" ] } ], "source": [ "# Put the non-dimensional derivatives into the local namespace\n", "locals().update(B747_lat_ders)\n", "\n", "Yv = q * S / m / U0 * C_y_b\n", "Yp = 0\n", "Yr = 0\n", "Lv = q * S * b / Ixx / U0 * C_l_b\n", "Lp = q * S * b**2 / 2/ Ixx / U0 * C_l_hp\n", "Lr = q * S * b**2 / 2/ Ixx / U0 * C_l_hr\n", "Nv = q * S * b / Izz / U0 * C_n_b\n", "Np = q * S * b**2 / 2 / Izz / U0 * C_n_hp\n", "Nr = q * S * b**2 / 2 / Izz / U0 * C_n_hr\n", "\n", "\n", "print(f\"Yv = {Yv:1.4f}\")\n", "print(f\"Yp = {Yp:1.4f}\")\n", "print(f\"Yr = {Yr:1.4f}\")\n", "print(f\"Lv = {Lv:1.4f}\")\n", "print(f\"Lp = {Lp:1.4f}\")\n", "print(f\"Lr = {Lr:1.4f}\")\n", "print(f\"Nv = {Nv:1.4f}\")\n", "print(f\"Np = {Np:1.4f}\")\n", "print(f\"Nr = {Nr:1.4f}\")\n", "\n", "# Do the sadme with the control terms \n", "# 'C_l_da' : 0.0461, 'C_n_da' : 0.0064, 'C_y_dr' : 0.175,\n", "# 'C_l_dr' : 0.007, 'C_n_dr' : -0.109}\n", "Yda = 0\n", "Ydr = q * S / m * C_y_dr\n", "Lda = q * S * b / Ixx * C_l_da\n", "Ldr = q * S * b / Ixx * C_l_dr\n", "Nda = q * S * b / Izz * C_n_da\n", "Ndr = q * S * b / Izz * C_n_dr\n", "\n", "\n", "# \n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And put into the system matrix for lateral-directional motion\n", "\n", "$$\\boldsymbol{A}=\\begin{bmatrix}Y_v & 0 & -U_0 & g\\cdot\\cos\\theta_0 & 0\\\\L_v^* & L_p^* & L_r^* & 0 & 0\\\\N_v^* & N_p^* & N_r^* & 0 & 0\\\\0 & 1 & \\tan\\theta_0 & 0 & 0 \\\\ 0 & 0 & \\sec\\theta_0 & 0 & 0\\end{bmatrix}$$" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "scrolled": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The system matrix for lateral-directional motion is \n" ] }, { "data": { "text/latex": [ "$\\displaystyle [A_{lat}] = \\left[\\begin{matrix}-0.099723 & 0 & -278.49 & 32.174 & 0\\\\-0.0057348 & -1.0909 & 0.28442 & 0 & 0\\\\0.0014622 & -0.039417 & -0.24488 & 0 & 0\\\\0 & 1.0 & 0 & 0 & 0\\\\0 & 0 & 1.0 & 0 & 0\\end{matrix}\\right]$" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "The eigenvalues are: [ 0. +0.j -1.2285872 +0.j -0.08023692+0.74204291j\n", " -0.08023692-0.74204291j -0.04646794+0.j ]\n", "The coefficients of the CE are: [1. 1.43552897 0.81876855 0.71945085 0.03180283 0. ]\n" ] } ], "source": [ "# Making the starred terms\n", "Imess = Ixx * Izz / (Ixx * Izz - Ixz**2)\n", "I2 = Ixz / Ixx\n", "\n", "# Lstarred terms\n", "Lvstar = Imess * (Lv + I2 * Nv)\n", "Lpstar = Imess * (Lp + I2 * Np)\n", "Lrstar = Imess * (Lr + I2 * Nr)\n", "Ldrstar = Imess * (Ldr + I2 * Ndr)\n", "Ldastar = Imess * (Lda + I2 * Nda)\n", "\n", "# Nstarred terms\n", "I2 = Ixz / Izz\n", "Nvstar = Imess * (Nv + I2 * Lv)\n", "Npstar = Imess * (Np + I2 * Lp)\n", "Nrstar = Imess * (Nr + I2 * Lr)\n", "Ndrstar = Imess * (Ndr + I2 * Ldr)\n", "Ndastar = Imess * (Nda + I2 * Lda)\n", "\n", "Alat = np.matrix([[Yv, 0, -U0, g*np.cos(theta_0), 0],\n", " [Lvstar, Lpstar, Lrstar, 0, 0],\n", " [Nvstar, Npstar, Nrstar, 0, 0],\n", " [0, 1, np.tan(theta_0), 0, 0],\n", " [0, 0, 1/np.cos(theta_0), 0, 0]])\n", "\n", "print(\"The system matrix for lateral-directional motion is \")\n", "display(Math('[A_{lat}] = ' + sp.latex(sp.Matrix(Alat).evalf(5))))\n", "\n", "print(\"The eigenvalues are:\", np.linalg.eig(Alat)[0])\n", "\n", "print(\"The coefficients of the CE are:\", np.poly(Alat))" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "8.467904726657125" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.pi*2/.7420" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Transient Longitudinal Response\n", "\n", "Note that we've not made either of the control matrices yet - and there's a good reason for that. If we consider \n", "the linearised model of the aircraft created, the _stability_ of the aircraft (_i.e.,_ its response to a perturbation) will be governed by the system matrix.\n", "\n", "In the first instance, it is desired to understand the _stick fixed stability_ and understand what will happen to the aircraft if no inputs are made.\n", "\n", "Consider what the two matrices developed actually show us; if the aircraft is disturbed from a given condition, then if the disturbance is a _symmetric_ state variable ($u, w, q, \\theta$), then the output is a time rate of change of _all_ other symmetric variables. \n", "\n", "The control matrices tell us how the aircraft responds due to control input. Therefore, if we're interested in aircraft response in the absense of pilot or control-system input, then we're only interested in the A matrices themselves. The A matrices give the _open loop response_ of the aircraft.\n", "\n", "Since these are ODEs, we could write a means to time-march through the solution - but there are some handy inbuilt tools in Python to look at the responses of state-space systems to different inputs.\n", "\n", "You will probably need to install these to work on your system - [look here](https://python-control.readthedocs.io/en/0.8.3/intro.html#overview-of-the-toolbox).\n", "\n", "But if we create the $B$ matrix, it can be used to excite the aircraft with unity impulse input. For this reason, it makes sense to make the two control matrices with derivatives that have been scaled by a factor of $\\frac{\\pi}{180}$ to allow a unity _degree_ control input." ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "text/html": [ "