NASA Technical Memorandum 110437

Finite Element Aircraft

Simulation of Turbulence

R. E. McFarland, Ames Research Center, Moffett Field, California

February 1997

National Aeronautics and

Space Administration

Ames Research Center

Moffett Field, California 94035-1000

NOMENCLATURE

b

f(s),g(s),h(s)

F(z),G(z),H(z)

V

L., L,,, &,,

l(s),m(s),n(s)

L(z),M(z),g(z)

T

dp

dq

dr

X,Y

x,y

0i

P

Pa, qa, ra

Pc, qc, rc

The span length

Dryden u, v, and w turbulence filters in Laplace form

Dryden u, v, and w turbulence filters as z-transforms

Limited horizontal velocity

Turbulence characteristic lengths

Conventional model p, q, and r filters in Laplace form

Conventional model p, q, and r filters as z-transforms

The cycle time

Distance between right and left wing centers of pressure along span

Distance from fuselage center of pressure to the center of pressure of the

horizontal tail

Distance from fuselage center of pressure to the center of pressure of the

vertical tail

Independent random number sequences with Gaussian distribution for use in

the correlation algorithm (X = 03, Y = 04)

Dependent number sequences with Gaussian distribution from correlation

algorithm

Uncorrelated X and Y filtered by H(z)

Right and left wing center of pressure vertical turbulence velocities as

developed by the correlation algorithm, x and Y filtered by H(z)

Unity-variance, zero-mean independent random noise sequences with

Gaussian distribution

The spanwise vertical correlation coefficient

Rotational sums of stochastic turbulence and discrete gusts

Conventional model rotational stochastic turbulence components

iii
Ps, qs, rs

UM, VM, WM

LIa , V a _ Wa

I_ic , Vc _ W c

Us, Vs, Ws

Ug, Vg, Wg

pg, qg, rg

