16.333: Lecture # 14

16.333: Lecture # 14 Equations of Motion in a Nonuniform Atmosphere Gusts and Winds 1 Fall 2004 16.333 12–2 Equations of Motion • Analysis to date has assumed that the atmosphere is calm and fixed – Rarely true since we must contend with gusts and winds – Need to understand how these air motions impact our modeling of the aircraft. • Must modify aircraft equations of motion since the aerodynamic forces and moments are functions of the relative motion between the aircraft and the atmosphere, and not of the inertial velocities. – Thus the LHS of the dynamics equations (F� = m�a) must be written in terms of the velocities relative to the atmosphere. – If u is the aircraft perturbation velocity (X direction), and ug is the gust velocity in that direction, then the aircraft velocity with respect to the atmosphere is ua = u − ug • Now rewrite aerodynamic forces and moments in terms of aircraft velocity with respect to the atmosphere (see 4–11) ΔX = ∂X ∂X ∂X (w − wg ) + (u − ug ) + (w˙ − ẇg ) ˙ ∂U ∂W ∂W ∂X ∂X ∂X g + (q − qg ) + . . . + θ + ... + θ + ΔX c ∂Q ∂Θ ∂Θ – The gravity terms ∂X g ∂Θ and control terms ΔX c = Xδe δe + Xδp δp stay the same. Fall 2004 16.333 12–3 • The rotation gusts pg , qg , and rg are caused by spatial variations in the gust components ⇒ rotary gusts are related to gradients of the vertical gust field pg = ∂wg ∂wg and qg = ∂y ∂x -wg Gust field creating an effective rolling gust. +wg y Equivalent distribution -wg +w g Equivalent to velocity created by a pitching motion Figure 1: Gust Field creating an effective pitching gust. Fall 2004 16.333 12–4 • The next step is to include these new forces and moments in equations of motion ⎤⎡ ⎤ ⎡ ⎡ u̇ 0 −mg cos Θ0 m 0 0 0 Xu Xw ⎥⎢ ⎥ ⎢ ⎢ ⎢ 0 m − Zẇ 0 0 ⎥⎢ ẇ ⎥ ⎢ Zu Zw Zq + mU0 −mg sin Θ0 ⎥⎢ ⎥ = ⎢ ⎢ Mq 0 ⎣ 0 −Mẇ Iyy 0 ⎦⎣ q̇ ⎦ ⎣ Mu Mw 0 0 0 1 0 0 1 0 θ̇ ⎡ ⎤ Xδe Xδp � � ⎢ ⎥ ⎢ Zδe Zδp ⎥ δe + ⎢ ⎥ ⎣ Mδe Mδp ⎦ δt 0 0 ⎤ ⎡ ⎡ ⎤ −Xu −Xw 0 u ⎢ ⎥ 0 ⎥⎣ g ⎦ ⎢ −Zu −Zw + ⎢ ⎥ wg ⎣ −Mu −Mw −Mq ⎦ qg 0 0 0 the ⎤⎡ u ⎥⎢ ⎥⎢ w ⎥⎢ ⎦⎣ q θ ˆ + Bˆ uu + Bˆ w w ⇒Eẋ = Ax • Multiply through by E −1 to get new state space model ẋ = Ax + Buu + Bw w which has both control u and disturbance w inputs. – A similar operation can be performed for the lateral dynamics in terms of the disturbance inputs vg , pg , and rg . • Can now compute the response to specific types of gusts, such as a step or sinusoidal function, but usually are far more interested in the response to a stochastic gust field ⎤ ⎥ ⎥ ⎥ ⎦ Fall 2004 16.333 12–5 Atmospheric Turbulence • Can develop the best insight to how aircraft will behave with gust disturbances if we treat the disturbances as random processes. – What is a random process? Something (signal) that is random so that a deterministic description is not practical!! – But we can often describe the basic features of the process (e.g. mean value, how much it varies about the mean). • Atmospheric turbulence is a random process, and the magnitude of the gust can only be described in terms of statistical properties. – For a random process f (t), talk about the mean square � T 1 f 2(t) = lim f 2(t) dt T →∞ T 0 as a measure of the disturbance intensity (how strong it is). • Signal f (t) can be decomposed into its Fourier components, so can use that to develop a frequency domain measure of disturbance strength – Φ(ω) ≈ that portion of f 2(t) that occurs in the frequency band ω → ω + dω – Φ(ω) is called the power spectral density • Bottom line: For a linear system y = G(s)w, then Φy (ω) = Φw (ω)|G(jω)|2 ⇒ Given an input disturbance spectral density (e.g. gusts), quite sim­ ple to predict expected output (e.g. ride comfort, wing loading). Fall 2004 16.333 12–6 Implementing a PSD in Matlab • Simulink has a Band­Limited White Noise block that can be used for continuous systems. – Primary difference from the Random Number block is that this block produces output at a specific sample rate, which is related to the correlation time of the noise. • Continuous white noise has a correlation time of 0 ⇒ flat power spectral density (PSD), and a covariance of infinity. – Non­physical, but a useful approximation when the noise distur­ bance has a correlation time that is small relative to the natural bandwidth of the system. – Can simulate effect of white noise by using a random sequence with a correlation time much smaller than the shortest time con­ stant of the system. • Band­Limited White Noise block produces such a sequence where the correlation time of the noise is the sample rate of the block. – For accurate simulations, use a correlation time tc much smaller than the fastest dynamics fmax of the system. 1 2π tc ≈ 100 fmax Fall 2004 16.333 12–7 • Power spectral densities for Von Karman model (see MIL­F­8785C) �−5/6 2Lu � 1 + (1.339LuΩ)2 Φug (Ω) = σu2 π where σu is the intensity measure of the disturbance, Ω is the spatial frequency variable, and Lu is a length scale of the disturbance. – All scale parameters depend on day, altitude, and turbulence type. – Different models available for each direction. – Scale length parameter L (in feet) varies with altitude – MIL­F­ 8785C model valid up to 1000 feet Lu = h (0.177 + 0.000823h)1.2 – Turbulence intensity for low altitude flight MIL­F­8785C as σu = 0.1W20 (0.177 + 0.000823h)0.4 3 where W20 is the wind speed as measured at 20 ft 3 W20 < 15 knots is classified as “light” turbulence W20 ≈ 30 knots is “moderate” W20 > 45 knots is “heavy” Fall 2004 16.333 12–8 – Standard approach: assume turbulence field fixed (frozen) in space ⇒ aircraft is just flying into it. Vertical gust Longitudinal gust Lateral gust 10 ft/sec 10 ft/sec Atmospheric gust profiles 10 ft/sec 3 Then can relate turbulence spatial frequency Ω to the temporal frequency ω that the aircraft would feel ω Ω≡ ⇒ ω = ΩU0 U0 – Used to model how the aircraft will behave in bumpy air. 300 ft on radio altimeter Figure 2: Wind gust examples 5 sec Touchdown Fall 2004 16.333 12–9 Wind Shear • Wind shear is of special interest because it involves significant local changes in the vertical and horizontal velocity (e.g. downdraft) – Particularly important near airports during landings. • Simple analysis – consider the case shown, where the change in (hor­ izontal) wind velocity is represented by � � du ug = Δh dh where – Δh corresponds to changes in altitude (what we previously just called h) – Value du/dh gives magnitude of wind shear, which is 0.08–0.15s−1 for moderate and 0.15–0.2s−1 for strong. Mean wind Flight path Glide slope Flare Figure 3: Aircraft descending into a horizontal windshear • But can write the changes in altitude (use ḣ = Δh/Δt) as ḣ = U0(θ − α) Fall 2004 16.333 12–10 • Now have a coupled situation that is quite interesting – Forces on the aircraft change with altitude (h) because height changes the gust velocity – Changes in the forces on the aircraft impact the height of the airplane since both θ and α will change. • Analyze the coupling by looking at the longitudinal equations � �T x = u w q θ h � � h = 0 0 0 0 1 x = Chx with controls fixed (zero) and input ug ˜ + B˜ w (:, 1)ug = Ax ˜ + B˜ w (:, 1) du h ẋ = Ax dh ˜ + du B˜ w (:, 1)Chx = Ax dh � � du = Ã + B˜ w (:, 1)Ch x dh ⇒ dynamics have been modified because of the coupling between the change in altitude and the change in forces (with ug ). • Closed this loop on the B747 dynamics to obtain du dh Phugoid Poles 0 ­0.0033 ± 0.0672i 0.08 ­0.0014 ± 0.1150i 0.15 0.0002 ± 0.1442i 0.2 0.0014 ± 0.1619i – Clearly this coupling is not good, and an unstable Phugoid mode is to be avoided during landing operations. Fall 2004 16.333 12–11 Wind Code 1 2 3 % Gust modeling % 16.333, Fall 2004 % Jonathan P. How 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 Xu=­1.982e3;Xw=4.025e3;Zu=­2.595e4;Zw=­9.030e4;Zq=­4.524e5;Zwd=1.909e3; Mu=1.593e4;Mw=­1.563e5;Mq=­1.521e7;Mwd=­1.702e4; % g=9.81;theta0=0;S=511;cbar=8.324; U0=235.9;Iyy=.449e8;m=2.83176e6/g;cbar=8.324;rho=0.3045; Xdp=.3*m*g;Zdp=0;Mdp=0; Xde=­3.818e­6*(1/2*rho*U0^2*S);Zde=­0.3648*(1/2*rho*U0^2*S); Mde=­1.444*(1/2*rho*U0^2*S*cbar);; % Ehat=[m 0 0 0;0 m­Zwd 0 0;0 ­Mwd Iyy 0;0 0 0 1]; Ahat=[Xu Xw 0 ­m*g*cos(theta0);[Zu Zw Zq+m*U0 ­m*g*sin(theta0)]; [Mu Mw Mq 0];[ 0 0 1 0]]; Bhat=[Xde Xdp;Zde Zdp;Mde Mdp;0 0]; % % form the gust input matrix Bwhat=[­Xu ­Xw 0;­Zu ­Zw 0;­Mu ­Mw ­Mq; 0 0 0]; % % add height state Ehat(5,5)=1; % \dot h state not coupled Ahat(5,5)=0;Ahat(5,[1:4])=[0 ­1 0 U0]; % add height state Bhat(5,2)=0; % noinput Bwhat(5,3)=0; % no input 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 % form the full model % E \dot x = \hat A x + \hat B u_controls + \hat B_w w % ==> \dot x = A x + B u_controls + B_w w % A=inv(Ehat)*Ahat; B=inv(Ehat)*Bhat; Bw=inv(Ehat)*Bwhat; % % set u_controls=0 % assume that w_g=q_g=0 and % u_g = (du/dh) h % du_dh=[0 .08 .15 .2]; A1=A+Bw(:,1)*[0 0 0 0 1]*du_dh(1); A2=A+Bw(:,1)*[0 0 0 0 1]*du_dh(2); A3=A+Bw(:,1)*[0 0 0 0 1]*du_dh(3); A4=A+Bw(:,1)*[0 0 0 0 1]*du_dh(4); 45 46 47 48 ev1=eig(A1);ev2=eig(A2);ev3=eig(A3);ev4=eig(A4); plot([ev1 ev2 ev3 ev4],’x’) [ev1 ev2 ev3 ev4]