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]