(Yu, (Yv, (Yw

FEAST model rotational stochastic turbulence components

Fuselage center of pressure u, v, and w translational velocities

Translational sums of stochastic turbulence and discrete gusts

Conventional model translational stochastic turbulence components

FEAST model translational stochastic turbulence components

Translational discrete gust components

Rotational discrete gust components

Turbulence dispersions

iv
FINITE ELEMENT AIRCRAFT SIMULATION OF TURBULENCE

R. E. McFarland

Ames Research Center

SUMMARY

A turbulence model has been developed for realtime aircraft simulation that accommodates

stochastic turbulence and distributed discrete gusts as a function of the terrain. This model is

applicable to conventional aircraft, V/STOL aircraft, and disc rotor model helicopter simulations.

Vehicle angular activity in response to turbulence is computed from geometrical and temporal

relationships rather than by using the conventional continuum approximations that assume uniform

gust immersion and low frequency responses. By using techniques similar to those recently

developed for blade-element rotor models, the angular-rate filters of conventional turbulence models

are not required.

The model produces rotational rates as well as air mass translational velocities in response to both

stochastic and deterministic disturbances, where the discrete gusts and turbulence magnitudes may

be correlated with significant terrain features or ship models.

Assuming isotropy, a two-dimensional vertical turbulence field is created. A novel Gaussian

interpolation technique is used to distribute vertical turbulence on the wing span or lateral rotor disc,

and this distribution is used to compute roll responses. Air mass velocities are applied at significant

centers of pressure in the computation of the aircraft's pitch and roll responses.

DISCUSSION

The Simulation of Rotor Blade Element Turbulence (SORBET) model was developed (ref. 1) in

order to accommodate stochastic and discrete gust velocity components at elements along the rotor

blades of a helicopter simulation. This model received improved pilot comments over an

implementation of the conventional MIL SPEC body-fixed formulation (ref. 2) for low-speed and

low-altitude flight. The introduction of air mass velocities at the blade elements obviates the angularrate filters, because the moment arms to the aerodynamic points of application are implemented.

This feature has the added benefit that gusts may be input as a function of terrain, so that angular as

well as translational disturbances may then be correlated with the terrain profile. This is quite

important for the simulation of low speed flight of various vehicles in the vicinity of buildings,

ships, etc.

Recent research has shown that air wake and turbulence modeling require improvement. "A key

issue with these turbulence models is the way in which they are applied to the aircraft. In all cases
thiswasthroughtheaircraftcenterof gravity.Morecomplexmodelswouldrequireamoreeffective

interfacewith theaircraft"(ref. 3). Insofarascommunicatingterrain-relatedinformationto the

aircraftmodel,two optionsareavailable."Whenairflow gradientsarepresentthataresignificant

overdistancesthataresimilar to themaindimensionsof theaircraft,thenit becomesnecessaryto

considerintroducingforcesandmomentsdueto gradientseitherby (1)defininggradientsatthe

singlepoint andusingappropriatederivativesof forcesandmoments,or (2)calculatinglocal forces

andmomentsondifferentpartsof theaircraftandintegratingtheseto obtaintotal forces"(ref.4).

In this material,atmosphericdisturbancesareseparatedintostochasticturbulenceanddiscretegusts.

Stochasticturbulenceis treatedin somedetail,anddifferencesfrom theconventionalformulation

arenotedin angularresponses.Turbulenceisdevelopedfrom theconventionalatmospheric

relationships(ref.2), exceptthatrotationaltermsaredevelopedby a distributionof translational

velocitiesoverthespecificvehicle'seffectiveaerodynamicsurfacesusingtransportdelays;anovel

Gaussiandistributionalgorithmis usedin thedevelopmentof roll rateasa functionof thetwodimensionalvertical velocityfield. A comparisonis madeto theconventionalsystem.

Discretegustsandturbulencemagnitudesmaybedefinedfrom acomputationalfluid dynamics

analysisof flow overaterrainprofile. It is assumedherethatsuchadatabaseis available;the

remainingproblemis to applythisdatabaseto anaircraftsimulationmodelin realtime, thereby

correlatingaircraftresponsesto turbulenceasafunctionof theterrain.

Fora distributedaerodynamicmodeltheappropriaterotationalbehaviorof thevehiclearises

naturallyfrom forcedifferentialsof thevariousaircraftcomponents,providedonly thatthe

translationalatmosphericdisturbancesareproperlydistributedto theappropriatecentersof pressure

(CP)of thesecomponents.Suchwasthecasein theSORBETmodel,whichusedarotatingbladeelementmodelof theUH60 helicopter,with anequal-annulidistributionof aerodynamiccenters

alongtheblades.However,for total forceandmomentmodels,typicallyusingaerodynamic

coefficients,thevehicle'srotationalbehaviorin responseto discretegustsrequiresgeometrical

considerations,which arereviewedherein.

Therealtimetransferof databaseinformationto thesimulationmodel,generallyincludingvector

fieldsof translationaldiscretegustvelocitiesandturbulencedispersions,mayalsoinclude

atmosphericgradientdata.This impliesasingle-pointtransfer(probablythevehiclecenterof

gravity) for eachcomputecycle.Thisformulationappearsto becomputationallyefficient.However,

atmosphericgradientsderivedfrom agrid of translationalvelocitydifferencesareindependentof a

specificvehicle'sgeometryandaerodynamics.In thedevelopmentof vehiclerotationalratesthis

mustbeconsidered(seethediscussion,ref.5). Furthermore,usinggradients,theimplied

interpolationbetweendatapointsis likely to smoothsignificantdisturbanceparametersunlessavery

large(dense)databaseis available.Spacingassmallasoneor two metershavebeensuggestedin

regionsof [flow field] instability (ref.4). Considerthatin ahelicoptersimulationof landingona

ship,"changesin meanairflow in themodelappearedtoogradual,particularlyaroundtheedgesof

thehanger"(ref.4)."It shouldbementionedherethattheaerodynamiceffectsof gustgradientsare

not quitethesameastheeffectsof theairplaneangularvelocities"(ref.5).

In this paperasimilar techniqueto thatdevelopedin reference1is usedto distributeturbulence

componentsto thecentersof pressureof aconventionalaircraft.A modelis createdcalledFinite

2
ElementAircraft Simulationof Turbulence(FEAST).With minor modificationthesetechniques

maybeappliedto helicoptersimulationsthatusearotordiscmodel.

Thecompellingreasonto usea modellike this isto accommodateacomponentlevelmodelof the

airframeaerodynamics.Thedistributionof translationalturbulencecomponentsaswell asdiscrete

guststo thevariousaerodynamiccenterswouldthenof itself inducethecorrectrotationalbehavior

of thevehicle.For acoefficientmodelthatdoesnotdistributetheaerodynamics,however,wemay

approximatetherotationalbehaviorof airframesby consideringonly thevehicle-specificmoment

armsto theeffectiveaerodynamiclocations.Generatingrotationaltermsfor arbitraryaircraft is then

reducedto a geometryproblem.

Thecorrectmomentarmsto usein sucha formulationwouldseemto bethosethatreplicatethe

outputsof theconventionalturbulencemodel(ref.2) in standardflight regimes.It will beseen,

however,thatforcing suchacorrespondencedoesnotproducerealisticmomentarmsfor usein the

FEASTmodel,whichis developedfrom elementarymathematicalrelationships.Rather,whenthe

correctseparationsof centersof pressureareutilizedin theFEASTmodel,therotationaltreatment

in MIL SPECis broughtintoquestion,in particularforthehigherfrequencyregions.This is

examinedin detail.

In theFEASTmodelit ismaintainedthat"this idealizationof theturbulencefield permits

idealizationof theairframe,for calculationof the(generalized)aerodynamicforcesandmoments,as

atwo-dimensionalifting surface(in thex-y plane)fortheverticalcomponentWg and as onedimensional force distributions (along the x axis) for the horizontal components Ug and Vg" (ref. 5).

An exponential correlation form is assumed between vertical turbulence velocities. As is shown, this

form produces considerably more roll activity at low altitudes than does the conventional model. The

exponential correlation form is certainly not optimal, but improvements in this form are beyond the

scope of this paper. Referring to scale lengths of a couple of hundred feet, Etkin stated that at "such

small scales, the variation in gust velocity over the airplane becomes important, and analytical

methods of some refinement and complexity are indicated" (ref. 6).

The techniques developed in this paper for stochastic turbulence are applicable to either a gradient or

distributed finite element model. Terrain related gusts are treated using the specific vehicle

geometry, and such a formulation should collapse to the gradient method in the mathematical limit,

providing that the specific vehicle geometry is considered.

Geometry and Transport Delays

The FEAST model uses the MIL SPEC translational filters for stochastic turbulence, but does not

use the angular filters (ref. 2). These filters are replaced by using the geometrical separation of

centers of pressure, and a novel spanwise correlation factor, as developed from a two-dimensional

vertical turbulence field. In obviating the pitch and yaw angular filters, the model uses transport

delays for the emersion rate into the turbulent field.

The use of a two-dimensional vertical turbulence field permits the development of roll activity as a

function of the spanwise distribution of velocities. "...when structural modes are significant, the

exact distribution of the turbulence velocity components over the airframe should be considered ...
For purposesof MIL-F-8785B, it is acceptable to consider ug and vg as being one-dimensional

functions of x, but Wg must be considered two-dimensional, a function of both x and y, for

evaluation of the aerodynamic forces and moments ... The spanwise (y) variation of wg often

produces significant lateral responses of the airplane and is, therefore, always to be considered,

except possibly for very slender configurations when approved by the procuring activity." (Ref. 5.)

Similar to the techniques used in SORBET, a turbulence "onset line" is defined perpendicular to the

aircraft body axis in the lateral plane, along which the outputs of the translational turbulence filters

are computed. For convenience, this line is placed at the fuselage center of pressure. Stochastic

turbulence components are developed in the body frame, and the onset line is perpendicular to the

longitudinal body axis. This line is useful in considering the propagation of stochastic turbulence to

the horizontal and vertical tail centers of pressure as a function of the vehicle velocity. Discrete

terrain related gusts are later shown to be treated differently from stochastic contributions because

they are instantaneously applicable at distributed locations.

Turbulence velocities u M and vM are developed (M=middle line) at the fuselage CP. The vertical

velocity at the fuselage CP, w M, is produced from two uncorrelated Gaussian inputs given by X and

Y, as will be shown. Additionally, from these uncorrelated inputs, correlated Wr and we at the right

and left wing centers of pressure (along the span) are developed. Only four independent noise

sources are required, as in the conventional MIL SPEC (ref. 2) formulation, and a distribution occurs

to the points of application to account for both the geometry and transport lags.

As the vehicle travels through a frozen field of turbulence, the history of the turbulence encountered

at a downstream point is equivalent to the history at the fuselage CP, with a time of application

equivalent to a transport delay. This delay is inversely proportional to the horizontal aerodynamic

velocity 1 of the aircraft V, and is given by tq = dq/V for the vertical field that causes pitch rate, and

by tr = dr/V for the horizontal field that causes yaw rate. For example, the lateral turbulence

velocity applied to the vertical tail's center of pressure is given by VM(tr). This is merely a

restatement of Taylor's frozen field hypothesis, which states that temporal and spatial relationships

are proportional.

The yaw rate (due to the tail's side velocity) is developed by the difference vM(tr) -- VM(O) divided

by the fuselage-to-vertical-tail separation d r. The sign convention results from the fact that the

translational turbulence and gust velocities are air mass velocities. The pertinent continuum

relationships are developed in an early MIL SPEC document (ref. 5) which produced the partial

derivatives:

pg = -- OqWg_

qg =

"V" is actually the torward velocity of the aircraft, with a lower limit discussed in appendix 1.

4
If the velocity of the air mass is positive downward on the right wing (and absent at the left wing),

then the rotationally invariant vehicle appears to roll left. If the air mass velocity at the C.G. is

positive downward (and it is absent at the horizontal tail) then the rotationally invariant vehicle

appears to be pitching up. If the velocity of the air mass at the C.G. is positive rightward (and it is

absent at the vertical tail) then the rotationally invariant vehicle appears to be yawing left. The

translational turbulence and gust velocities of the air mass are subtracted from the vehicle velocities

in simulation, whereas these rotational terms are added to certain aerodynamic terms (ref. 5).

The rotational terms may be applied to a C.G.-based aerodynamic model. For a distributed

aerodynamic model the translational velocities at various aerodynamic centers will induce the proper

rotations by the development of force differentials.

For the computation of roll rate, two correlated vertical velocities are produced along the right and

left wings, at their centers of pressure. The distance dp is the lateral separation of these two points.

The correlated vertical velocities at the wing CPs are developed from two uncorrelated random

sequences. As will be shown, these uncorrelated vertical velocities are also used for the computation

of the fuselage CP vertical velocity WM.

The generalization of this geometry to a helicopter model requires different points of application.

These points could include four locations about the hub, typically at radii of about seventy-five

percent of the rotor radius (ref. 7). In this fashion the rotational gradients due to turbulence on the

rotor disc may be computed, as well as the influence on the C.G. and tail rotor.

The Conventional Stochastic Model

The conventional (MIL SPEC) turbulence model is presented in block diagram form in figure 1.

n,(s) _ uc(s) ._

v(s)

n3(s) _ we(s)

n,(s) v

qc(s)

p¢(s)

Figure 1. Conventional model.
The sums and differences of random variables in this figure are used for convenience. The basis for

these operations is given in the next section.

The Dryden form for the longitudinal, lateral, and vertical turbulence filters (that produce Uc, Vc, Wc)

are given in Laplace form as follows (refs. 1 and 2):

s+ V/l 

g(s) = 3V_--]--_[s + V/ff-3_,]

[s+ Wz ]2

h(s)

+

[s + 2

and the roll, pitch, and yaw rotational filters (that produce Pc, qc, rc) are given by

t(s) (E/4b)7/6 0_.8V

/-_3(s + 7cV/4b)

m(s) - rcs/4b

s + tcV/4b

n(s) - -Irs/3b

s + zrV/3b

where the quantity b is the wing span. The conventional model's rotational filters have a history of

being brought into question. For example, the quantity b has even been replaced by the total aircraft

length in the pitch and yaw filters (ref. 8).

In the conventional model, it should be noticed that roll rate is independent of vertical activity at all

flight conditions (one-dimensional field). The pitch and yaw rates consist of differentiated vertical

and lateral velocities, using first-order filters with gains that are the inverse of the aerodynamic

velocity. The z-transform solutions to the above Laplace filters are given in appendix 2.

The translational turbulence velocities are the same in both the conventional and FEAST

formulations, where the point of application is the fuselage center of pressure. The rotational

formulations are different. In particular, FEAST's roll rate uses certain properties of random

variables, as examined in the following section. Two correlated vertical velocity sequences are

developed at the wing centers of pressure, and an exponential form is assumed for the correlation

between the velocities. The moment arm between these aerodynamic centers is then used to develop
roll rate. Pitch and yaw rates in the FEAST model are developed from transport delays acting on the

vertical and lateral velocity histories.

Figure 2 summarizes all of the relationships in z-transform notation for the conventional model

given above, where its outputs have the subscript c. In addition, this figure introduces the FEAST

model, where its outputs have the subscript s.

Appendix 1 is provided to show an implementation of the transport delays that are indicated in the

FEAST paths of this figure.

Ol(Z)

r/2(z)

u,.(z)= u,(z)

r

re(z)

r.,.(z)

v,.(z) = v.,,(z)

03(0

04(z)

, z ,Tt

w,(z) = w,(z)

qc(Z)

q,( z)

pc(z)

v

v

2.4__ p) ] p,.(z) ,,_

dp ]

Figure 2. The conventional model and the FEAST model.

Correlating Variables

In order to create two correlated sequences, two independent, zero mean Gaussian sequences {X}

and {Y} are first developed, both with the same standard deviation, O'x = cry = cr. These may be

created in realtime simulation by using software 2 such as XNORM ©. Two correlated sequences {x}

and {y} with the same standard deviations o"x = O'y = O" (and all other moments) may then be

created as follows. The correlation coefficient for the sequences {x} and {y} is defined by

O'xy (0</9<1)

O'xO'y

2The copyrighted Ames Research Center standard simulation software is available to industry.

7
whereonly directrelationshipsarepermitted.This is accomplishedby excludinginverse

relationshipsin theselectedcorrelationcoefficient,andthelowerlimit of zeroisproperratherthan

negativeunity (ref. 9).

Giventhecorrelationcoefficient p as shown above, two symmetrically correlated sequences {x,y}

may be computed from two uncorrelated sequences {X, Y} with identical statistical properties by the

functions:

and their correlation coefficient is p.

The proof is as follows: E(XY)=O for independent variables, and E(X2) = E(y2)= cr 2. The

expected values for the dependent variables are given by:

=_E{[I_+_/x+/,_- _t_12)

--1⁄4E{t ,,_+_/2x2+l,,_-_t2Y2}

= l_2[(l+p)+(l_p)+(l+p)+(l_p) ] = _2

E{y2} = 1E { [( _1_/ - I_- p )X + ( I_---P + _/-_ ) Y ]2 }

= 1⁄4o2[(l+p)+(1-p)+(l+p)+(1-pt]= o2

E{xy} = 1E{[(_+ ]._-p-p)X+(I_/_-I._Z-p)Y][(1._-_-I_J_-p)X+( 1..v/_+ 1.,/_-p)Y]}

: 1⁄4E{(_p+_,/___)( 1/1/__1/_-p)_+(,_,/_-,_;_-_)( ,,/_ +_),,=}

1 2

Thus

E{xy} o-_,, (0___p _<l)

p - 0" 2 O-xO'y

8
In the FEAST model this result is used to produce correlated turbulence velocities at the wing's

spanwise centers of pressure by relating 3 the vertical velocities x _ Wr and y _ we as linearly

filtered random variables. These are computed from the uncorrelated X _ wR and Y _ WL as

shown below. An exponential correlation form is used:

E_w(ro)w(r 0 + r)] = crZe -_

where for this model L is the characteristic length 4 of the vertical component of turbulence and r is

the distance between the two points of interest (the difference between the wing centers of pressure,

dp). Thus, the correlation coefficient becomes

p = e- dp//fLw

and the criterion for only positive correlation is satisfied due to the lower limit of the MIL SPEC

definition (ref. 2) for low-altitude turbulence models

'%12

10 h<lOfi

= h lOft<h<lOOOft

1000 h > 1000 ft

where h is the vehicle altitude. For a center of pressure difference dp = 32.17 ft, the correlation coefficient would be as shown in figure 3.

1.0 i ..

0.8 ...................... ',....................... i......................... _ ................... _.......................

o.G....................... i.......................... i........................... i............................ i..........................

o.4 ........................ i..................................................... i........................................ 5

=_o.2 ....................................................... ............

o

u 0.0

0 200 400 600 800 1000

Altitude (or Characteristic Length) - ft

Figure 3. Wing CP correlation coefficient.

3The random sequences first pass through the H(z) atmospheric filter. Gaussian variables that pass through linear filters

remain Gaussian variables.

The characteristic length is the spatial separation of two points where the correlation of velocity components is

assumed to vanish.

9
Theverticalturbulencevelocity at thefuselageCPmayalsobecomputedfrom theuncorrelated

variables.This quantity w M results from the correlation algorithm by setting p = 1, a condition

called for when the distance between the two points of interest vanishes (r = 0). Using the uncorrelated quantities X and Y we see that this produces the linear combination for the centerline vertical

turbulence velocity

Y+ X] wM : WL+WR]

