Aircraft Stability Analysis

Aircraft Stability Analysis

Aircraft Stability Analysis: Axes, the Derivative Set, the Equations of Motion, the Modes and the Criteria

This page gives the model behind the tool, its conventions and what each result means. The model is the classical linear one: small disturbances about steady, straight flight of a rigid aeroplane, with the longitudinal motion (speed, angle of attack, pitch) separate from the lateral-directional motion (sideslip, roll, yaw). It is the model of the stability-derivative reports of the 1960s and 1970s and of every flight-dynamics textbook, and it is a study aid for preliminary design and teaching. It does not model stall, large angles, flexible structure, control-system dynamics, power effects beyond one thrust derivative, or Mach effects beyond the derivatives that are entered.

Axes and signs

Forces, moments, velocities and rates are in body-fixed stability axes: origin at the centre of gravity, x forward along the flight path of the trimmed condition, y out of the right wing, z down. The axes are fixed to the aeroplane once the trim is chosen. Rolling moment and roll rate p are positive right wing down, pitching moment and pitch rate q positive nose up, yawing moment and yaw rate r positive nose right. Sideslip β = v/V is positive with the wind on the right cheek. A control deflection is positive in the sense its derivatives are given in; the usual convention (elevator and aileron trailing edge down, rudder trailing edge left) gives a negative C_mδe and a negative C_nδr. Angles are radians inside the calculation; the form takes and shows degrees where it says so.

Positions in a derivative set (the moment reference point, the centre of gravity, the neutral point) are in the geometry axes of the lattice that produced them: x positive aft, z positive up.

The derivative set

The tool reads the "anowysis-derivative-set" (version 1), a JSON object written by the Vortex Lattice Analysis tool or typed by hand. Its blocks are reference (S, c, b and the moment reference point), condition (V, ρ, Mach, α, β, C_L, C_D, C_D0, e, altitude), mass (m, I_xx, I_yy, I_zz, I_xz, c.g. position), static (C_L0, C_m0, neutral point, static margin), derivatives and controls. Everything is in SI units and radians. The derivative block holds, for each of α, β and the rates p̂ = p b/(2V), q̂ = q c/(2V), r̂ = r b/(2V), the six coefficients C_L, C_D, C_Y, C_l, C_m, C_n per unit of the variable, plus C_L and C_m per unit α̇ c/(2V) and C_L, C_D, C_m per unit u/V. Each control has a kind (elevator, aileron, rudder, flap, other) and its six coefficients per radian. A value of null means "not provided".

What the tool does with a set, in order:

• A null speed, density, mass or inertia is taken from the flight-condition group of the form, and the page says which. Check that those form values belong to your aircraft: the form opens with the numbers of the Navion.
• If the moment reference point is not the centre of gravity, the moments are transferred. With (r_x, r_z) the forward and downward offsets of the c.g. from the reference point in stability axes, C_m gains r_z C_D/c − r_x C_L/c, C_l gains r_z C_Y/b and C_n loses r_x C_Y/b, column by column. The rate derivatives are first corrected for the change of axis of rotation: a pitch rate about the c.g. is the same pitch rate about the old point plus a plunge velocity, so every q derivative loses (2a/c) times the α derivative, every r derivative gains (2a/b) times the β derivative and every p derivative loses (2d/b) times the β derivative, where a and d are the forward and downward offsets of the old point from the c.g. The small speed term of that correction is neglected.
• Inertias are used as given (stability axes) or, by the "Inertia Axes" choice, rotated from body reference axes through the angle of attack: I_xx,s = I_xx cos²α + I_zz sin²α − I_xz sin 2α, I_zz,s = I_xx sin²α + I_zz cos²α + I_xz sin 2α, I_xz,s = (I_xx − I_zz) sin α cos α + I_xz cos 2α.
• Derivatives a vortex lattice cannot give are estimated when null, and the page marks them: C_Lu = M² C_L/(1 − M²) (Prandtl-Glauert), C_Du = 0, C_mu = 0; C_Lα̇ = 2 η V_H a_t dε/dα and C_mα̇ = −2 η V_H a_t (l_t/c) dε/dα from the tail data when they are entered (the lag of the wing downwash in reaching the tail), otherwise zero.

Typed derivatives and the two conventions

