{"id":174,"date":"2017-02-18T20:08:06","date_gmt":"2017-02-18T20:08:06","guid":{"rendered":"http:\/\/www.science-now.org\/giampaolo\/?p=174"},"modified":"2018-08-25T12:37:28","modified_gmt":"2018-08-25T11:37:28","slug":"flight-dynamics-interactively","status":"publish","type":"post","link":"https:\/\/www.science-now.org\/giampaolo\/2017\/02\/18\/flight-dynamics-interactively\/","title":{"rendered":"Flight dynamics, interactively!"},"content":{"rendered":"

It’s possible to match two interesting things: flight dynamics and interactive programming with Python. The magic tool we’re going to use is called Jupyter Notebook. Basically, it’s Python programming done with our favourite web browser.
\n<\/p>\n

As we seen in the previous post about flight dynamics, the coupling of rigid body motions and aerodynamics leads to a peculiar aircraft behaviour. Using a Jupyter notebook, we’re going to integrate the equations of motion as function of time. We’ll produce some plots (like the one below!) showing the short and long term response of the airplane.<\/p>\n

\"\"<\/p>\n

The rigid rotation is shown in this plot. The short term oscillations disappear after few seconds since the short period mode is strongly damped. However, the long term response is dominated by the phugoid mode which has low frequency and little damping.<\/p>\n

I’ve included the Jupyter notebook I used to produce this plot here below. The implementation is quite technical but I tried to include some comments to make it easier to read. However, if you are interested, please feel free to contact me!<\/p>\n

Fligth dynamics, interactively!<\/h1>\n

Let’s import some Python modules. We’re going to use some modules from standard library (time, math) but also NumPy and SciPy<\/p>\n

\r\nimport time\r\nimport math\r\n\r\n# Numpy and Scipy\r\nimport scipy as sc\r\nimport numpy as np\r\n\r\nimport scipy.integrate as si\r\n\r\nimport scipy.linalg as sl\r\nimport numpy.linalg as nl\r\n\r\nfrom numpy import pi as pi\r\nfrom scipy.optimize import fsolve as fsolve\r\n\r\n# Matplotlib\r\n%matplotlib inline\r\nimport matplotlib.pyplot as plt\r\nimport matplotlib.pylab as pylab\r\n\r\npylab.rcParams['figure.figsize'] = (8, 5)\r\npylab.rcParams['figure.dpi'] = 96\r\npylab.rcParams['savefig.dpi'] = 96\r\n<\/pre>\n
Parameters<\/h1>\n
The test case we’re going to implement is taken from literature, specifically from the book “Flight Dynamics Principles” – Cook, M.V., 2nd edition at page 106.<\/p>\n
\r\n# Units: Imperial\r\n\r\n#-------------------------------- General parameters\r\n#- Parameters at sea level: air density and thrust\r\nTSL = 1000\r\nrho0 = 1.10\r\n\r\n#- Initial conditions: speed, angle of attack,\r\n#  Euler angle and air density at given altitude\r\nV0 = 305\r\nalpha0 = 0.\r\ntheta0 = 0.\r\nzE0  = 0.\r\nrho = 0.00238\r\n\r\n#- Mass and inertia properties\r\nmass = 746\r\nIy = 65000\r\ngravity = 32.2\r\n\r\n#-------------------------------- Stability derivatives\r\n# These derivatives represent how the aerodynamic forces\r\n# change due to the flight dynamics unknowns. For example,\r\n# Zw represents the change in the vertical force due to the\r\n# vertical velocity w\r\n\r\n#- Derivatives - Z direction\r\nZu = -159.64\r\nZu_dot = 0.\r\n\r\nZw = -328.24\r\nZw_dot = 0\r\n\r\nZq = 0\r\nZq_dot = 0\r\n\r\nZdelta = -16502.0\r\n\r\n#- Derivatives - X direction\r\nXu = -26.26\r\nXu_dot = 0.\r\n\r\nXw = 79.82\r\nXw_dot = 0\r\n\r\nXq = 0\r\nXq_dot = 0\r\n\r\nXdelta = 0\r\n\r\n#- Derivatives - Rotation\r\nMu = 0\r\nMu_dot = 0\r\n\r\nMw = -1014.0\r\nMw_dot = -36.4\r\n\r\nMq = -18135.0\r\nMq_dot = 0\r\n\r\nMdelta = -303575.0\r\n\r\n#-------------------------------- Simulation options\r\n# We are going to simulate the first 60 seconds. The solution to\r\n# the differential equations is calculated each 0.1 seconds\r\n\r\nsimTime = [0, 60]\r\nsimStep = 0.1\r\n\r\n# Horizontal and vertical velocity at equilibrium condition\r\nUe = V0*np.cos(alpha0)\r\nWe = V0*np.sin(alpha0)\r\n\r\n# Initial conditions. We're going to simulate an initial disturbance\r\n# in the angular velocity q\r\nstate_vector0 = [0, 0, 0.1, 0]\r\nelevator0 = 0\r\nthrottle0 = 0\r\n<\/pre>\n
Equations of motion<\/h1>\n
These differential equations express the evolution of the system with respect to time. This is the set of equations that we’ll integrate to get the flight dynamics response of the system.<\/p>\n
$$\\begin{align}
\n\\left[ \\begin{array}{cccc}
\nm-X_{\\dot u} & -X_{\\dot w} & -X_{\\dot q} & 0\\\\
\n-Z_{\\dot u} & m-Z_{\\dot w} & -Z_{\\dot q} & 0 \\\\
\n-M_{\\dot u} & -M_{\\dot w} & I_y-M_{\\dot q} & 0\\\\
\n0 & 0 & 0 & 1
\n\\end{array} \\right]
\n\\left[ \\begin{array}{c}
\n{\\dot u} \\\\
\n{\\dot w} \\\\
\n{\\dot q} \\\\
\n{\\dot \\theta} \\\\
\n\\end{array} \\right] = &
\n\\left[ \\begin{array}{cccc}
\nX_{u} & X_{w} & X_{q} & 0 \\\\
\nZ_{u} & Z_{w} & Z_{q} & 0 \\\\
\nM_{u} & M_{w} &M_{q} & 0 \\\\
\n0 & 0 & 0 & 0
\n\\end{array} \\right]
\n\\left[ \\begin{array}{c}
\n{ u } \\\\
\n{ w } \\\\
\n{ q } \\\\
\n{\\theta} \\\\
\n\\end{array} \\right] + \\\\ &
\n\\left[ \\begin{array}{cccc}
\n0 & 0 & – m W_e & – m g \\cos \\theta_0 \\\\
\n0 & 0 & m U_e & – m g \\sin \\theta_0 \\\\
\n0 & 0 & 0 & 0 \\\\
\n0 & 0 & 1 & 0
\n\\end{array} \\right]
\n\\left[ \\begin{array}{c}
\n{ u} \\\\
\n{ w} \\\\
\n{ q} \\\\
\n{\\theta} \\\\
\n\\end{array} \\right] + \\\\ &
\n\\left[ \\begin{array}{c}
\nX_\\delta \\\\
\nZ_\\delta \\\\
\nM_\\delta \\\\
\n0
\n\\end{array} \\right]
\n\\delta +
\n\\left[ \\begin{array}{c}
\nT_{SL} \\frac{\\rho}{\\rho_0} \\\\
\n0 \\\\
\n0 \\\\
\n0
\n\\end{array} \\right] \\eta
\n\\end{align}
\n$$<\/p>\n
Matrices and eigenvalues<\/h2>\n
\r\ndef getMatrices(Ue, We):\r\n    # This function return the Jacobian matrix of the system\r\n    # and the other vectors related to control surfaces and thrust\r\n    # that we need for the integration\r\n\r\n    # ---------\r\n    I1 = np.array([[mass-Xu_dot, -Xw_dot, -Xq_dot, 0.],\r\n                   [-Zu_dot, mass-Zw_dot, -Zq_dot, 0.],\r\n                   [-Mu_dot, -Mw_dot, Iy-Mq_dot, 0.],\r\n                   [      0,       0,        0, 1]])\r\n\r\n    In = nl.inv(I1)\r\n\r\n\r\n    # ---------\r\n    An = np.array([[Xu, Xw, Xq, 0],\r\n                   [Zu, Zw, Zq, 0],\r\n                   [Mu, Mw, Mq, 0],\r\n                   [0,0,0,0]])\r\n\r\n    Kn = np.array([[0, 0, - mass * We, - mass * gravity * np.cos(theta0)],\r\n                   [0, 0, + mass * Ue, - mass * gravity * np.sin(theta0)],\r\n                   [0, 0, 0, 0],\r\n                   [0, 0, 1, 0]])\r\n\r\n    A = An + Kn\r\n\r\n    # ---------\r\n    D = np.array([Xdelta, Zdelta, Mdelta, 0])\r\n    T = np.array([TSL * rho \/ rho0, 0, 0, 0])\r\n\r\n    return np.dot(In, A), np.dot(In, D), np.dot(In, T)\r\n<\/pre>\n
Using the function we have written, we can calculate the eigenvalues of the system. We’re expecting two eigenvalues corresponding to short period and phugoid!<\/p>\n
\r\nJac, Jac_D, Jac_T = getMatrices(V0, 0.)\r\n\r\nprint "Jacobian matrix = "\r\nprint Jac\r\nprint\r\n\r\nprint "Eigenvalue analysis: "\r\nevals = nl.eig(Jac)[0]\r\nfor ev in evals:\r\n    print "\\t", ev\r\nprint\r\n<\/pre>\n
\r\nJacobian matrix =\r\n[[ -3.52010724e-02   1.06997319e-01   0.00000000e+00  -3.22000000e+01]\r\n [ -2.13994638e-01  -4.40000000e-01   3.05000000e+02   0.00000000e+00]\r\n [  1.19836997e-04  -1.53536000e-02  -4.49800000e-01   0.00000000e+00]\r\n [  0.00000000e+00   0.00000000e+00   1.00000000e+00   0.00000000e+00]]\r\n\r\nEigenvalue analysis:\r\n\t(-0.44586984133+2.16437192886j)\r\n\t(-0.44586984133-2.16437192886j)\r\n\t(-0.0166306948626+0.147431081456j)\r\n\t(-0.0166306948626-0.147431081456j)\r\n<\/pre>\n
Eigenvalues are complex conjugate since they represent harmonic motions. The imaginary part of the eigenvalue is related to oscillation frequency. The eigenvalue with higher frequency is the short period whereas the low-frequency mode is the phugoid.<\/p>\n
Literature comparison<\/h3>\n
The solution to the eigenvalue problem can be derived analitically. The final result is copied here from literature as a comparison with the evaluation we made. In oarticular, transfer function is given as
\n$$
\nTF(s) = \\frac{-4.664 (s+0.135) (s+0.267)}{(s^2+0.033 s + 0.022) (s^2 + 0.893 s + 4.884)}
\n$$<\/p>\n
with stability modes given by the solution of the denominator
\n$$
\ns_1 = -0.0165 \\pm 0.1474 \\mathrm{i}
\n$$<\/p>\n
and<\/p>\n
$$
\ns_2 = -0.4465 \\pm 2.164 \\mathrm{i}
\n$$<\/p>\n
As you can see, the numeric evaluation gives very similar results.<\/p>\n
Integration<\/h2>\n
At this point, it is time to integrate the equations and have some nice plots of the flight dynamics response.<\/p>\n
\r\ndef integrate(y0, tSpan, dt, update=None, callback=None, maxInnerSteps=10000):\r\n    '''\r\n    Integrate a set of first order differential equation. This function acts\r\n    similarly to MATLAB's ode45 function.\r\n\r\n    t, y, timing = integrate(y0, tSpan, dt, update, callback, maxInnerSteps)\r\n\r\n    INPUT\r\n      y0 = initial conditions (y at tBeg)\r\n      tSpan = [tBeg, tEnd]\r\n      dt = step size\r\n\r\n      update = fn(t, y)\r\n      callback = fn(t, y) called at every succeful integration\r\n\r\n      maxInnerSteps = max number of inner iterations (default: 10000)\r\n\r\n    OUTPUT\r\n      t = list of timesteps\r\n      y = solutions at timestep (one line per timestep)\r\n      timing = elapsed time for the integration of the timestep\r\n    '''\r\n\r\n    tBegin = min(tSpan)\r\n    tEnd = max(tSpan)\r\n\r\n    history_t = []\r\n    history_y = []\r\n\r\n    integr = sc.integrate.ode(update)\r\n    integr.set_integrator('dopri5', nsteps=maxInnerSteps)\r\n    integr.set_initial_value(y0, tBegin)\r\n\r\n    t_elapsed = []\r\n    tic = time.time()\r\n\r\n    history_t.append(tBegin)\r\n    history_y.append(np.array(y0))\r\n\r\n    while integr.t &lt; tEnd:\r\n        values = integr.integrate(integr.t + dt)\r\n\r\n        if not integr.successful():\r\n            raise Exception("Integration not successful")\r\n\r\n        history_t.append(integr.t)\r\n        history_y.append(values)\r\n\r\n        if callback is not None:\r\n            try:\r\n                callback(integr.t, values)\r\n            except:\r\n                print "Error in the callback function"\r\n\r\n        toc = time.time() - tic\r\n        t_elapsed.append(toc)\r\n\r\n    return np.array(history_t), np.array(history_y), np.array(t_elapsed)\r\n\r\ndef eqsOfMotion(time, state_vector, elevator, throttle, Ue, We):\r\n    """\r\n    This function return the first order derivative of the state vector at each timestep.\r\n    The derivative is calculated using the Jacobian matrix evaluated earlier.\r\n    """\r\n\r\n    if np.any(np.isnan(state_vector)):\r\n        raise ValueError("NaN!")\r\n\r\n    Jac, Jac_D, Jac_T = getMatrices(Ue, We)\r\n\r\n    state_vector_dot = np.dot(Jac, state_vector) \\\r\n        + Jac_D * elevator \\\r\n        + Jac_T * throttle\r\n\r\n    return state_vector_dot\r\n<\/pre>\n
It’s good to output the initial conditions and the initial external disturbance<\/p>\n
\r\nprint "Initial state vector:"\r\nprint state_vector0\r\nprint\r\n\r\nprint "Elevator = ", elevator0 * (180.\/np.pi)\r\nprint "Throttle = ", throttle0\r\n\r\nfn_integrate = lambda t, x: eqsOfMotion(t, x, elevator0, throttle0, Ue, We)\r\nt_time, t_values, t_elapsed = integrate(np.array(state_vector0), simTime, simStep, fn_integrate, None)\r\n\r\nprint "Integration performed!"\r\n<\/pre>\n
Results<\/h1>\n
In this section the integration is performed and four plots are produced as result of the integration. Short period and phugoid are clearly visible in the response of the system.<\/p>\n
\r\nnames = ["u", "w", "q", "theta"]\r\ntitles = ["Horizontal speed", "Vertical speed", "Angular speed", "Rotation"]\r\nunits = ["[ft\/s]", "[ft\/s]", "[deg\/s]", "[deg]"]\r\nconvertion = [1, 1, (180.\/np.pi), (180.\/np.pi)]\r\n\r\nplt.figure()\r\n\r\nfor index in range(4):\r\n\r\n    plt.figure(index +1)\r\n    plt.plot(t_time, t_values[:,index]*convertion[index])\r\n    plt.ylabel("{} {}".format(names[index], units[index]))\r\n    plt.xlabel("Time [s]")\r\n    plt.title(titles[index])\r\n    plt.grid('on')\r\n\r\nplt.show()\r\n<\/pre>\n
<\/p>\n
<\/p>\n
<\/p>\n
<\/p>\n
Calculating flight path<\/h2>\n
The flight path can be calculated as post-processing step.<\/p>\n
\r\n# The derivatives of the horizontal and vertical position in the inertial frame are calculated\r\nxE_dot = +(Ue + t_values[:,0]) * np.cos(t_values[:,3]) + (We + t_values[:,1]) * np.sin(t_values[:,3])\r\nzE_dot = +(Ue + t_values[:,0]) * np.sin(t_values[:,3]) - (We + t_values[:,1]) * np.cos(t_values[:,3])\r\n\r\n# When these derivatives are needed at time steps not sampled, we have to use interpolation.\r\n# This is needed because the integration routine can choose time steps for which no evaluation of\r\n# xE_dot or zE_dot is available\r\ndef int_path(t, x):\r\n    xE_dot_temp = np.interp(t, t_time, xE_dot)\r\n    zE_dot_temp = np.interp(t, t_time, zE_dot)\r\n    return np.array([xE_dot_temp, zE_dot_temp])\r\n\r\npath_time, path_values, _ = integrate(np.array([0,0]), simTime, simStep, int_path, None)\r\n\r\nplt.figure()\r\nplt.plot(path_values[:,0], path_values[:,1])\r\nplt.xlabel('Horizontal translation [ft]')\r\nplt.ylabel('Vertical translation [ft]')\r\nplt.title('Flight path')\r\nplt.grid()\r\nplt.show()\r\n<\/pre>\n
<\/p>\n","protected":false},"excerpt":{"rendered":"
It’s possible to match two interesting things: flight dynamics and interactive programming with Python. The magic tool we’re going to use is called Jupyter Notebook. Basically, it’s Python programming done with our favourite web browser.<\/p>\n","protected":false},"author":1,"featured_media":0,"comment_status":"open","ping_status":"open","sticky":false,"template":"","format":"standard","meta":{"footnotes":""},"categories":[21,9],"tags":[22,26,10],"class_list":["post-174","post","type-post","status-publish","format-standard","hentry","category-aeronautical-engineering","category-programming","tag-flight-dynamics","tag-jupyter","tag-python"],"_links":{"self":[{"href":"https:\/\/www.science-now.org\/giampaolo\/wp-json\/wp\/v2\/posts\/174","targetHints":{"allow":["GET"]}}],"collection":[{"href":"https:\/\/www.science-now.org\/giampaolo\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/www.science-now.org\/giampaolo\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/www.science-now.org\/giampaolo\/wp-json\/wp\/v2\/users\/1"}],"replies":[{"embeddable":true,"href":"https:\/\/www.science-now.org\/giampaolo\/wp-json\/wp\/v2\/comments?post=174"}],"version-history":[{"count":7,"href":"https:\/\/www.science-now.org\/giampaolo\/wp-json\/wp\/v2\/posts\/174\/revisions"}],"predecessor-version":[{"id":224,"href":"https:\/\/www.science-now.org\/giampaolo\/wp-json\/wp\/v2\/posts\/174\/revisions\/224"}],"wp:attachment":[{"href":"https:\/\/www.science-now.org\/giampaolo\/wp-json\/wp\/v2\/media?parent=174"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/www.science-now.org\/giampaolo\/wp-json\/wp\/v2\/categories?post=174"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/www.science-now.org\/giampaolo\/wp-json\/wp\/v2\/tags?post=174"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}