draft

Dynamic Stability

In this course, the equations of motion have been developed in a Newtonian framework - these, in conjunction with the relationship between body rates and the time rate of change of the Euler angles, allow us to describe full unconstrained flight of a six degree of freedom aircraft.

Regarding static stability, some basic static stability properties of fixed-wing aircraft were explored - what the immediate tendency of the aircraft is following a disturbance from trim. In the previous module we developed the equations of motion into two sets of uncoupled linearised equations of motion, representing the longitudinal and lateral/directional motions of fixed wing aircraft.

In this final module, we look at the dynamic stability characteristics of an aircraft - how the aircraft actually responds with time, and how the response to basic pilot inputs may be determined. The stability of any dynamic system is governed by consideration of its free (unforced) motion. Its response to an input (an aircraft control) is determined by the response of the system to forcing.

A simple example that helps us understand dynamic systems is the mass/spring/damper system as shown:

Mass Spring Damper System
Fig. 1:

Mass Spring Damper System

This system is governed by the equation of motion:

mx¨+cx˙+kx=F(t)m\ddot{x} + c\dot{x} + kx = F(t)

where xx is the displacement, mm is the mass, cc is the damping, kk is the spring stiffness, and F(t)F(t) is the input forcing. For free motion, F(t)=0F(t)=0:

mx¨+cx˙+kx=0m\ddot{x} + c\dot{x} + kx = 0

which we can represent as:

x¨+2ζωnx˙+ωn2=0(1)\ddot{x} + 2\zeta\omega_n\dot{x} + \omega_n^2 = 0\tag{1}

where ωn\omega_n is the undamped natural frequency of the system (=km=\sqrt{\frac{k}{m}} for this case) and ζ\zeta is the damping ratio of the system (=c2km=\frac{c}{2\sqrt{km}} for this case).

The stability of this system is given by solution of Eq (1) directly. The general solution is found by assuming that the solution will have the form x=Aeλtx=Ae^{\lambda_t}, hence x˙=Aλeλt,x¨=Aλ2eλt\dot{x}=A\lambda e^{\lambda t}, \ddot{x} = A\lambda^2 e^{\lambda t}. Substitution of this into Equation (1) yields the characteristic equation:

λ2+2ζωnλ+ωn2=0(2)\lambda^2 + 2\zeta\omega_n\lambda + \omega_n^2 = 0\tag{2}

the roots of which are

λ=ωn(ζ±ζ21)(3)\lambda = \omega_n\left(-\zeta \pm \sqrt{-\zeta^2-1}\right)\tag{3}

by superposition, the general solution is

x(t)=A1eλ1t+A2eλ2tx(t)=A1eωn(ζ+ζ21)t+A2eωn(ζζ21)t(4)\begin{aligned} x(t) &= A_1\cdot e^{\lambda_1t} + A_2\cdot e^{\lambda_2t}\\ x(t) &= A_1\cdot e^{\omega_n\left(-\zeta + \sqrt{-\zeta^2-1}\right)\cdot t} + A_2\cdot e^{\omega_n\left(-\zeta - \sqrt{-\zeta^2-1}\right)\cdot t}\end{aligned}\tag{4}

where

ωd= the damped natural frequency =ωn1ζ2\omega_d = \text{ the damped natural frequency } = \omega_n\sqrt{1-\zeta^2}

Because 0ζ0\leq\zeta\leq\infty, the radicand may be positive, negative, or zero. These three possibilities give rise to three different types of damped motion, and the difference in behaviour may be classified according to the table below

ζ\boldsymbol{\zeta}Type of Motionλ1,λ2\boldsymbol{\lambda_1,\lambda_2}Motion Characteristics
>1>1’Overdamped’Real, negative, distinctMotion decays to zero over a period of time, tt, with zero oscillation.
=1=1‘Critically Damped’Real, negative, equalMotion decays to zero over a period of time tt, where tt is the shortest possible time to return to equilibrium for this system.
<1<1‘Underdamped’ComplexMotion returns toward, and beyond zero in time tt where tt is smaller than the time for the critically damped case. Motion will overshoot and exhibit harmonic oscillation.

different values of ζ\zeta are demonstrated for free vibration of the mass-spring-damper system, with the maxima/minima of the underdamped cases shown, to show the change in the damped natural frequencies of the system. Furthermore, the period of the oscillatory systems are used to show the damped natural frequency which is then compared with the value from the formula. Hopefully this should elucidate the difference between the two sorts of natural frequencies.

Interactive demo

Most aircraft problems that we will look at exhibit underdamped motion ζ<1\zeta < 1 and Eq (4)1 can be re-written as:

x(t)=(A1ei1ζ2ωnt+A2ei1ζ2ωnt)eζωntx(t) = \left(A_1e^{i\sqrt{1-\zeta^2}\omega_n t} + A_2e^{-i\sqrt{1-\zeta^2}\omega_nt}\right)e^{-\zeta\omega_nt}