Typed derivatives are lines of name = value. Nondimensional names follow the set (CLa for C_Lα, Cmq, Clb, Cndr and so on). Dimensional names are those of the stability-derivative reports: X_u, X_w, Z_u, Z_w, Z_ẇ (Zwd), Z_q, M_u, M_w, M_ẇ (Mwd), M_q, the control terms Xde, Zde, Mde; Y_v, Y_p, Y_r, L_β (Lb), L_p, L_r, N_β, N_p, N_r, Lda, Nda, Ldr, Ndr, and the side-force control terms either as accelerations (Yda, Ydr) or divided by the speed as the reports print them (Ysda, Ysdr). The lines units = ft or units = m and primed = 1 or 0 say whether lengths are feet or metres and whether the lateral moments already contain the product of inertia. The conversion, with q̄ = ρV²/2 and k = q̄ S/(m V), is printed on the results page:

X_u = −k (2 C_D + C_Du) + X_Tu, Z_u = −k (2 C_L + C_Lu), X_w = k (C_L − C_Dα), Z_w = −k (C_Lα + C_D), Z_ẇ = −k (c/2V) C_Lα̇, Z_q = −k (c/2) C_Lq, M_u = q̄ S c C_mu/(I_yy V), M_w = q̄ S c C_mα/(I_yy V), M_ẇ = q̄ S c² C_mα̇/(2 I_yy V²), M_q = q̄ S c² C_mq/(2 I_yy V), Y_v = k C_Yβ, L_β = q̄ S b C_lβ/I_xx, L_p = q̄ S b² C_lp/(2 I_xx V), N_β = q̄ S b C_nβ/I_zz, N_r = q̄ S b² C_nr/(2 I_zz V), and the same pattern for the rest.

Dimensional derivatives given in body axes whose x axis lies at an angle α₀ above the flight path (the "Axis Angle" of the form; NASA CR-2144 uses the fuselage reference line) are rotated to stability axes before they are made nondimensional. The rotation moves a little of Z_ẇ, Z_q and M_ẇ into terms the set has no place for (X_ẇ, X_q, M_u̇); they are dropped and the page says so. Their effect is confined to the phugoid damping. With dimensional input the drag coefficient cannot be separated from its speed derivative: unless a line CD = value is given, the tool takes C_D = −X_u/(2k), which reproduces X_u exactly. The thrust option is not applied to dimensional data, which already contain their thrust terms.

For nondimensional data the thrust is taken equal to the drag (plus the weight component in a climb). Constant thrust adds nothing to X_u; constant power adds X_Tu = −T/(m V), and with a thrust line a distance z_T below the c.g. the moment derivative M_u gains X_Tu m z_T/I_yy.

Equations of motion

The equations are those of Appendix C of NASA CR-2144, written for stability axes (U₀ = V, W₀ = 0, θ₀ equal to the flight-path angle γ₀):

u̇ = X_u u + X_w w − g cos θ₀ θ + X_δ δ
(1 − Z_ẇ) ẇ = Z_u u + Z_w w + (U₀ + Z_q) q − g sin θ₀ θ + Z_δ δ
q̇ = M_u u + M_w w + M_ẇ ẇ + M_q q + M_δ δ
θ̇ = q

v̇ = Y_v v + Y_p p + (Y_r − U₀) r + g cos θ₀ φ + Y_δ δ
ṗ = L′_β v/V + L′_p p + L′_r r + L′_δ δ
ṙ = N′_β v/V + N′_p p + N′_r r + N′_δ δ
φ̇ = p + r tan θ₀

The primed derivatives absorb the product of inertia so that each rate equation holds one acceleration: L′ = (L + I_xz N/I_xx) G and N′ = (N + I_xz L/I_zz) G with G = 1/(1 − I_xz²/(I_xx I_zz)). The ẇ equation is solved for ẇ and substituted in the q̇ equation, which gives the two 4 × 4 state matrices printed on the page.

Eigenvalues, modes and eigenvectors

Each matrix has four eigenvalues, the roots of a quartic characteristic equation. They are found twice. The first method is the QR algorithm: the matrix is reduced to Hessenberg form by elementary similarity transformations and iterated with Francis's double shift, the method of the EISPACK routines. The second is the closed solution of the quartic: the characteristic polynomial is formed by the Faddeev-LeVerrier recursion, depressed, and split into two quadratics through the positive root of Ferrari's resolvent cubic. The page prints both and their difference, which is of the order of 10⁻¹³ for ordinary aircraft.