which is the Gaussian centroid, and all statistical moments are retained, as in the case of the

turbulence at the wing span positions. X and Y are filtered to produce Wg and wL.

To complete the picture for characteristic lengths, the MIL SPEC restricted lower limits for the

horizontal components are also used, and it is noted (ref. 2) that the off-axis dispersions are determined once the vertical component (y,,, is selected.

h

(0.177 + 0.000823h) 12

(10 ft <_h <_1000 ft)

O't/ _ O" V

0" W

(0.177 + 0.000823h) °4

(10 J_ _<h _ 1000 ft)

Using the relationships of this section, the spanwise center-of-pressure vertical velocities Wr(Z) and

we(z) for the right and left wings, and the midline vertical velocity WM(Z), are created from two

independent noise sources. An example time history of these velocities is presented in figure 4,

where the correlation coefficient used was t9 = 0.6525.

j:\,, ^_ / _ ,'_.

.... " ' ¢ - . .- _" /me

d / -_"'":'":"":"""'""' - ,ry.,I- "w..J v¢ _!i!iii"_'"_" ":"":""''"""''::"''"" w _/_": 5":.

Figure 4. Vertical turbulence at wing and fuselage CPs.

The velocities in the above figure were created at an altitude of 40 feet, an aerodynamic velocity

of 5 ft/sec, and d e = 17.08 ft. The CP vertical velocity histories were used to produce the roll rate

given in figure 5.

10
o 0.05 -

0.00-

-0.05 -

rr -0.10 -

O

rc -0.15--

0

I I

5 10 15 20

Time - sec

Figure 5. Resultant roll rate history.

Note that correlation between the vertical velocities at the wing CP locations is only a statistical

measure. It does not provide any indication of the direction of roll with respect to the fuselage

vertical velocity. Indeed, roll is completely uncorrelated with the fuselage vertical velocity, as the

following expected value shows:

E{psWM}_ E{[w,,-wR][w,,+wR]} = o_-o_ = o

Vertical Velocity

Consider the vertical velocities w L and w R, which are uncorrelated Gaussian random variables with

the same dispersion (7, and zero mean value. Such variables with identical statistics have the

following properties:

E(w,w_} = 0

_{Iw_+w.l_}: _{w_+w_}: _

_{iw,.-w_l _}= _{w_+w;}: _o_

e{[wL+w,,IwL-wR]}= e{w?.- w_} = 0

The expected values of both the sums and differences produce new random variables that have twice

the variance of the originals. Hence, we may create new random variables with the same dispersion

cy by dividing by the square root of two. That is, in terms of discrete, unity-variance random

variables

ns(z) FI4(Z)-1"O3(Z )

4_

04(Z)-- O3(Z)

T/6(Z ) :

11
and from the above properties the sum represented by r/5(z ), the Gaussian centroid, is uncorrelated

with the difference represented by r/6(z ). For practical purposes the random turbulence vertical

velocity may be created in z-transform notation from the sum

w(Z) = cyw_H(z)[r]4(z)_2-_3(z)l = cr,.,._H(z)rl5(z )

where the system function for vertical velocity in Laplace notation is

w(,t= o_,,h(sl-- Ow_ i7;_

The autospectrum or power spectral density (one-sided PSD) of the Dryden form of the vertical

response is given by (ref. 2)

3VO_w(O,2 +v_/31_)

,_(o,_+_/,z)_

"The total area under the autospectral density function is the variance of the data plus the square of

the mean value of the data" (ref. 9). Since the mean value of the random input is zero, the variance

of this autospectrum is given by the frequency integral from zero to infinity.

- 3V°'2'[ n" 2(V//_)2 _ ]

n'L,v 2(V//_,,) 3 4(V/I__)3

v_/z_ ]_

-_ 2/aO) (o,2+v_/,z)j

= O_w

Rotational Variances

In appendix 3, the conventional formulation's rotational relationships are reviewed and the rotational

variances are shown to be

2 - \_

0.4awrC ( rc _

12
°q c O-2-2(3fCLw + 2/

32b2( 4--_ + 1) 2

=

1862I +l)2

Also, in appendix 4 the FEAST model's rotational relationships are developed. The rotational

variances are shown to be

0-2s _ 20"2(1-t9)

wE/ a/- q1

0-2

0-gs- dq2 -2 e '+2

0-2s- 0-2[ (df//Lv-2)e-dr//Lv-_r +2]

Equating Variances

Equating the variances between the conventional and FEAST models would seem to be an excellent

technique for quality control. In pitch and yaw the sharper responses of the FEAST model could be

used to justify the lower-frequency continuum responses of the conventional model, and in roll the

influence of the vertical correlation relationships could be examined to determine the altitude where

behavior was similar.

Unfortunately, as determined in this section, the technique of equating the variances produces

completely unrealistic CP separation distances.

For comparison of the two formulations with equal variances, the distance between the wing centers

of pressure for use in the FEAST model may be determined if the span (b) and the characteristic

length (/__,) are known. Solving for dp in the roll dimension, this requires the solution to the

following transcendental equation, obtained by equating o'2 c = 0-p2s.

1 -- e -d/_v

0.2rc_ --_-) -

13
Similarly, equatingo'gc = O'qs and o",%= O'2s also produces transcendental equations

3 + 2(4b/rd-,w)

2(4b/M_,,)[1 + (4b//d_,_)]2

3 + 2(3b/:rc/_,)

2(3b/a'L_,)[1 + (3b/;rd_v)]2

, 2

2 +[(dr/L_,)- 2]e -'4r/L*

( d r ,/L_, )2

which may be used to solve for the other two CP differences. These three equations are independent

of both variances and vehicle velocity. They are easily solved using Newton's iteration (ref. 10).

Very few cycles are required.

For b=32.17 ft and Lw=250 ft (Lu=791.48 ft) Newton's iteration produces the values given in

figure 6 for the center-of-pressure separation distances dp, dq, and dr.

100

80

C

C

0

40

Q.

20

Q)

0

" dq = 81.07

................. :..... d! 6!39

- ,"" dp= 33.39

' ' ' ' I ' _ ' ' _ _ J _ l _ _ _ l J _ , _ I , , J _ l , _ , _ I _ _ _ _

0 2 4 6 8

Iteration Count

Figure 6. Computing the required CP differences using equal variances.

As shown in figure 7, the separation distances vary somewhat when equal variances are used

(b=32.17 ft).

14
i

o'1

U

t-"

o3

tO

r_

ta

(.9

100

8O

o

60

40

20£

0

0

r i

................. !....................I

u u i i i u i i u I I s u u i i i i i I I u l I I I I I I l i i U l I I u n I I i I u u I i I

200 400 600 800 1000

Altitude (or Vertical Characteristic Length) - ft

Figure 7. Separation distances vs. altitude for equal variances.

The CP separation distances for pitch and yaw are fairly constant for high-altitude flight, although

for spectral correspondence with the conventional formulation the required CP values for the FEAST

model are too high (not physically credible) for typical aircraft geometry. The required separation

distance for roll correspondence is also too high, except that it appears to approach a reasonable

value for high-altitude flight.

Although the variances are independent of vehicle velocity, the autospectra are not. Selecting

both a low and moderate velocity case for illustration, the following six figures show the spectral

differences between the conventional and FEAST models, where the CP separation distances

have been computed from the equal-variance relationships. In these examples where Lw = 250 ft

(/-a, = 791.48 ft), the selected value for o"w = 1 ft/sec, which produces or,, = 1.468 ft/sec. The

first three figures are for a velocity of 5 ft/sec.

_6

O)

e-_

O

<

O

t'r"

l I u I i Ill u I

10.2

' ' '''"i ' ' ' =''"I ' ' ' ' ''"

10"1 100 101

Frequency- rad/sec

Figure 8. Roll autospectra (V = 5 ft/sec) using equal variances.

15
Q.

89

o

<

(..,

Q.

10.2 - i

10"3- '_ .............

10., ........................ J .................

10.5 -- . J ............... :-.. _. ,

106 -- _ -- 32"17 ft' d" = 81"07 ft I .....

[L w= 250 ft, aq = .0110 rad/secJ

10 .7

10.3 10.2 101 10° 101

Frequency- rad/sec

Figure 9. Pitch autospectra (V = 5 ft/sec) using equal variances.

10"2

10 .3

10.4 -

0

10.5

>.-

10 .6

I I I I I I II t I I I I I III I I I 1 I I | II

10.2 101

Frequency- rad/sec

_r_ e)

l I l I I I II

10 0 101

Figure 10. Yaw autospectra (V = 5 ft/sec) using equal variances.

The second set of three figures is for a velocity of 100 ft/sec.

09

e'_

¢a')

O

<

O

IV"

10 .3 --

10.5 ............................... ([Dpc(e))

_0.-I_--_ _,,,_o--_._o. ..................... _-... ....

Lw= 250 ft, Op = .0150 rad/sec _

108- ' [ ....... / ........ I ........ 1 ......

10 .2 10 "1 10 0 10 1 10 2

Frequency- rad/sec

Figure 11. Roll autospectra (V = 100 fdsec) using equal variances.

16
1 0 .4

O

10.6 --.

_L 10 .7 - "

10 .6 -

10 .2 10 1 10 o 101 102

Frequency- rad/sec

Figure 12. Pitch autospectra (V = 100 ft/sec) using equal variances.

lO-6...................... _ ...... i

_ 10.6

. . b = 32.17 ft, dr= 61.39 ft ......................

>- 10 .7

106 ....... I '/ I I I I I I II I I I I I I II

10 .2 10 1 100

Frequency- rad/sec

(I)rs((l)) ......

10 _ 10 2

Figure 13. Yaw autospectra (V = 100 ft/sec) using equal variances.

Using Reasonable CP Separation Distances

For an aircraft with a span of 32.17 ft, the computed dp, dq, and dr required to obtain equal

variances are clearly too large. An example of an aircraft that has this particular wing span is the

XV-15. Data is available for the center-of-pressure distances for this aircraft. Considering the CP of

the fuselage given by 24.416 ft, for the XV-15 the differences in these distances are

dp = 17.08ft

dq - 46.666-24.416 = 22.25ft

dr = 47.501-24.416 = 23.085ft

b = 32.17 ft

We may investigate the changes in FEAST's autospectra using these values. As above, both low and

moderate speed cases are shown in the following six figures. The first three figures are for a velocity

of 5 ft/sec.

17
Q.

to

o

<

IT10"2

10 .3

10 .4

10.5

10 .6

10 .7

m

Ib=32"7"'do--'7"°8"'_ _o" Ik . "%....

.......... I°oc=°'_° raOse°'_ ='02'3raO"ecl ............................ "%

10 .3 10 -2 101 100 101

Frequency - rad/sec

Figure 14. Roll autospectra (V = 5 ft/sec) using CP distances.

m

Q_

co

o

<

e13._

10.2 -

10.3 ....................................... i .................

_o-,............... s!'!,

1o-_........._ ._qc_.._ . "4/_ _

_o°t/l_=%::o,t,Z:.:::,.%E::°,secl .......... _._

1°7-4, ........ i ........ i ........ i ........

10 .3 10 .2 10 "1 10 0 101

Frequency- rad/sec

Figure 15. Pitch autospectra (V = 5 ft/sec) using CP distances.

£).

to

O

"5

,,(

t_

>.-

10 .2 ,o 1..................................... i.............................

10 .4

10.5

10 .6

10.7

m

m

m

10 .3

j ................... _r_,................... -..._ _.........

Jb = 3217 ft' dr= 23"09 ft' Lw= 250 ft I- : ./l "_A"

| i ' ' ' ill] ! | | I I Ill I ] l I I I III I I I ! I I I II

10.2 10"I 100 101

Frequency - radlsec

Figure 16. Yaw autospectra (V = 5 ft/sec) using CP distances.

18
The second set of three figures is for a velocity of 100 ft/sec.

¢j

09

o'J

<

O

EE

10 .3 --

10 .4 --

10-_ -

10 .6 -

10 .7 --

10 .8 -

10 .2

................. i........... _c(_O)

...... iOoo=.o_ora_,sec, oo=0_1_ra_secl ........................................ _

u u i I IIII I i I I 1 I III i I I I I I III I I I I I I I II

10 1 10 ° 101 10 2

Frequency- rad/sec

Figure 17. Roll autospectra (V = 100 ft/sec) using CP distances.

10 s _

0

10.6

10 .7

108 ' ' ' ' .... i ' ..... "1 ' ' ' '''"1 ' ' ' .....

10 .2 10 1 10 ° 101 10 2

Frequency- rad/sec

Figure 18. Pitch autospectra (V = 100 ft/sec) using CP distances.

¢-_

o

<

>-

1 0 .4

10 s _

10 .6 --

1 0 .7 --

10 .8 -

10 .2

_tb_-32,7., _,_-23o,,t _w:2_o. I........ i............... tt_

I°'c: %--ra "ct VV

I I I I lliS I I I l i i'll i I I i a lili I s I I I I I II

10 "_ 10 0 101 10 2

Frequency - rad/sec

Figure 19. Yaw autospectra (V = 100 ft/sec) using CP distances.

19
Fromthesefiguresweseethatthepitchandyawspectrain theFEASTformulationareextendedin

thefrequencydomain,andtheyhavelargerdispersionsbecauseof this. Theroll spectrumhasless

bandwidth,buta largerdispersion.Thesephenomenaarewell within apilot's bandwidth,andwould

thusimpactsimulationresults.Theresultantrotationaldispersionsarecomparedto theconventional

dispersionsin thefollowing threefigures,o"c is the conventional model's dispersion and crs is the

FEAST model's dispersion.

0.06

0.04"

i

¢-

0

m o.o2

{3.

(/)

o 0.00

0

'll''l'l''il'''l'W,,I,,,,I,,,_l,,,,I,,,,I,,,li,,,,

200 400 600 800 1000

Altitude(orVerticalCharacteristic Length)-ft

Figure 20. Roll dispersion.

0.06

, 0.04

c0

_ 0.02

(D

O_

IIIIIIIIIiiiiiiiiii

0

JllalllllJll_Jllll]JltlJlJtll

200 400 600 800 1000

Altitude (or Vertical Characteristic Length) -

Figure 21. Pitch dispersion.

0.06

, 0.O4

E

0

_ o.o2 G.)

Q.

u'}

D

0.00

0

,llll,_,,i,,,llllWl ,',lll,l,l,l,,lll,,I,,,,11111

200 400 600 800 1000

Altitude (or Vertical Characteristic Length) - ft

Figure 22. Yaw dispersion.

20
Step Responses through the Turbulence Models

The conventional rotational turbulence model uses only the wing span parameter, where the value

b = 32.17 ft corresponds to the XV-15 aircraft. When the FEAST model is used, and it is

statistically required to produce the same angular dispersions as the conventional model, unrealistic CP differences were computed, given by dp = 33.39 ft, dq = 81.07 ft, and dr = 61.39 ft.

For the XV-15 aircraft, the smaller values of dp = 17.08 ft, dq = 22.25 ft, and dr = 23.09 ft are

known, where the span length is also given by b = 32.17 ft. These three configurations produce

three different step responses, given in the following three figures. These responses were created

at an altitude of 250 feet, and a vehicle velocity of 5 ft/sec.

U

i

o

0.020 - : :

:,...................... _,. i p_ (d=17.08)

0.01 5 .................. "...................... ""'. ......................... _ P .....

,_z,., ......... z i

0.010 ......... '" _- '" Ps(dp =33.39) .....

.......

0.000 -- ,,,,i .... I .... I .... I .... i .... .... I .... I .... I ' ;'i']"

0 20 40 60 80 100

Time - sec

Figure 23. Roll rate, step input.

t9

tO

"O

I

e_

t'-

0.08 -

0.06-

0.04

0.02

0.00 -

-0.02

0 20 40 60 80 100

Time - sec

Figure 24. Pitch rate, step input.

21
0.02 -

_.i

Q; 0.00

-0.02

I .

-0.04

_,,

-0.06

(_

>_

-0.08

'1' : °, ........ i i

i .............

rc / '""-_--------_r s (d=23.09)

Illlll JllllllllWllllllJllllllllllllllgl lllllllll

0 20 40 60 80 1O0

Time - sec

Figure 25. Yaw rate, step input.

From an examination of these figures the conventional model is shown to be a low frequency

approximation to the behavior computed in the FEAST model. The artifact of setting the dispersions

equal and then computing the required CP separation distances produces somewhat similar pitch and

yaw responses, but FEAST's roll response is larger.

In all cases it is seen that the conventional model has less rotational activity than the FEAST model,

and this is especially true if the correct CP distances are used.

Although step responses are not particularly revealing in terms of responses to random data, these

figures do indicate that the conventional model's rotational responses may only be valid for low

frequency behavior.

Time Response Comparisons

In the following material the XV-15 aircraft's CP separation distances are used in the FEAST model,

and the XV-15 aircraft's span length is used for the conventional model. The velocity is 5 ft/sec, and

time histories are presented for two different altitudes. Identical random sequences were used to

drive both formulations.

22
U

Ill

"O

i

o -0.1

-0.2

°21 II Pc

0.1 -- Ps

0.0

(_) Ih = 20 ft I

I''''1''''1''''1''''

0 5 10 15 20

Time - sec

0.06

O.04

_. 0.02

' 0.00

'_" -0.02

o

-0.04

-0.06

b h f

0 5 10 15 20

Time - sec

U

(/)

i

cr,

0.2

0.1

0.0

-0.1

-0.2

0

........... i J qc

Ih--2o,t]

-r'rrq-rr'rr -rrrq-rrr-r

S 10 1.5 20

Time - sec

0.06 qc

0.04

o.oz

i

o.oo .........._

-0.02

_- -0.04 -0.06

0 S 10 15 20

Time - sec

0.2

U

0.1

-0

s

0.0

-0.1

-0.2

I

(e)

i i i I

0

i -- rc

h 20 ft

I I I I I I I I I ] I I I I

5 10 15 20

Time - sec

0.06

_- °.°' 0.02 1..................... ...............

o.oo - _r,

-0.02 .................................. ""

>-

-0.04

-0.06

I''''1''''1''''1''''

0 5 10 15 2O

Time - sec

Figure 26. Rotational turbulence.

23
It shouldbenoticedthatthescalesusedin thesefiguresfor thelow altitudecasesaremuchlarger

thanthoseusedfor thehighaltitudecases.

TheFEASTmodelroll behaviortendsto approachtheconventionalmodelroll behaviorfor high

altitudeflight. Pitchandyawresponses,however,aregenerallylargerfor theFEASTmodel.

In reference1,a techniquewasusedto distributerandomturbulenceovertheUH60helicopter's

rotor disc.It ensuredthatthedispersionwasidenticalatall rotorbladeelementsby usingan

algorithmthatis equivalentto theonedevelopedhereinfor roll activity.It did not assumeany

particularcorrelationform.BecausetheUH60 simulationis adistributedaerodynamicmodel,

individual elementforceswerecomputedfrom extensiveaerodynamictables,andmomentswere

inducedfrom therotatingframegeometricalrelationships.Theresultsof thisSORBETsimulation

model,which hadextensivepilotedtesting,maybeusedto revealafeaturein commonwith the

FEASTmodel.Althoughpilot opinionwasgenerallyfavorablefor theSORBETmodel,the

following quoteindicatesthattheturbulencedispersionshouldbereducedfor low velocity flight

(whichimplies low altitude).Fromreference1:

"Pilots generallyagreedthattheturbulencemagnitudeshouldbea functionof velocity.

Above40knots,thesubjectivevaluesfor light andmoderateturbulencewereselectedas

crw = 5 ft/sec and 8 ft/sec, respectively. Near hover, pilots selected standard deviations that

were half these values. The reason for this is unknown, and should be investigated further."

The FEAST model's increased angular dispersions of figures 20-22 are thus an approximation to

what would occur if a distributed aerodynamic model were used, and the objectionable magnitudes

may be isolated to the translational spectra at low velocities and altitudes. The integrals of these

spectra (the translational variances), however, are selected, and they are independent of velocity and

altitude. Hence, it is suggested that the translational variances themselves should have some

functional decrease with altitude (in the low altitude region). This is strongly suggested (ref. 2) in a

graph of eYw versus altitude, which is unfortunately applicable to "medium/high altitude only."

Clearly, further research is required to determine the proper functionality for low altitude flight.

Discrete Gust Contributions

The literature is replete with papers dealing with mapping the environment around ships (refs. 11-13).

Also, current projects at Ames Research Center include terrain mapping of significant features of a

visual data base to correlate gusts to these features using computational fluid dynamics (ref. 14). The

mapping functions are exceedingly complex; these complexities are avoided here by an example using

a simple terrain related gust model. The FEAST model produces complex atmospheric/vehicle motion

in response to a simple terrain related model.

Gusts are assumed to be created in a flame of reference corresponding to the visual model, and

transformed to the instantaneous distributed centers of pressure. Terrain correlated gusts may be

formulated to produce representative angular activity if they are applied at the pertinent CPs. A

system that uses this geometry requires the transmittal of seven gust components from the

atmospheric model to the aircraft model. These consist of one u-gust, two v-gusts, and four w-gusts.

The fuselage CP and the CPs of the tail and wing are the pertinent points. The quantities required

24
eachcycletime are: ug(O), vg(O), vg(dr), wg(O), Wg(+ l dp), Wg(-l dp) and wg(dq). The centerof-gravity translational discrete gusts are then the values with the arguments of zero (the CG is very

close to the CP). The rotational contributions are then given by the equations

Pg = dp

qg

w ,(dq)

dq

rg

v (O)

dr

An example using these discrete gust equations is given in appendix 5.

The fuselage CP vertical discrete gust could be approximated by the average of the wing CP gusts,

which would reduce the required number of gust transmitted quantities to six.

Wg

+

2

The total contributions to the atmospheric model are given by sums of the stochastic turbulence and

discrete gust velocities, such as Ua(Z)= Us(Z) + us(z) and Pa(Z) = Ps(Z) + pg(z).

CONCLUSIONS

A generic turbulence and gust model has been developed for aircraft simulation, where both the

translational and rotational components are computed. This model, called FEAST, accommodates

stochastic turbulence and distributed discrete gusts as a function of the terrain, while the required

aircraft-specific parameters are related to the vehicle geometry. FEAST delivers rotational stochastic

components that are extended in the frequency domain beyond those delivered by the conventional

model.

New discrete techniques for the computation of rotational rates are introduced. A correlation

coefficient based upon the physical separation of aerodynamic centers, and a two-dimensional

vertical turbulence field are used in the computation of roll rate. Filtered random noise sources are

propagated in time by use of interpolated transport delays. For the computation of pitch rate, the

Gaussian centroid vertical velocity of the field influences the centers of pressure of the fuselage and

the horizontaltail as a function of vehicle velocity. Similarly, the vertical tail's center of pressure is

used for the computation of yaw rate.

25
The FEAST model should be useful for low speed and low altitude flight simulation. It uses the

specific vehicle geometry in the computation of rotational responses to turbulence, and gusts as a

function of the terrain. The state-space techniques given for the discrete implementation of this

model are computationally efficient, and deliver stable outputs for all stable inputs regardless of the

required cycle time.

26
APPENDIX 1

Interpolation and the Circular Buffers

The pitch center-of-pressure difference dq may be larger than the yaw center of pressure difference

dr, as is assumed in this appendix. Regardless, the larger value is used to compute the minimum

aerodynamic velocity that must be imposed in order to retain all of the required vertical and lateral

velocity data in buffers of size N. In order to accomplish this we have

V > dq

(N- 1)T

The velocity used throughout this paper is assumed to conform to this inequality, resulting in a

minimum velocity of a couple of feet per second using reasonable buffer length and cycle time

parameters. We may then develop an integer index kq with the following limits5:

with remainder flq given by

dq kq < 1

O < flq - VT

The lower velocity limit assures us that if the first ( n = 0) cell of a buffer contained the current time

point's fuselage CP data, then we could interpolate within the buffer to find the horizontal tail's CP

vertical velocity data by use of

W-vy ( dq ) = [1- flq]w(kq)+ flqw(kq + 1)

and find the required end-point values, regardless of the aerodynamic velocity (above the minimum).

Similarly, if the yaw CP difference is smaller than the pitch CP difference, then the same minimum

velocity is used such that the lateral inequalities become

dr k r < l

O <- flr - VT

and the vertical tail's CP lateral velocity data may be obtained from

5 Using the "floor" or "intier" operation (least integer).

27
 -Vfdr ] = [1-- _r]v(kr)+ _rv(kr + l)

J

The buffers are initially full of velocity data. This may be accomplished by a negative-indexing

initializing scheme. Using the filtered random inputs, first the last cell ( N - 1) is filled followed by

( N - 2), etc., until the index becomes (0). In simulation this requires N passes through the

turbulence equations (initial condition mode) before realtime operations are begun (operate mode).

The first storage cell in realtime is then also 6 N - 1, and thereafter, a decrement continues each cycle

time for the entire run, with a modulus of N.

The "right arrow" is here read "becomes"; at the transition point this produces

0 _ (O-1)Mod(N)= (O-1)-N[_-_] =-1-N[--_-J = N-1

effectively completing the circle. The circular buffers each require only one store operation every

cycle time, into the cell designated as n. The next higher indexed cell (ModN) is always the cell

from the previous deposit, so this scheme is sufficient for linear interpolation of the delayed vertical

and lateral velocity values.

From the above preliminaries, the entire procedure is given below in twelve steps for depositing the

current filter outputs into the circular buffer (where these values are the fuselage CP velocities), and

retrieving the vertical and lateral velocity at the horizontal and vertical tail centers of pressure.

Step(l) n _ n-1 N-1 n>l<

Step (2) Deposit current filter outputs into w(n) and v(n).

Step(3) kq = [---_J

dq

Step (4) _q - VT kq

I n+kq I

Step (5) mq =- n + kq - g[-------_-J

6 Not really required as long as the buffer is full and the correct index is retained in transferring to operate mode.

Initializing just requires at least N cycles.

28
Step(6)

Step (7)

Step(8)

Step (9)

Step(lO)

Step (l l)

Step(12)

mq

kr

J_F

m r

0 q + 1

dr kr

VT

mq<N-1

mq=N-1

m r or+l mr m <N-1 r = N-1

= [1-flq]V(mq)+flqV(m;)

= [1-_r]V(mr)+_rV(m*r)

29
APPENDIX 2

z-Transform Solutions

For discrete implementation of the translational filters using either of the two models a zero-order

hold formulation is typically used, as shown by the three z-transform relationships:

F(z) = Z{

I _T] '-eSTI 1 _7,1-eL,,

s f(s) = _ _rvl z-_-E-vTJ

G(z) = z g(s)}:

s

+C2)

( 12 V T

Z--e L_,

H(z) fll_<ST,, t

II 2

= Z s )h(s) = v T

Z--e Lw

where the coefficients are:

V V

CI : l-e L:'T+(w/3--1)_-_e -U,T

C2 = -e L,,vT[1-e _,_T L_ + -1

V V

C3 = l-e LwT+(.f3--1)--_e-_ T

G = -e L_ 1-e Lw + -1

The difference equations that result from the z-transforms are driven by zero-mean, unity variance

discrete Gaussian noise sequences Ok (k= 1, 2, 3, 4) with amplitude range {-3.45, 3.45}, which

preserves 99.94% of the continuous total probability density function. These derivations assume the

30
utilization of the realtime XNORM © Gaussian generation software, with zero mean and unity

standard deviation.

The conventional model's rotational terms are developed below to make a comparison with angular

outputs of the FEAST model. The rotational terms require three z-transforms. Whereas the roll L(z)

term uses the zero-order hold formulation, as was used in the development of the translational ztransforms, th_ pitch M(z) and yaw N(z) formulations use the triangular hold, because they are

computed, in sequence, after functions that use the zero-order hold. This preserves the proper

input/output phase relationships between all variables.

M(z) {ll:x lae ---f-m(s)

l _ e-nVT_4b l(z_ l)

L,z, = t

(Zc/4b)l//6ff-_II_e-rCVT/"4b)

 lz- I

N(z) Z ('I-_-_T1"e-ST _n(s)I ] =

For discrete implementation a normalizing factor is required to preserve the total power integrals

(cy 2). Where the rli(z ) are unity-variance, zero-mean Gaussian variables, the translational ztransforms are given by

Uc(Z)=Cru F(z)rh(z)

Vc(Z)= Crv_-_yO(z)o2(z)

Wc(Z ) = _w_-_ n(z)O5(z) '- CYw_z z n(z)[O4(z)+ 03(z)]

and, with a sampling switch between translational and rotational filters, the rotational rates become:

31
qc(z) = aw H(z)M(z)Os(z)=aw H(z)M(z)[O4(z)+03(z)]= M(z)wc(z)

rc(Z) = a,,# C(z)g(z),72(z)= g(z)Vc(Z)

Note that these rotational filters are not used in the FEAST model.

Difference equations are easily created from these z-transforms. For example, the longitudinal

turbulence velocity becomes

Uc(n + l) = e

L. uc(n)+Cru___V____l_ e r/l (n)

State space solutions to the turbulence filter equations are preferred because they produce stable

outputs from all stable inputs, independent of the size of the cycle time parameter, and the inputs to

these filters have considerable frequency content. Additionally, the zero-order hold solution

produces an exact solution to a step input, and the triangular hold produces an exact solution to a

linear input. Note that the discrete coefficients must be recomputed whenever velocity or

characteristic length change. Due to the fact that these are "slowly varying" quantities, however, a

practical implementation technique is to distribute the workload of the coefficient-update operation

over a number of computer cycles. In this way the workload becomes competitive with that of less

robust algorithms, and the entire coefficient set is updated in a few computer cycles. This technique

has been used at Ames Research Center for 20 years.

32
APPENDIX 3

Conventional Model Rotational Relationships

The system functions, autospectra, and variance are given for the conventional model.

Conventional Model Roll Relationships

The conventional (MIL SPEC) formulation's roll response may be written in terms of the

independent unity-variance random variable r/6(z )

pc(z) = Ow yyP(z) 6(z)

For this conventional formulation the roll system function is

pc(S)

Crw O_.8V(zc/4b) 7/66

/-_3(s + :rrV/4b)

The roll autospectrum given by

_ pc((O)

0.8cr2 V( zr,/4b )7/3

/4_[092 +(_V/4b) 2]

and its variance is computed by the frequency integral from zero to infinity

O-2c

= IO _ pc 0.8Cr2wVQr/4b) 7/3I009 do) (r-o)d (-o = _- 2 +(zrV/4b)2

oso wVC _azv ] O.4cr2zr (rr)_

- I_ _.4bJ [2(n:V/4b)J = L_ _.4"-bJ

Conventional Model Pitch and Yaw Relationships

Where a sampling switch is assumed between the linear and rotational filters, pitch and yaw for the

conventional model are given in z-transform notation by

33
qc(zl-- w(z)M(z)os(z)

re(z) = v(z)U(z)O2(z)

= CTw HzMz It It,5llZ

= crv_[-_G(z)N(zlrl2(z )

The system functions for these rates may be written

qc(s) = CrwW(S)m(s )

r_(s) = cr,,v(s),,(s)= -

4b _zcl__, '

(s + 7cV/4b)(s + V/Lw) 2

Ov,_f_t, , +[3v v.4_Z.)s

(, +7rW3b)(s+ V/L,,)2

The autospectra are given by:

f_)rc ((D)

3°gv'_o,2(,o2+v2/3_,,,)

16bZLw

((.02 + E2V2/16b2)((_02 + V2/L2v) 2

°'2V_ 0)2(0)2 + V2/BL_,)

3bZ/w

(09 2 + yt2V2/9b2)((.o 2 + V2/L_) 2

and from these expressions the variances are:

1_2c = I_o¢_qc(OJ)do)

(y2c = I_Otf_rc (O) )d o)

Ow,_ 2 2/_3/r/--w (q;-+ 2)

2

v \b J

34
APPENDIX 4

FEAST Model Rotational Relationships

The system functions, autospectra, and variance are given for the FEAST model.

FEAST Model Roll Relationships

Using the properties of random variables, an expression for roll rate in the FEAST model may be

developed. Where dp is the difference between the centers of pressure for the left and right wings,

vertical velocities applied at these points produce a roll rate given by

Wg -- W r

Ps --

dp

These velocities at the centers of pressure are correlated. However, using the previously developed

correlation coefficient the roll rate may be expressed as a function of the two uncorrelated vertical

velocities wL and w k

Ps

w -wr 14 p[wL-w ]

n

dp dp

Hence, FEAST's roll rate in z-transform notation becomes

(Yw_Jl- p P

dp dp

and the system function for roll rate is given by

ps(s) = dp G

-- (s+ WL ,)2

FEAST's roll autospectrum is thus a function of the vertical autospectrum

2(l_p)_w(o_) :

- d7

6cr2 V(1- p)(o_ 2 + V2/31__.)

/rJ_,wd2((.o 2 + V2/L_,) 2

although roll response is uncorrelated with vertical response. The roll variance is similarly computed

by a frequency integral from zero to infinity.

(712,s = Io_P PS(Og)dw -2(l_/2P);o_W(Og)d°9 = 2crzw(1-p)d2

35
FEAST Model Pitch and Yaw Relationships

It can be assumed that w(z) is created at the fuselage center of pressure without loss of generality if

transport delays are then imposed, which are equivalent to the distances to the horizontal stabilizer

(dq, for pitch rate) and vertical stabilizer (dr, for yaw rate). The time shifts are a function of the

vehicle velocity

t q = d q/vV

tr = dr//V

which, in discrete simulation, are equivalent to selecting and interpolating filtered velocity values

from circular buffers as shown in appendix 1. The pitch and yaw rates in z-transform notation are

then given by,

rs(Z, = dr °vdr _T_Z /7"_llG(z)Ti2(z)

These produce the system functions

(1- e-_'_) crwh(s)

qs(S) - dq

e -st_ - 1)

r_(_) - dr _g(s)

O'w(1-e -stq )_fV///;rcl_, (s + V/'q_ l-.w)

aq(S+ v/z,,) 2

dr(sWW/Lv) 2

which may be written generically

cAs)

O'(1-e-Sr)ff3a/,/_(s +a/if-33)

d(s+a) 2

with a generic autospectrum given by

,:,_(co)

3acr2(2_e-J °Jr -eJ_°r)(co 2 + a2/3)

_(o_: +a:)_

6ao-2(1- coscor)(co2 + a2/3)

,_(o_ +a_)_

The individual autospectra thus become

36
6a2V(1-coscodq/V)(co 2 + V2/3_,,)

/l:'Lwd2((.o 2 + V2/L_) 2

6f2V(1-cosOJdr/V)(09 2 + V2/3L_)

The generic variance is given by the integral

- 6af2 - cosgor) .o22 1 _5__ 2.|do;

f2s = s((-O) d('O ;_/2 +a 2 ((.02+ a2) j

/rd2 (_02+a2 2,, 3 ((.02+a2)2 0)2co_ + a2 +--52,.so.) (0.)2+ a2)2 do.)

- 6,...{= ruJ2 Ta-3 2..r,.,q L4a 3_[ ,-.-.-0" 2a +g2..,[- L4a" i,+..)e-°_-]}

f2[(a'r-2)e -at +2]

From this expression the FEAST model's pitch and yaw variances are given by

<:,-Ts= -a7-,'Lt/_- +2l

0..2= _Lt/,_.,._2)e-";_,, f_ Pro./ +2]

37
APPENDIX 5

Discrete Gust Example

To illustrate the influence of the distributed points of application in relationship to a terrain model, a

simplified canyon (or ship deck) is constructed. This model uses two separate terrain altitudes, in a

Cartesian plane, as shown in figure 27.

h(terrain) = 50 ft

PLATEAU

V

[xAO),yA0)]

k[ : i(,[erriii)-0:fti i i

: : :6Ah,_oh: : : :

v

Figure 27. Terrain elevation.

An aircraft travels at a constant speed (V), at an azimuth angle ( gt ), and at a constant altitude. The

altitude of the aircraft above the plateau is 50 ft, and over the canyon it is 100 ft. For illustrative

purposes, neither the aircraft trajectory nor attitude is influenced by the gusts. The disturbance model

is constructed on this plane such that a vertical (upward) gust of 5 ft/sec is encountered when the

aircraft arrives at the canyon wall. This occurs at xi(t) = O, where i = (F, HT, VT, RW, LW) denoting

the fuselage CP, the horizontal tail CP, the vertical tail CP, the right wing CP, and the left wing CP.

As a function of time (t) the points of application are given by

XF(t ) = XF(O ) + Vtsin _

XHT(t) = XF(t) - dq sin I/t

xvT-(t) = xF(t) -- d r sin Ikt

XRw(t ) = XF(t ) + 1⁄2dp coslg

XLW(t ) = XF(t ) - 1⁄2dp cosl//

yi(t) values are not required in this example because the y-gust is assumed to be a function of the wgust, which is only a function of xi(t ). For convenience, the vertical gust decays linearly to zero

38
when.x,i(t) = 100 ft. Also, a positive y-gust of one-half the magnitude of the vertical gust is assumed

to exist that is independent of the yi(t) value.

Wg(i) = 100 J

xi <0

0 < xi < 100fi

X i > 100ft

Ug(i ) = -- _ Wg(i ) COS I]./

1

Vg(i ) = "_ Wg(i ) sin

The discrete gust profile given above leads to asymmetric responses when the canyon is approached

at an angle V = 45°" In order to investigate these responses the initial condition xF(O) "- -17.68 ft is

used, where V = 5 ft/sec. The rate of closure on the canyon becomes 3.536 ft/sec, so that it takes

exactly 5 seconds for the fuselage CP to arrive at the canyon. Due to the geometry, (1) the right wing

arrives first at t =3.29 sec, (2) the fuselage CP arrives at t = 5.00 sec, (3) the left wing arrives at

t=6.71 sec, (4) the horizontal tail arrives at t-9.45 sec, and (5) the vertical tail arrives at t=9.62 sec.

The vertical gust velocity that occurs at each of these aircraft locations is given as a function of time

in figure 28.

(.,;

¢J3

i

,i.a

(/)

(..9

t,)

>

-2 n

-4

-5

0

| |

: 1 J )

: I , i

: : , ]

: : i )

: : m i

; f ]

; ._r _ ,

ililill

10

_-_.,_"- Vertical Tail

_iv --_'- Horizontal Tail

, ,, ,i, ,,,I,,, ,i,,,,i,, , , i,,, ,I

20 30 40

Time - sec

Figure 28. Vertical gusts at aircraft locations.

The gusts influence the fuselage CP in translational space. The onset of these gusts thus occur at

t = 5 sec, and they are as shown in figure 29, where the u-gust and v-gust are half the magnitude

of the w-gust.

39
{-4_-z ,,',i.... w ......... vg.... i......... ............',,,

nn

0 10 20 30 40

Time - sec

Figure 29. Translational gusts at the center of gravity.

Induced rotations occur as the various aerodynamic centers become immersed in the gusts. Since

the right wing arrives at the canyon at t = 3.29 sec, the roll response begins at this time. Note that

partial recovery occurs at t = 6.71 sec, when the left wing also becomes immersed in the gust field.

The differential gradient occurs due to the different vertical gust magnitudes, out to a point where

neither wing remains in the gust field. The pitch and yaw responses occur almost simultaneously

because dq = dr. These responses are shown in figure 30.

"o

i

tO

.m

4,.-i

t_

ft, -0.3 -

.3 u ,-._

o.z........E "'_p_ ................!..........

....... i ....... _.......... .... rg .............. ................. : ................

o.o- 1 -o..--L-_::.:_::.:_`.__::_:_::_:_:_::::_i_._:_-_"_.-_._!_._;_._._

] ..... , ................ . ..................... .

-0.1 ................ ! ........... i .......................... ..................................

J

I I ! I I I l l t I i i ! I I l J J l I i I i ! I ; l J I [ l l l i I l l l I

0 10 20 30 40

Time - sec

Figure 30. Rotational response to discrete gusts.

From this example we see that the FEAST model produces atmospheric variations based upon the

specific vehicle geometry in response to terrain-related discrete gusts. This is important in evaluating

responses of different aircraft to the same terrain related (or ship) turbulence and gust profile.

Refinements to the FEAST model await the conversion of an area of CGI-terrain data to an

appropriate form for analysis using computational fluid dynamics, and the application of the

resultant flow-field data to various aircraft simulations.

40
REFERENCES

1. McFarland, R. E.; and Duisenberg, Ken: Simulation of Rotor Blade Element Turbulence. NASA

TM-108862, Jan. 1995.

2. "United States Department of Defense: Flying Qualities of Piloted Aircraft." Military

Specification MIL-F-8785C, Nov. 1980.

. Tate, S. J.: A Dynamic Challenge: Helicopter/Ship Interface Simulation-Development,

Integration and Application. AGARD Conference Proceedings 577, Braunschweig,

Germany, May 1995.

4. Woodfield, A. A.; and Tomlinson, B. N.: Ship Airwakes-A New Generic Model for Piloted

Simulation. AGARD Conference Proceedings 577, Braunschweig, Germany, May 1995.

. Chalk, C. R.; Neal, T. P.; Harris, T. M.; Pritchard, F. E.; and Woodcock, R. J.: Background

Information and User Guide of MIL-F-8785B(ASG), "Military Specification-Flying

Qualities of Piloted Airplanes," Wright-Patterson Air Force Base Technical Report AFFDLTR-69-72, Aug. 1969.

. Etkin, B.: Theory of the Flight of Airplanes in Isotropic Turbulence-Review and Extension.

NATO Advisory Group for Aeronautical Research and Development, Report 372, April

1961.

. Prasad, J.V.R.; Riaz, Jamshed; Gaonkar, Gopal H.; and Yingyi, Dang: Real Time

Implementation Aspects of a Rotorcraft Turbulence Simulation Method, presented at the

American Helicopter Society 49th Annual Forum, St. Louis, Missouri, May, 1993.

. Schattenmann, W.: Simulation of Atmospheric Turbulence, in Particular its Influence on Aircraft

Yaw Motion. NASA Technical Translation F-15,802, Deutsche Gesellschaft fuer Luft- und

Raumfahrt, Cologne, Germany, Report DLR-Mitt-70-12, Nov. 1970.

9. Bendat, Julius S.; and Piersol, Allan G.: Engineering Applications of Correlation and Spectral

Analysis. John Wiley & Sons, New York, 1980.

10. Conte, S. D.; and de Boor, Carl: Elementary Numerical Analysis, An Algorithmic Approach.

McGraw-Hill Book Co., New York, 1965.

11. Nave, Ronald L.: Development and Analysis of a CVA and a 1052 Class Fast Frigate Air Wake

Model. Naval Air Development Center, Warminster, Penn. Report WF-41-400-000, P. E.

No. 62241N, Sept. 1978.

12. Fortenbaugh, R. L.: Mathematical Models for the Aircraft Operational Environment of DD-963

Class Ships. Vought Corp. Paper 2-55800/8R-3500, Sept. 1978.

41
13. Garnett, Theodore S. Jr.: Investigation to Study the Aerodynamic Ship Wake Turbulence

Generated by a DD963 Destroyer. Boeing Vertol Co. for the Naval Air Development Center,

Warminster, Penn., NADC Report No. N62269-78°C-0097, March 1980.

14. "Atmospheric Turbulence Modeling for Real-Time Simulation of Nap-of-the-Earth (NOE)

Flight," Systems Technology, Inc. contract NAS2-14105 to NASA, Ames Research Center.

42

Form Approved

REPORT DOCUMENTATION PAGE oueNo.0704.0788

Public reporling burden for this collection of information is estimated to average 1 hour per response, including the time for reviewing instructions, searching existing data sources,

gathering and maintaining the data needed, and completing and reviewing the collection of information Send comments regarding this burden estimate or any other aspect of this

collection of inlormation, including suggestions for reducing this burden, to Washington Headquarters Services, Directorate for information Operations and Reporls, 1215 Jefferson

Davis Highway, Suite 1204, Arlinglon, VA 22202-4302, and to the Office of Management and BudgeL Paperwork Reduction Project (0704-0188). Washington, DC 20503.

1. AGENCY USE ONLY (Leave blank) 2. REPORT DATE 3. REPORT TYPE AND DATES COVERED

February 1997 Technical Memorandum

4. TITLE AND SUBTITLE

Finite Element Aircraft Simulation of Turbulence

6. AUTHOR(S)

R. E. McFarland

7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES)

Ames Research Center

Moffett Field, CA 94035-1000

SPONSORING/MONITORING AGENCYNAME(S)ANDADDRESS(ES)

National Aeronautics and Space Administration

Washington, DC 20546-0001

5. FUNDING NUMBERS

505-64-84

8. PERFORMING ORGANIZATION

REPORT NUMBER

A-975955

10. SPONSORING/MONITORING

AGENCY REPORT NUMBER

NASA TM- 110437

11. SUPPLEMENTARY NOTES

Point of Contact: R.E. McFarland, Ames Research Center, MS 243-1, Moffett Field, CA 94035-1000

(415) 604-3863

12a. DISTRIBUTION/AVAILABILITY STATEMENT 12b. DISTRIBUTION CODE

Unclassified-Unlimited

Subject Category 05

Available from the NASA Center for AeroSpace Information,

800 Elkridge Landing Road, Linthicum Heights, MD 21090; (301) 621-0390

13. ABSTRACT (Maximum 200 words)

A turbulence model has been developed for realtime aircraft simulation that accommodates stochastic turbulence and distributed discrete gusts as a function of the terrain. This model is applicable to conventional aircraft,

V/STOL aircraft, and disc rotor model helicopter simulations. Vehicle angular activity in response to turbulence is

computed from geometrical and temporal relationships rather than by using the conventional continuum approximations that assume uniform gust immersion and low frequency responses. By using techniques similar to those

recently developed for blade-element rotor models, the angular-rate filters of conventional turbulence models are

not required.

The model produces rotational rates as well as air mass translational velocities in response to both stochastic

and deterministic disturbances, where the discrete gusts and turbulence magnitudes may be correlated with

significant terrain features or ship models.

Assuming isotropy, a two-dimensional vertical turbulence field is created. A novel Gaussian interpolation

technique is used to distribute vertical turbulence on the wing span or lateral rotor disc, and this distribution is

used to compute roll responses. Air mass velocities are applied at significant centers of pressure in the computation of the aircraft's pitch and roll responses.

14. SUBJECTTERMS

Keywords: Simulation, turbulence, transport delay, autospectra

17. SECURITY CLASSIFICATION

OF REPORT

Unclassified

NSN7S40-O_-aBO-_SO0

18. SECURITY CLASSIFICATION

OF THIS PAGE

Unclassified

19. SECURITY CLASSIFICATION

OF ABSTRACT

15. NUMBER OF PAGES

48

16. PRICE CODE

AO0

20. LIMITATION OF ABSTRACT

Standard Form 298 (Rev. 2-89)

Prescribed by ANSI Sld Z39-18

298-102