16.333: Lecture # 6 Aircraft Longitudinal Dynamics ? Typical aircraft open-loop motions ? Longitudinal modes ? Impact of actuators ? Linear Algebra in Action! ? ? Fall 2004 16.333 5–1 Longitudinal Dynamics ? Recall: X denotes the force in the X-direction, and similarly for Y and Z, then (as on 4–13) ?X , . . . X u ≡ ?u 0 ? Longitudinal equations (see 4–13) can be rewritten as: mu˙ = X u u + X w w ? mg cos Θ 0 θ + ΔX c m( ˙ w ? qU 0 ) = Z u u + Z w w + Z w˙ w˙ + Z q q ? mg sin Θ 0 θ + ΔZ c I yy q˙ = M u u + M w w + M w˙ w˙ + M q q + ΔM c There is no roll/yaw motion, so q = θ ˙ .? ? Control commands ΔX c , ΔZ c , and ΔM c have not yet been speci?ed. Fall 2004 16.333 5–2 ? Rewrite in state space form as ???????? mu˙ X u X w 0 ?mg cos Θ 0 u ΔX c ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ΔZ c (m ? Z w˙ ) ˙w Z u Z w Z q + mU 0 ?mg sin Θ 0 M u M w M q 0 w q + = ΔM c ?M w˙ w˙ + I yy q˙ θ ˙ 0 0 1 0 θ 0 ???? m 0 0 0 u˙ ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? w˙ q˙ 0 m ? Z w˙ 0 0 0 ?M w˙ I yy 0 θ ˙ 0 0 0 1 ?????? X u X w 0 ?mg cos Θ 0 u ΔX c ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ΔZ c Z u Z w Z q + mU 0 ?mg sin Θ 0 M u M w M q 0 w q + = ΔM c 0 0 1 0 θ 0 ˉ E ˙ = AX + ? descriptor state space form cX ˉ = E ?1 (AX + ?c) = AX + c?X ˙ ? ? ? ? ? ? ? ? ? ? ? ? ? ? Fall 2004 16.333 5–3 ? Write out in state space form: ? X u X w 0 ?g cos Θ 0 m m Z q + mU 0 ?mg sin Θ 0 Z u Z w m ? Z w˙ m ? Z w˙ m ? Z w˙ m ? Z w˙ I ?1 [M u + Z u Γ] I ?1 [M w + Z w Γ] I ?1 [M q + (Z q + mU 0 )Γ] yy mg sin Θ 0 Γ?I ?1 yy yy yy 0 0 1 0 ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? A = M w˙ Γ = m ? Z w˙ ? Note: slight savings if we de?ned symbols to embed the mass/inertia ? ? X u = X u /m, Z ? u = Z u /m, and M q = M q /I yy then A matrix collapses to: ?? ? ? 0 ?g cos Θ 0 X u X w A ? = ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? Z ? u Z ? w ?g sin Θ 0 Z ? q + U 0 1 ? Z ? w˙ 1 ? Z ? w˙ 1 ? Z ? w˙ 1 ? Z ? w˙ ? ? ? ? ? M u + Z ? u Γ M w + Z ? w Γ M q + ( Z ? q + U 0 ) ?? Γ ?g sin Θ 0 Γ 0 0 1 0 ? ? M w˙ Γ = 1? Z ? w˙ ? Check the notation that is being used very carefully ? To ?gure out the c vector, we have to say a little more about how the control inputs are applied to the system. Fall 2004 16.333 5–4 Longitudinal Actuators ? Primary actuators in longitudinal direction are the elevators and thrust. – Clearly the thrusters/elevators play a key role in de?ning the steady-state/equilibrium ?ight condition – Now interested in determining how they also in?uence the aircraft motion about this equilibrium condition de?ect elevator → u(t), w(t), q(t), . . . ? Recall that we de?ned ΔX c as the perturbation in the total force in the X direction as a result of the actuator commands – Force change due to an actuator de?ection from trim ? Expand these aerodynamic terms using same perturbation approach ΔX c = X δ e δ e + X δ p δ p – δ e is the de?ection of the elevator from trim (down positive) – δ p change in thrust – X δ e and X δ p are the control stability derivatives ? ? ? ? Fall 2004 16.333 5–5 Now we have that ? ???? ΔX c ? ? ? ? = E ?1 ? ? ? ? X δ e X δ p Z δ e Z δ p M δ e M δ p ? ? ? ? c = E ?1 ? ? ? ? ΔZ c ΔM c δ e = Bu δ p 0 0 0 ? For the longitudinal case ?? X δ e X δ p m m Z δ e Z δ p m ? Z ˙w m ? Z ˙w B = ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? I ?1 I ?1 [M δ e + Z δ e Γ] yy yy M δ p + Z δ p Γ 0 0 ? Typical values for the B747 X δ e = ?16.54 X δ p = 0.3mg = 849528 Z δ e = ?1.58 · 10 6 Z δ p ≈ 0 M δ e = ?5.2 · 10 7 M δ p ≈ 0 ? Aircraft response y = G(s)u X ˙ = AX + Bu → G(s) = C(sI ? A) ?1 B y = CX ? We now have the means to modify the dynamics of the system, but ?rst let’s ?gure out what δ e and δ p really do. ? ? ? ? Fall 2004 16.333 5–6 Longitudinal Response ?? Final response to a step input u = u/s, y = G(s)u, use the FVT u? lim y(t) = lim s G(s) s 0 s t→∞ → lim y(t) = G(0)? u? t→∞ u = ?(CA ?1 B)? ? Initial response to a step input, use the IVT u? y(0 + ) = lim s G(s) = lim G(s)?u s→∞ s s→∞ – For your system, G(s) = C(sI ? A) ?1 B + D, but D ≡ 0, so lim G(s) → 0 s→∞ – Note: there is NO immediate change in the output of the motion variables in response to an elevator input ? y(0 + ) = 0 ? Consider the rate of change of these variables y˙(0 + ) y˙(t) = C ˙ = CAX + CBu X and normally have that CB =? 0. Repeat process above 1 to show that y˙(0 + ) = CBu?, and since C ≡ I, y˙(0 + ) = Bu? ? Looks good. Now compare with numerical values computed in Mat- lab. Plot u, α, and ?ight path angle γ = θ ? α (since Θ 0 = γ 0 = 0 – see picture on 4–8) CB 1 Note that C(sI ? A) ?1 B + D = D + + CA ?1 B + ... 2 s s Fall 2004 16.333 5–7 Elevator (1 ? elevator down – stick forward) ? See very rapid response that decays quickly (mostly in the ?rst 10 seconds of the α response) ? Also see a very lightly damped long period response (mostly u, some γ, and very little α). Settles in >600 secs ? Predicted steady state values from code: 14.1429 m/s u (speeds up) -0.0185 rad α (slight reduction in AOA) -0.0000 rad/s q -0.0161 rad θ 0.0024 rad γ – Predictions appear to agree well with the numerical results. – Primary result is a slightly lower angle of attack and a higher speed ? Predicted initial rates of the output values from code: -0.0001 m/s 2 u˙ -0.0233 rad/s α˙ -1.1569 rad/s 2 q˙ 0.0000 rad/s θ ˙ 0.0233 rad/s γ˙ – All outputs are zero at t = 0 + , but see rapid changes in α and q. – Changes in u and γ (also a function of θ) are much more gradual – not as easy to see this aspect of the prediction ? Initial impact Change in α and q (pitches aircraft) ? Long term impact Change in u (determines speed at new equilib- rium condition) Fall 2004 16.333 5–8 0 200 400 600 05 1015202530 u time 0 200 400 600 ?0.03 ?0.025 ?0.02 ?0.015 ?0.01 ?0.005 0 a (rad) time Step response to 1 deg elevator perturbation 0 200 400 600 ?0.1 ?0.05 0 0.05 0.1 g time 0 10 20 30 40 05 1015202530 u time 0 10 20 30 40 ?0.03 ?0.025 ?0.02 ?0.015 ?0.01 ?0.005 0 a (rad) time 0 10 20 30 40 ?0.1 ?0.08?0.06?0.04?0.02 0 0.02 g time Figure 1: Step Response to 1 deg elevator perturbation – B747 at M=0.8 Fall 2004 16.333 5–9 Thrust (1/6 input) ? Motion now dominated by the lightly damped long period response ? Short period motion barely noticeable at beginning. ? Predicted steady state values from code: 0 m/s u 0 rad α 0 rad/s q 0.05 rad θ 0.05 rad γ – Predictions appear to agree well with the simulations. – Primary result – now climbing with a ?ight path angle of 0.05 rad at the same speed we were going before. ? Predicted initial rates of the output values from code: 2.9430 m/s 2 u˙ 0 rad/s α˙ 0 rad/s 2 q˙ 0 rad/s θ ˙ 0 rad/s γ˙ – Changes to α are very small, and γ response initially ?at. – Increase power, and the aircraft initially speeds up ? Initial impact Change in u (accelerates aircraft) ? Long term impact Change in γ (determines climb rate) Fall 2004 16.333 5–10 0 200 400 600 ?15?10 ?5 05 1015 u time 0 200 400 600 ?0.02 ?0.015 ?0.01 ?0.005 0 0.005 0.01 0.015 0.02 a (rad) time Step response to 1/6 thrust perturbation 0 200 400 600 0 0.020.040.060.08 0.1 g time 0 10 20 30 40 01234567 u time 0 10 20 30 40 ?5 05 101520 x 10 ?4 a (rad) time 0 10 20 30 40 0 0.010.020.030.040.050.060.070.080.09 g time Figure 2: Step Response to 1/6 thrust perturbation – B747 at M=0.8 Fall 2004 16.333 5–11 Frequency Domain Response ? Plot and inspect transfer functions from δ e and δ p to u, w, and γ – See following pages From elevator: ? – Huge response at the phugoid mode for both u and γ (very lightly damped) – Short period mode less pronounced – Response falls o? very rapidly – Response to w shows a pole/zero cancelation (almost) of the phugoid mode. So the magnitude level is essentially constant out to the frequency of the short period mode Why would we expect that? From thrust: ? – Phugoid peaks present, but short period mode is very weak (not in u, low in γ, w). ? entirely consistent with the step response. – Thrust controls speed (initially), so we would expect to see a large response associated with the phugoid mode (speed variations are a key component of this mode) Fall 2004 16.333 5–12 10 ?2 10 ?1 10 0 10 0 10 1 10 2 10 3 10 4 |G u d e | Freq (rad/sec) 10 ?2 10 ?1 10 0 10 ?1 10 0 10 1 10 2 10 3 10 4 |G a d e | Freq (rad/sec) Transfer function from elevator to flight variables 10 ?2 10 ?1 10 0 10 ?1 10 0 10 1 10 2 10 3 10 4 |G g d e | Freq (rad/sec) 10 ?2 10 ?1 10 0 ?350?300?250?200?150?100 ?50 0 arg G u d e Freq (rad/sec) 10 ?2 10 ?1 10 0 ?350?300?250?200?150?100 ?50 0 arg G a d e Freq (rad/sec) 10 ?2 10 ?1 10 0 ?350?300?250?200?150?100 ?50 0 arg G g d e Freq (rad/sec) Figure 3: TF’s from elevator to ?ight variables – B747 at M=0.8 Fall 2004 16.333 5–13 10 ?2 10 ?1 10 0 10 ?2 10 ?1 10 0 10 1 10 2 10 3 |G u d p | Freq (rad/sec) 10 ?2 10 ?1 10 0 10 ?2 10 ?1 10 0 10 1 10 2 10 3 |G a d p | Freq (rad/sec) Transfer function from thrust to flight variables 10 ?2 10 ?1 10 0 10 ?2 10 ?1 10 0 10 1 10 2 10 3 |G g d p | Freq (rad/sec) 10 ?2 10 ?1 10 0 ?150?100 ?50 0 50 100150 arg G u d p Freq (rad/sec) 10 ?2 10 ?1 10 0 ?150?100 ?50 0 50 100150 arg G a d p Freq (rad/sec) 10 ?2 10 ?1 10 0 ?350?300?250?200?150?100 ?50 0 arg G g d p Freq (rad/sec) Figure 4: TF’s from thrust to ?ight variables– B747 at M=0.8 Fall 2004 16.333 5–14 ? Summary: – To increase equilibrium climb rate, add power. – To increase equilibrium speed, increase δ e (move elevator further down). – Transient (initial) e?ects are the opposite and tend to be more consistent with what you would intuitively expect to occur ? ? ? ? ? Fall 2004 16.333 5–15 Modal Behavior ? Analyze model of vehicle dynamics to quantify the responses seen. – Homogeneous dynamics of the form X ˙ = AX, so the response is X(t) = e At X(0) – a matrix exponential. ? To simplify the investigation of the system response, ?nd the modes of the system using the eigenvalues and eigenvectors – λ is an eigenvalue of A if det(λI ? A) = 0 which is true i? there exists a nonzero v (eigenvector) for which (λI ?A)v = 0 Av = λv? – If A (n ×n), typically get n eigenvalues/eigenvectors Av i = λ i v i – Assuming that eigenvectors are linearly independent, can form ? ? λ 1 0 A v 1 = v 1 ? ··· v n ··· v n ? . . . 0 λ n A T = T Λ T ?1 AT = Λ , A = T ΛT ?1 ? 1 – Given that e At = I + At + 2! (At) 2 + . . ., and that A = T ΛT ?1 , then it is easy to show that n X(t) = e At X(0) = T e Λt T ?1 X(0) = v i e λ i t β i i=1 – State solution is a linear combination of system modes v i e λ i t e λ i t – determines nature of the time response v i – gives extent to which each state participates in that mode β i – determines extent to which initial condition excites the mode Fall 2004 16.333 5–16 ? The total behavior of the system can be found from the system modes ? Consider numerical example of B747 ? ? A = ? ? ? ? ?0.0069 0.0139 0 ?9.8100 ?0.0905 ?0.3149 235.8928 0 0.0004 ?0.0034 ?0.4282 0 0 0 1.0000 0 ? ? ? ? which gives two sets of complex eigenvalues λ = ?0.3717 ± 0.8869 i, ω = 0.962, ζ = 0.387, short period λ = ?0.0033 ± 0.0672 i, ω = 0.067, ζ = 0.049, Phugoid - long period – Result is consistent with step response - heavily damped fast response, and a lightly damped slow one. ? To understand eigenvectors, must do some normalization (scales each element appropriately so that we can compare relative sizes) – u? = u/U 0 , α = w/U 0 , q? = q/(2U 0 /c) – Then divide through so that θ ≡ 1 Short Period Phugoid u? 0.0156 + 0.0244 i ?0.0254 + 0.6165 i α 1.0202 + 0.3553 i 0.0045 + 0.0356 i q? ?0.0066 + 0.0156 i ?0.0001 + 0.0012 i θ 1.0000 1.0000 ? Short Period – primarily θ and α = w? in the same phase. The u? and q? response is very small. ? Phugoid – primarily θ and u?, and θ lags by about 90 ? . The α and q? response is very small. ? Dominant behavior agrees with time step responses – note how initial conditions were formed. Fall 2004 16.333 5–17 1 2 30 210 60 240 90 270 120 300 150 330 180 0 0.96166 0.5 1 30 210 60 240 90 270 120 300 150 330 180 0 0.067282 0 5 10 15 ?0.5 0 0.5 1 Perturbation States u,w,q time (sec) u w q 0 100 200 300 400 500 600 ?1 ?0.5 0 0.5 1 Perturbation States u,w,q time (sec) u w q Figure 5: Mode Response – B747 at M=0.8 Fall 2004 16.333 5–19 ? Relative motion between aircraft and an observer ?ying at a constant speed U 0 t ? Motion of perturbed aircraft with respect to an unperturbed one ? Note phasing of the forward velocity x˙ e with respect to altitude z e – aircraft faster than observer at the bottom, slower at the top – The aircraft speeds up and slows down – leads and lags the ob- server. ? Consistent with ?ight path? ? Consistent with Lanchester’s approximation on 4–1? Fall 2004 16.333 5–20 Summary ? Two primary longitudinal modes: phugoid and short-period – Have versions from the full model – but can develop good approx- imations that help identify the aerodynamic features that deter- mine the mode frequencies and damping Impact of the various actuators clari?ed: – Short time-scale – Long time-scale Fall 2004 16.333 5–21 Matrix Diagonalization ? Suppose A is diagonizable with independent eigenvectors V = [v 1 , . . . , v n ] – use similarity transformations to diagonalize dynamics matrix x˙ = Ax x˙ d = A d x d ? ? ? λ 1 . . . ? Δ V ?1 AV = ? = Λ = A d λ n – Corresponds to change of state from x to x d = V ?1 x ? System response given by e At , look at power series expansion At = V ΛtV ?1 2 (At) 2 = (V ΛtV ?1 )V ΛtV ?1 = V Λ t 2 V ?1 ? (At) n = V Λ n t n V ?1 1 At e = I + At + (At) 2 + . . . ? 2 ? = V I + Λ + 1 Λ 2 t 2 + . . . V ?1 2 ? ? λ 1 t e . . . ? V ?1 = V e Λt V ?1 = V ? λ n t e Fall 2004 16.333 5–22 ? Taking Laplace transform, ? ? 1 ? s?λ 1 ? (sI ? A) ?1 = V ? . . . ? V ?1 1 s?λ n n ? R i = s ? λ i i=1 where the residue R i = v i w i T , and we de?ne ? ? T ? ? w 1 . V = v 1 . . . v n , V ?1 = ? . . ? T w n ? Note that the w i are the left eigenvectors of A associated with the right eigenvectors v i ? ? ? ? λ 1 λ 1 AV = V ? . . . ? V ?1 A = ? . . . ? V ?1 ? λ n λ n ? ? ? ?? ? T w 1 λ 1 w 1 T . . ? . ? A = ? . . . ?? . ? . . T T w λ n w nn T where w i A = λ i w T i ? ? Fall 2004 16.333 5–23 ? So, if x˙ = Ax, the time domain solution is given by n x(t) = e λ i t v i w T i x(0) dyad i=1 n ? x(t) = [w T i x(0)]e λ i t v i i=1 ? The part of the solution v i e λ i t is called a mode of a system – solution is a weighted sum of the system modes – weights depend on the components of x(0) along w i ? Can now give dynamics interpretation of left and right eigenvectors: Av i = λ i v i , w i A = λ i w i , w i T v j = δ ij so if x(0) = v i , then n T x(t) = (w i x(0))e λ i t v i i=1 λ i t = e v i ? so right eigenvectors are initial conditions that result in relatively simple motions x(t). With no external inputs, if initial condition only disturbs one mode, then the response consists of only that mode for all time. ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? Fall 2004 16.333 5–24 ? If A has complex conjugate eigenvalues, the process is similar but a little more complicated. ? Consider a 2x2 case with A having eigenvalues a± bi and associated eigenvectors e 1 , e 2 , with e 2 = eˉ 1 . Then ? a + bi 0 ? ?1 A = e 1 e 2 e 1 e 2 0 a ? bi ? a + bi 0 ? = e 1 eˉ 1 ?1 ≡ TDT ?1 0 a ? bi eˉ 1 e 1 Now use the transformation matrix 1 ?i M ?1 = 1 1 M = 0.5 1 i i ?i Then it follows that ? A = TDT ?1 = (TM)(M ?1 DM)(M ?1 T ?1 ) = (TM)(M ?1 DM)(TM) ?1 which has the nice structure: ? a b ? A = Re(e 1 ) Im(e 1 ) Re(e 1 ) Im(e 1 ) ?1 ?b a where all the matrices are real. ? With complex roots, the diagonalization is to a block diagonal form. ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?? ? ? ? ? ? ? ? Fall 2004 16.333 5–25 For this case we have that ? cos(bt) sin(bt) ? At e = Re(e 1 ) Im(e 1 ) e at Re(e 1 ) Im(e 1 ) ?1 ? sin(bt) cos(bt) ? Note that Re(e 1 ) Im(e 1 ) ?1 is the matrix that inverts ? Re(e 1 ) Im(e 1 ) ? ?1 ? ? 1 0 Re(e 1 ) Im(e 1 ) Re(e 1 ) Im(e 1 ) = 0 1 ? So for an initial condition to excite just this mode, can pick x(0) = [Re(e 1 )], or x(0) = [Im(e 1 )] or a linear combination. ? Example x(0) = [Re(e 1 )] ? cos(bt) sin(bt) x(t) = e At x(0) = Re(e 1 ) Im(e 1 ) e at ? sin(bt) cos(bt) · Im(e 1 ) ?1 [Re(e 1 )]Re(e 1 ) ? cos(bt) sin(bt) 1 = Re(e 1 ) Im(e 1 ) e at ? sin(bt) cos(bt) 0 ? cos(bt) at = e Re(e 1 ) Im(e 1 ) ? sin(bt) at = e (Re(e 1 ) cos(bt) ? Im(e 1 ) sin(bt)) which would ensure that only this mode is excited in the response Fall 2004 16.333 5–26 Example: Spring Mass System ? Classic example: spring mass system consider simple case ?rst: m i = 1, and k i = 1 M1 k1 k2 k3 k4 k5 M3 M2 Z1 Z3 Z2 ? ? x = z 1 ? z 2 z 3 ˙z 1 ˙z 2 ? ˙z 3 0 I A = ?M ?1 K ? 0 M = diag(m i ) ? K = ? k 1 + k 2 + k 5 ?k 5 ?k 2 ?k 5 k 3 + k 4 + k 5 ?k 3 ?k 2 ?k 3 k 2 + k 3 ? ? Eigenvalues and eigenvectors of the undamped system λ 1 = ±0.77i λ 2 = ±1.85i λ 3 = ±2.00i v 1 v 2 v 3 1.00 1.00 1.00 1.00 1.41 ±0.77i ±0.77i ±1.08i 1.00 ?1.41 ±1.85i ±1.85i ?2.61i ?1.00 0.00 ±2.00i ?2.00i 0.00 - - - Fall 2004 16.333 5–27 ? Initial conditions to excite just the three modes: x i (0) = α 1 Re(v i ) + α 2 Im(v 1 ) ?α j ∈ R – Simulation using α 1 = 1, α 2 = 0 ? Visualization important for correct physical interpretation ? Mode 1 λ 1 = ±0.77i M 1 M 3 M 2 – Lowest frequency mode, all masses move in same direction – Middle mass has higher amplitude motions z 3 , motions all in phase 0 1 2 3 4 5 6 7 8 9 10 ?0.8 ?0.6 ?0.4 ?0.2 0 0.2 0.4 0.6 time displacement V1 z1 z2 z3 Fall 2004 16.333 5–28 ? Mode 2 λ 2 = ±1.85i M 1 M 3 M 2 -  - – Middle frequency mode has middle mass moving in opposition to two end masses – Again middle mass has higher amplitude motions z 3 0 1 2 3 4 5 6 7 8 9 10 ?0.4 ?0.3 ?0.2 ?0.1 0 0.1 0.2 0.3 0.4 time displacement V3 z1 z2 z3 Fall 2004 16.333 5–29 ? Mode 3 λ 3 = ±2.00i M 1 M 3 M 2 - 0  – Highest frequency mode, has middle mass stationary, and other two masses in opposition 0 1 2 3 4 5 6 7 8 9 10 ?0.4 ?0.3 ?0.2 ?0.1 0 0.1 0.2 0.3 0.4 time displacement V2 z1 z2 z3 ? Eigenvectors with that correspond with more constrained motion of the system are associated with higher frequency eigenvalues