A complex pair σ ± jω_d is an oscillation with undamped natural frequency ω_n = √(σ² + ω_d²), damping ratio ζ = −σ/ω_n, period 2π/ω_d, time to half amplitude ln 2/|σ| (to double if σ is positive) and cycles to half amplitude equal to that time over the period. A real root λ is an exponential with time constant −1/λ. The longitudinal pairs are named by frequency: the faster is the short period (angle of attack and pitch rate at nearly constant speed), the slower the phugoid (an exchange of speed and height at nearly constant angle of attack). The lateral-directional roots are the Dutch roll pair, the large real roll subsidence root and the small real spiral root, which may be slightly unstable. When the pattern is different (two real short-period roots, a coupled roll-spiral oscillation) the page names the roots accordingly.

The eigenvector of a mode gives the relative size and phase of the states in it. The phasor tables show them with the velocities divided by V, so that u/V, α and β are comparable with the angles. The ratio |φ/β| of the Dutch roll comes from that table.

Classical approximations

• Short period: the (w, q) equations with the speed held constant, ω² = [Z_w M_q − (U₀ + Z_q) M_w]/(1 − Z_ẇ). Usually within a few per cent.
• Phugoid: Lanchester's result ω = √2 g/V and ζ = 1/(√2 L/D), and the two-degree form ω² = −g Z_u/V, 2ζω = −X_u. The frequency is usually within 10 to 15 %; the damping can be far off, because the exact value depends on small terms such as C_Dα, M_u and the thrust derivative.
• Roll: 1/T_R = −L′_p. Spiral: the ratio e/d of the last two coefficients of the quartic. Dutch roll: ω² = N′_β + Y_v N′_r and 2ζω = −(Y_v + N′_r), which ignores the rolling of the mode and can miss the damping badly for aircraft with strong dihedral effect.

Flying-quality Levels of MIL-F-8785C

The specification (5 November 1980) rates each mode as Level 1 (clearly adequate), Level 2 (adequate with increased pilot workload) or Level 3 (controllable with excessive workload) for an airplane Class (I small light, II medium, with II-L land based and II-C carrier based, III large heavy, IV high manoeuvrability) and a Flight Phase Category (A rapid manoeuvring and precise tracking, B gradual manoeuvres such as climb and cruise, C terminal phases). The tool uses these limits:

The acceleration sensitivity is n/α = q̄ S C_Lα/(m g), in g per radian, and ω_sp²/(n/α) is the control anticipation parameter. When ω_d²|φ/β| exceeds 20 (rad/s)², Table VI raises the required ζω by 0.014, 0.009 and 0.005 times the excess for Levels 1, 2 and 3. Class III airplanes need no more than ζ = 0.7. In Category C, Figure 3 also draws floors on ω_sp and n/α that carry no printed numbers; the tool uses 0.87 and 0.70 rad/s with n/α of 2.7 and 2.0 for Level 1, and 0.6 and 0.4 rad/s with 1.8 and 1.0 for Level 2 (Classes I, II-C, IV and Classes II-L, III), read from the graph. Two real short-period roots are rated by the damping ratio of their equivalent pair. The specification's requirements on equivalent systems, control forces, roll performance, residual oscillations and atmospheric disturbances are not checked, and the Class IV combat row of Table VI is not used.

Static stability, longitudinal

• Static margin −C_mα/C_Lα about the c.g.; positive is stable. The neutral point lies that fraction of the chord behind the c.g.
• Manoeuvre margin, stick fixed: static margin − ρ S c C_mq/(4m). The pitch damping of a pull-up adds to the stability, so the manoeuvre point lies behind the neutral point. Because C_mq itself changes with the c.g., the manoeuvre point is exact to first order in the shift.
• Trim needs the lift and moment at zero angle of attack and zero deflection (C_L0, C_m0 of the set). The two equations C_L = C_L0 + C_Lα α + C_Lδ δ and 0 = C_m0 + C_mα α + C_mδ δ give α and δ.
• Elevator per g: in a steady pull-up at constant speed the pitch rate is (n − 1) g/V. The increments of angle of attack and elevator for one g follow from the lift and moment equations with the q terms included; a negative (trailing edge up) elevator per g is the stable sense.
• Speed stability: dδ_e/dV = [−C_mα/(C_mδ C_Lα − C_mα C_Lδ)] (−2 C_L/V), compressibility left out; positive (more nose-down elevator to fly faster) is stable.
• c.g. limits: aft, the entered minimum static margin ahead of the neutral point; forward, the c.g. at which the entered elevator limit just trims the entered lift coefficient.

Static stability, lateral-directional