we introduce a parameter ωd=ωn1ζ2\omega_d=\omega_n\sqrt{1-\zeta^2}, which represents the \textsl{damped natural frequency}, which allows us to express this as

x(t)=(A3cosωdt+A4sinωdt)eζωnt(5)x(t) = \left(A_3\cos\omega_dt + A_4\sin\omega_dt\right)e^{-\zeta\omega_nt}\tag{5}

Where A3=A1+A2A_3=A_1+A_2 and A4=A1A2A_4 = A_1-A_2. Inspection of (5) shows an oscillation with angular frequency ωd\omega_d, and a decaying amplitude with time. This is demonstrated below:

To determine the system’s forced response to an input, we set F(t)F(t) to a given value, and solve. We will explore this via examples for real aircraft.

Determining Stability - Eigenvalues

We have already shown the characteristic equation of the mass spring damper system to be:

λ2+2ζωnλ+ωn2\lambda^2 + 2\zeta\omega_n\lambda + \omega_n^2

Note that if we had taken the Laplace Transform of the equation of motion with zero initial conditions, we would have arrived at the characteristic equation in terms of the Laplace Operator, ss

s2+2ζωns+ωn2=0(6)s^2 + 2\zeta\omega_n s + \omega_n^2 = 0\tag{6}

Clearly equations (2) and (6) are equivalent, hence the poles and eigenvalues are equivalent. Since we know that the poles of the system give information on its stability, then we can further surmise that its eigenvalues do the same. Furthermore, we have shown that the eigenvalues are given by the solution of the characteristic equation:

λ=ωn(ζ±ζ21)\lambda = \omega_n\left(-\zeta \pm \sqrt{-\zeta^2-1}\right)

hence we can see that the system’s eigenvalues and therefor its stability is a function of both the natural frequency and the damping.

For mechanical systems, damping can only be positive, and hence the term ζωn\zeta\omega_n will always be positive. This is a consequence of physical damping being a removal of kinetic energy - hence this always causes oscillations to decrease in amplitude. For aircraft, however, the interplay of aerodynamic and propulsive forcing causing a conversion between potential and kinetic energy can give rise to negative damping, whereby energy is added to a system. This causes oscillations which increase with time. For solution of the aircraft equations of motion, the eigenvalues may be real or complex, and have a positive or negative real component. This gives rise to four possibilities as shown in the table below:

λ\lambdaDynamic Response
Real and NegativeThis denotes an exponential decay or a convergence. This is indicative of static stability.
Real and PositiveThis denotes an exponential increase or a divergence. This is indicative of static instability.
Complex, (λ)<0\Re(\lambda) <0This denotes an oscillatory motion with a decreasing amplitude. This is indicative of dynamic stability.
Complex, (λ)>0\Re({\lambda})>0This denotes an oscillatory motion with an increasing amplitude. This is indicative of dynamic instability.

you will see from the above that the existence of a complex root gives rise to dynamic behaviour, and its sign does not matter.

Period and time to half amplitude

If we consider the general eigenvalue and decompose into its imaginary and real parts:

λ=(λ)+i(λ)\lambda = \Re(\lambda) + i\cdot\Im(\lambda)

Which gives the response

x=x0eλt=xee((λ)+i(λ))tx = x_0e^{\lambda t} = x_ee^{\left(\Re(\lambda) + i\cdot\Im(\lambda) \right)t}

For the oscillatory modes, the angular frequency is given by the imaginary part of the eigenvalue - hence the period of any oscillations, TT, is given by

T=2π(λ)T = \frac{2\pi}{|\Im(\lambda)|}

with a frequency given by

f=(λ)2πHzf= \frac{|\Im(\lambda)|}{2\pi} \text{Hz}

The decay rate of the oscillatory behaviour is given by the real part of the eigenvalue - hence this determines the damping. When looking at aircraft modes of motion, we tend to talk about the ‘time to half amplitude’ for convergent modes, or the ‘time to double amplitude’ for a divergent mode. Since the amplitude of the oscillation is given:

x=x0e(λ)t|x| = x_0e^{\Re(\lambda)t}

hence for a stable system, the time to half amplitude, t.5t_{.5}, is given by:

x02=x0e(λ)t.5    ln0.5=(λ)t.5t.5=0.69(λ)\begin{gathered} \frac{x_0}{2} = x_0e^{\Re(\lambda)t_{.5}}\\ \implies \ln|0.5| = \Re\left(\lambda\right)\cdot t_{.5}\\ \therefore t_{.5} = \frac{-0.69}{\Re({\lambda})}\end{gathered}

Note that as the time to half amplitude is only calculated for a stable system where the eigenvalue will be negative, the time will be a positive value. Similarly for an unstable system, the time to double amplitude t2t_2 is found from

