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 ?xed – 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. F– Thus the LHS of the dynamics equations ( ? = m?a) must be written in terms of the velocities relative to the atmosphere. – If u is the aircraft perturbation velocity (X direction), and u g is the gust velocity in that direction, then the aircraft velocity with respect to the atmosphere is u a = u ? u g ? Now rewrite aerodynamic forces and moments in terms of aircraft velocity with respect to the atmosphere (see 4–11) ?X ?X ?X ΔX = (u ? u g ) + (w ? w g ) + ? ˙ ( ˙ w˙ g ) ?U ?W W w ? ?X ?X ?X g + (q ? q g ) + . . . + θ + . . . + θ + Δ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 p g , q g , and r g are caused by spatial variations in the gust components ? rotary gusts are related to gradients of the vertical gust ?eld ?w g ?w g p g = and q g = ?y ?x g g +w g +w g y Gust field creating gust. -w -w an effective rolling Equivalent distribution Equivalent to velocity created by a pitching motion Figure 1: Gust Field creating an e?ective pitching gust. ? ? Fall 2004 16.333 12–4 ? The next step is to include these new forces and moments in the equations of motion ???????? 0 0 0 u˙ X u X w 0 ?mg cos Θ 0 m u ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? 0 m ? Z ˙w 0 0 0 ?M ˙w I yy 0 0 0 0 1 w˙ q˙ Z u Z w Z q + mU 0 ?mg sin Θ 0 0M u M w M q w q = θ ˙ 0 0 1 0 θ ?? X δ e X δ p Z δ e Z δ p M δ e M δ p ? ? ? ? δ e δ t ? ? ? ? + 0 0 0?X u ?X w ?? ? ? ? ? ? ? ? ? ? ? ? u g 0?Z u ?Z w ?M u ?M w ?M q ? + w g q g 0 0 0 ? B u u + ? Ex˙ = Ax + ? B w w? ? Multiply through by E ?1 to get new state space model x˙ = Ax + B u u + B w 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 v g , p g , and r g . ? Can now compute the response to speci?c types of gusts, such as a step or sinusoidal function, but usually are far more interested in the response to a stochastic gust ?eld ? ? ? ? ? ? ? ? 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 0 T →∞ 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 di?erence from the Random Number block is that this block produces output at a speci?c sample rate, which is related to the correlation time of the noise. ? Continuous white noise has a correlation time of 0 ? ?at power spectral density (PSD), and a covariance of in?nity. – 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 e?ect 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 t c much smaller than the fastest dynamics f max of the system. 1 2π t c ≈ 100 f max ? Fall 2004 16.333 12–7 ? Power spectral densities for Von Karman model (see MIL-F-8785C) ?5/6 Φ u g (Ω) = σ 2 2L u ? 1 + (1.339L u Ω) 2 u π where σ u is the intensity measure of the disturbance, Ω is the spatial frequency variable, and L u is a length scale of the disturbance. – All scale parameters depend on day, altitude, and turbulence type. – Di?erent models available for each direction. – Scale length parameter L (in feet) varies with altitude – MIL-F- 8785C model valid up to 1000 feet h L u = (0.177 + 0.000823h) 1.2 – Turbulence intensity for low altitude ?ight MIL-F-8785C as 0.1W 20 σ u = (0.177 + 0.000823h) 0.4 3 where W 20 is the wind speed as measured at 20 ft 3 W 20 < 15 knots is classi?ed as “light” turbulence W 20 ≈ 30 knots is “moderate” W 20 > 45 knots is “heavy” Fall 2004 16.333 12–8 – Standard approach: assume turbulence ?eld ?xed (frozen) in space ? aircraft is just ?ying into it. 3 Then can relate turbulence spatial frequency Ω to the temporal frequency ω that the aircraft would feel ω ω = ΩU 0 Ω ≡ U 0 ? – Used to model how the aircraft will behave in bumpy air. profiles 10 ft/sec altimeter 10 ft/sec Longitudinal gust Atmospheric gust Touchdown 10 f t/sec 300 ft on radio 5 sec Vertical gust Lateral gust Figure 2: Wind gust examples ? ? Fall 2004 16.333 12–9 Wind Shear ? Wind shear is of special interest because it involves signi?cant 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 u g = Δ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. Glide slope Flight path Flare Mean wind Figure 3: Aircraft descending into a horizontal windshear But can write the changes in altitude (use h ˙ = Δh/Δt) as ? h ˙ = U 0 (θ ? α) ? ? ? 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 = C h x with controls ?xed (zero) and input u g x˙ = Ax + ? ? B w (:, 1) ? B w (:, 1)u g = Ax + ? du h dh du ? ? = Ax + B w (:, 1)C h x ? dh du ? = A ? + B w (:, 1)C h x dh ? dynamics have been modi?ed because of the coupling between the change in altitude and the change in forces (with u g ). ? Closed this loop on the B747 dynamics to obtain du Phugoid Poles dh 0 0.08 0.15 0.2 -0.0033 ± 0.0672i -0.0014 ± 0.1150i 0.0002 ± 0.1442i 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 % Gust modeling 2 % 16.333, Fall 2004 3 % Jonathan P. How 4 5 Xu=-1.982e3;Xw=4.025e3;Zu=-2.595e4;Zw=-9.030e4;Zq=-4.524e5;Zwd=1.909e3; 6 Mu=1.593e4;Mw=-1.563e5;Mq=-1.521e7;Mwd=-1.702e4; 7 % 8 g=9.81;theta0=0;S=511;cbar=8.324; 9 U0=235.9;Iyy=.449e8;m=2.83176e6/g;cbar=8.324;rho=0.3045; 10 Xdp=.3*m*g;Zdp=0;Mdp=0; 11 Xde=-3.818e-6*(1/2*rho*U0^2*S);Zde=-0.3648*(1/2*rho*U0^2*S); 12 Mde=-1.444*(1/2*rho*U0^2*S*cbar);; 13 % 14 Ehat=[m 0 0 0;0 m-Zwd 0 0;0 -Mwd Iyy 0;0 0 0 1]; 15 Ahat=[Xu Xw 0 -m*g*cos(theta0);[Zu Zw Zq+m*U0 -m*g*sin(theta0)]; 16 [Mu Mw Mq 0];[ 0 0 1 0]]; 17 Bhat=[Xde Xdp;Zde Zdp;Mde Mdp;0 0]; 18 % 19 % form the gust input matrix 20 Bwhat=[-Xu -Xw 0;-Zu -Zw 0;-Mu -Mw -Mq; 0 0 0]; 21 % 22 % add height state 23 Ehat(5,5)=1; % \dot h state not coupled 24 Ahat(5,5)=0;Ahat(5,[1:4])=[0 -1 0 U0]; % add height state 25 Bhat(5,2)=0; % noinput 26 Bwhat(5,3)=0; % no input 27 28 % form the full model 29 % E \dot x = \hat A x + \hat B u_controls + 30 % ==> \dot x = A x + B u_controls + B_w w 31 % 32 A=inv(Ehat)*Ahat; 33 B=inv(Ehat)*Bhat; 34 Bw=inv(Ehat)*Bwhat; 35 % 36 % set u_controls=0 37 % assume that w_g=q_g=0 and 38 % u_g = (du/dh) h 39 % 40 du_dh=[0 .08 .15 .2]; 41 A1=A+Bw(:,1)*[0 0 0 0 1]*du_dh(1); 42 A2=A+Bw(:,1)*[0 0 0 0 1]*du_dh(2); 43 A3=A+Bw(:,1)*[0 0 0 0 1]*du_dh(3); 44 A4=A+Bw(:,1)*[0 0 0 0 1]*du_dh(4); 45 \hat B_w w 46 ev1=eig(A1);ev2=eig(A2);ev3=eig(A3);ev4=eig(A4); 47 plot([ev1 ev2 ev3 ev4],’x’) 48 [ev1 ev2 ev3 ev4]