• Weathercock stability C_nβ positive; dihedral effect C_lβ negative.
• Spiral criterion C_lβ C_nr − C_nβ C_lr positive for a stable spiral in level flight. It is the sign of the constant term of the lateral quartic; the mode table gives the root itself.
• Steady sideslip (a crosswind landing with the wing down): for the entered β the aileron and rudder come from C_l = 0 and C_n = 0, and the bank from the side-force balance C_Y + C_L sin φ = 0. The crosswind component is V sin β.
• Steady roll: p b/(2V) = −C_lδa δ_a/C_lp, the single-degree-of-freedom value for a rigid wing, compared with an entered target (0.07 is the classical figure for transports and trainers and 0.09 for fighters; MIL-F-8785C itself specifies bank angle in a given time).
• Engine out: the entered yawing moment is balanced at zero sideslip by the rudder, with the aileron holding the roll and a small bank balancing the rudder side force.
• Adverse yaw: the sign of C_nδa/C_lδa.

Transfer-function numerators

The response of an output to a control is a ratio of polynomials whose denominator is the characteristic polynomial. The numerator is c adj(sI − A) b, which the Faddeev-LeVerrier recursion gives at once; its roots are the zeros. The page lists them for pitch attitude to elevator (1/T_θ1, 1/T_θ2), roll rate to aileron and bank angle to aileron (ζ_φ, ω_φ; the ratio ω_φ/ω_d shows how much the Dutch roll appears in a roll).

Time histories

The response to a control step, pulse or doublet, to an initial angle of attack, sideslip or bank, or to a gust is computed without numerical integration error: over each interval of constant input the state is advanced by the exponential of a matrix that contains the system and the input. The 1 − cos gust, w_g = (U/2)(1 − cos 2πt/T), is produced by two extra states that generate its cosine. A vertical gust is an upward air velocity and enters through the w column of the matrix (the aeroplane responds to w − w_g); a side gust is air moving to the right. The gust is applied to the whole aeroplane at once: its gradual penetration and unsteady lift are not modelled.

Root loci

The c.g. sweep moves the moment reference of the derivative set along x and solves the longitudinal matrix again at each static margin, with the inertias held. It shows the short-period frequency falling as the c.g. moves aft, a real root reaching the origin when the c.g. is at the neutral point, and the short-period stiffness vanishing near the manoeuvre point. The speed sweep changes V at constant weight and altitude: C_L follows 1/V², the drag follows the parabolic polar when the set gives C_D0 and e (otherwise C_D is held), and the nondimensional derivatives are held.

Rate dampers

A pitch damper feeds pitch rate to the elevator and a yaw damper feeds yaw rate to the rudder, each with a constant gain: the closed-loop matrix is A + B K. There are no actuator or sensor dynamics and no washout, so a yaw damper also opposes a steady turn and shifts the spiral root. The entries are signed so that a positive gain adds damping.

Built-in aircraft and what the reports print

Five flight conditions of three aircraft are built in, from two public NASA contractor reports that give derivatives and the factored transfer functions computed from them. The Jetstar in the power approach comes from the nondimensional Table VII-1 of NASA CR-2144, whose lateral derivatives and inertias are in body axes and are rotated to stability axes. The Jetstar cruise and the two Boeing 747 conditions come from the report's dimensional body-axis tables. The Navion comes from the nondimensional tables of NASA CR-96008; its Table X-A prints S = 180 ft², but its dimensional derivatives follow from the 184 ft² of its Figure X-1, which the tool uses, and its table omits C_mδe, which is computed from the printed M_δe. The results page compares the modes with the printed factors, both from the set and directly from the printed dimensional derivatives. CR-2144 prints three significant figures, so the short period, Dutch roll and roll root agree within 0.7 %, the spiral within 2 % and the phugoid damping, a small difference of rounded terms, within 10 %; the Navion agrees within 0.15 %. From the Jetstar's two-figure nondimensional table, which carries no speed or thrust-line derivative, the phugoid frequency comes out 3 % low and the roll-rate numerator damping 13 % low; the other factors agree within 0.7 %.

Limits of the method

• Linear, small disturbances; longitudinal and lateral motions uncoupled; rigid aeroplane; controls fixed.
• The derivatives are taken as constants at one flight condition. A vortex-lattice set is inviscid and attached: stall, propeller slipstream and compressibility beyond the Prandtl-Glauert factor are not in it.
• Estimated derivatives (α̇, speed) are first approximations; the phugoid damping in particular depends on them and on the thrust model.
• The flying-quality Levels compare modal numbers with the limits of a 1980 specification written for military airplanes; civil certification rules (14 CFR Parts 23 and 25) are stated differently.