1

Phugoid Mode Dynamics Excited by Wind (Updated 11/20/15)

1.Response of the Phugoid Mode to Wind

The 2D state equation for the phugoid mode is:

.(1.1)

The transfer function for the command input is obtained from:

. Hence, .(1.2)

Now, it is easy to show that: .

Hence, the explicit expression for is:

(1.3)

where .

In a similar manner, for the wind gust input we have:

.(1.4)

2. Phugoid Response of the Boeing 747(@ M=.9) to a Wind Gust

The dimensional stability derivatives for this plane are:

Hence, (1.4) becomes:

(2.1)

where .

[Note: The above transfer function matrix can be obtained directly from the (A,B,C,D) state space system by using the Matlab command ‘ss2tf’.]

The u(t) FRF for input ug(t)- From (2.1) we have

.(2.2)

To obtain explicit expressions for the FRF magnitude, , and phase :

(i)The numerator of (2.2) is:

(ii)The denominator is:

Hence, (2.2) becomes:

.(2.3)

The θ(t) FRF for input ug(t)-

In exactly the same manner as above, we obtain:

.(2.4)

A note in relation to decibels(dB): For a magnitude, M, the conversion to dB is: . For a power, P, the conversion to dB is: . Because power and magnitude are typically related via , the two dB expressions are equivalent.

We can now use the Matlab ‘bode’ command to plot the FRFs (2.3) and (2.4)

Figure 1. Plot of the phugoid mode FRF’s and .

Response to a ‘sharp edge’ gust: By the term ‘sharp edge gust’, the authors mean a gust having the form of a step. We will now use Matlab to compute the responsesto , which is a unit step. However, first let’s try to figure out what we should expect to get.

From (2.1) we have: . (2.5)

The static gain of (2.5) is . We also have:

, or ; ;

Hence, the response, should converge to a value of 1.0 as . Moreover, this response will be essentially achived at . The damped natural frequency is , and so the damped natural period is . And so we should expect ~10 periods of oscillation prior to reaching the steady state condition. To obtain the unit step response we use the Matlab command ‘step(Hu)’:

Figure 2. The unit-step response of the system .

The response plotted in Figure 2 is accurately predicted by the above analysis.

In exactly the same mannar as above, we have the unit-step response of the system shown below.

Figure 3The unit-step response of the system .

QUESTION: What was the frequency content that was included in the unit-step input?

ANSWER: The Laplace transform of the unit step input is . The associated Fourier transform is . The Fourier transform gives the freqquecy-by-frequency description of the input . The Bode plot below shows this description.

Figure 3 The Fourier transform of the input

The Phugoid Mode Responses to Von Karmon Turbulence

Suppose that the plane encounters head wind that can be modeled using the approximate von Karman wind power spectral density (psd) model: where we have defined and . The function is called a shaping filter.

Even though is the distribution of the wind energy as a function of frequency, this energy has the same shape as the FRF of the fictitious system . The bandwidth of the approximate von Karman longitudinal wind psd model is: .

2.1 The Boeing 747 speed at M=0.9 and 40,000 ft. is . The natural frequency of the phugoid mode, , is so very low that it will almost certainly be within the turbulence bandwidth. To verify this, find the value of the spatial correlation length parameter such that the von Karmon psd bandwidth, , coincides with the phugoid resonance frequency.

Solution:

For a spatial correlation length greater than this, the phugoid mode will be less excited.

2.2 For a realistic spatial correlation length of 30[1] miles, find the turbulence bandwidth.

Solution: .

Since this frequency is far less than , the phugoid resonance frequency will be minimally excited by the lower high frequency energy in the turbulence.

For the turbulence, we have (assuming for convenience that ):

.

We will now simulate a von Karmon turbulence time series by running white noise through the turbulence shaping filter.

Figure 3. Plots of fictitious white noise input to the shaping filter, and of .

Figure 4. Plots of the phugoid mode and responses to turbulence .

Finally, from the data in Figure 4, we will use the Matlab ‘pwelch’ command to estimate the psd for .

Figure 5. Data-based estimate of the psd for. Note that the phugoid mode natural frequency in Hz is .

QUESTION: How would you compute the theoretical psd for ?

MATLAB CODE

%PROGRAM NAME: phResponse2Wind.m

%Simulate Boeing 747 response to von Karmon wind turbulence

% Plane Information:

u0=871; g=32.17; Xu=-.0076; Xw=-.00128; Zu=-.1053; Zw=-.3944;

A=[Xu,-g;-Zu/u0,0];

B=[-Xu;Zu/u0];

C=eye(2); D=zeros(2,1);

Hss=ss(A,B,C,D);

[Hn,Hd]=ss2tf(A,B,C,D,1);

Hu=tf(Hn(1,:),Hd);

Hth=tf(Hn(2,:),Hd);

wd=imag(roots(Hd)); wd=abs(wd(1,1))

figure(10)

wvec=logspace(-2,0,1000);

subplot(2,1,1),bode(Hu,wvec)

title('Phugoid FRF for u(t) in Response to ug(t)')

grid

subplot(2,1,2),bode(Hth,wvec)

title('Phugoid FRF for th(t) in Response to ug(t)')

grid

pause

figure(11)

subplot(2,1,1),step(Hu)

title('Phugoid u(t) in Response to unit step u_g(t)')

grid

subplot(2,1,2),step(Hth)

title('Phugoid th(t) in Response to unit step u_g(t)')

grid

pause

%------

% Turbulence Information:

Lu=30*5280; % Spatial Correlation Length of 30 miles

sigu2=pi/(2*Lu);

gs=(sigu2*2*Lu)/pi;

tau=1.339*Lu/u0;

wBW=1/tau;

Hug=tf(gs/tau,[1,wBW]);

%======

%TURBULENCE SIMULATION:

n=10000;

uwn=sqrt(sigu2)*randn(n,1);

Au=-wBW; Bu=1; Cu=1; Du=0;

dt=tau/1000;

t=0:n-1; t=dt*t';

sysu=ss(Au,Bu,Cu,Du);

ug=lsim(sysu,uwn,t);

figure(12)

subplot(2,1,1), plot(t,uwn)

title('White noise input to the shaping filter')

grid

subplot(2,1,2),plot(t,ug)

title('von Karmon turbulence')

grid

pause

%======

%SIMULATION OF PHUGOID RESPONSE TO TURBULENCE:

uth=lsim(Hss,ug,t);

u=uth(:,1);

th=uth(:,2);

figure(13)

subplot(2,1,1),plot(t,u)

title('u(t) Response to u_g(t)')

grid

subplot(2,1,2),plot(t,th)

title('th(t) Response to u_g(t)')

grid

pause

figure(14)

fd=wd/(2*pi)

fs=1/dt;

pwelch(u,[],[],[],fs);