2x0=x0e(λ)t2    ln2=(λ)t2t.5=ln(2)(λ)\begin{gathered} 2\cdot{x_0} = x_0e^{\Re(\lambda)t_{2}}\\ \implies \ln|2| = \Re\left(\lambda\right)\cdot t_{2}\\ \therefore t_{.5} = \frac{\ln(2)}{\Re({\lambda})}\end{gathered}

Again, since the time to double amplitude is only calculated for an unstable system where the eigenvalue will be positive, the time will always be a positive value.

Note that as presented above, we are assuming the dimensional equations of motion. If we are using them in their nondimensional form, then the reduced time is calculated and the answers must be multiplied by tt^*.

Damped and Undamped Natural Frequencies

The damped natural frequency, wdw_d, is the frequency at which the actual system will exhibit oscillatory behaviour if measured in the real world, subject to damping. The natural frequency, wnw_n, of the system is a frequency that the system would vibrate if the oscillatory behaviour were influenced by the system stiffness alone.

In short, we can get the two from:

ωd=(λ)\omega_d=\Im(\lambda) ωn=(λ)2+(λ)2\omega_n=\sqrt{\Re(\lambda)^2+\Im(\lambda)^2}

but it’s not always apparent why this is the case. We know that a second order system will have the general form of:

x¨+2ζωnx˙+ωn2=0\ddot{x} + 2\zeta\omega_n\dot{x} + \omega_n^2 = 0

and we have further shown that the eigenvalues found determine the modal characteristics of the system. Since these are the roots of the characteristic equation we can write the above as:

x¨+2ζωnx˙+ωn2=(x(Re(λ)+i(λ)))(x(Re(λ)i(λ)))\ddot{x} + 2\zeta\omega_n\dot{x} + \omega_n^2 = \left(x-\left(Re(\lambda)+i\Im(\lambda)\right)\right)\left(x-\left(Re(\lambda)-i\Im(\lambda)\right)\right)

the right hand side can be mutlipled out to give

λ22(λ)λ+(2(λ)+2(λ))\lambda^2 - 2\cdot\Re(\lambda){\lambda}+\left(\Re^2(\lambda)+\Im^2(\lambda)\right)

which, given the relationship between damped and undamped natural frequncies as

ωd=ωn1ζ2\omega_d=\omega_n\sqrt{1-\zeta^2}

may be readily rearranged to give the two expressions for natural frequency.

The stability of a multiple degree of freedom system

The theory presented above has been developed in a theory that is easily extensible to multiple degree of freedom systems. The general first order ODE can be represented as

x˙=Ax(7)\dot{\vec{x}} = \boldsymbol{A}\vec{x}\tag{7}

This equation represents a system of nn equations where nn is the number of states and degrees of freedom. Hence for the preceding analysis, n=1n=1, whereas it will be >1>1 for our aircraft problem. The general solution will be2:

xi=x0ieλti=1,,nx_i = x_{0_i}e^{\lambda t}\hspace{1cm}\forall i=1,\ldots,n

in vector form

x=x0eλt\vec{x}=\vec{x}_0e^{\lambda t}

similarly, differentiating and substituting into Eq. (7)

λx0eλt=Ax0eλt\lambda\vec{x}_0e^{\lambda t} = \boldsymbol{A}\vec{x}_0e^{\lambda t}

which is an eigenvalue (matrix/vector) problem and hence we can write

(λIA)x0=0(8)\left(\lambda\boldsymbol{I}-\boldsymbol{A}\right)\vec{x}_0=0\tag{8}

The trivial solution is zero displacement, x0=0x_0=0, but more useful solution is given by the determinant

λIA=0(9)\left|\lambda\boldsymbol{I}-\boldsymbol{A}\right|=0\tag{9}

Equation (9) is the characteristic equation of a multiple degree of freedom system, and will be a polynomial of order nn, yielding nn eigenvalues dictating the system stability. If any single value is unstable, then the aircraft is deemed dynamically unstable. Similarly, if any of the eigenvalues give a complex pair, then the characteristic motion will be oscillatory as the equations are coupled.

Once the eigenvalues have been solved, each of these may be substituted into Eq. (8) to solve for the vector of amplitudes to give nn vectors of nn amplitudes. The magnitudes of each respect state motion help us to determine which variables correspond to which eigenvalue (or which eigenvalues correspond to which mode of the aircraft).

Before we can utilise the eigenvalues to tell us information about the modes of the aircraft, we will require an understanding of what each of the aircraft modes of motion comprises. The following section will describe what each mode is, in a physical sense before we delve into the governing mathematics.

Footnotes

  1. ζ21=i1ζ2\sqrt{\zeta^2-1}=i\sqrt{1-\zeta^2} and eaeb=e(a+b)e^ae^b=e^{(a+b)}

  2. If you’ve not seen the symbol \forall it means ‘for each’ or ‘for all’ dependent upon usage. It’s an upside-down ‘A’. Not to be confused with what Charles Sanders Peirce denoted it to mean.