George Kopasakis

Glenn Research Center, Cleveland, Ohio

Atmospheric Turbulence Modeling for Aero 

Vehicles: Fractional Order Fits

NASA/TM—2010-216961

December 2010
NASA STI Program . . . in Profi le

Since its founding, NASA has been dedicated to the 

advancement of aeronautics and space science. The 

NASA Scientifi c and Technical Information (STI) 

program plays a key part in helping NASA maintain 

this important role.

The NASA STI Program operates under the auspices 

of the Agency Chief Information Offi cer. It collects, 

organizes, provides for archiving, and disseminates 

NASA’s STI. The NASA STI program provides access 

to the NASA Aeronautics and Space Database and 

its public interface, the NASA Technical Reports 

Server, thus providing one of the largest collections 

of aeronautical and space science STI in the world. 

Results are published in both non-NASA channels 

and by NASA in the NASA STI Report Series, which 

includes the following report types:

• TECHNICAL PUBLICATION. Reports of 

completed research or a major signifi cant phase 

of research that present the results of NASA 

programs and include extensive data or theoretical 

analysis. Includes compilations of signifi cant 

scientifi c and technical data and information 

deemed to be of continuing reference value. 

NASA counterpart of peer-reviewed formal 

professional papers but has less stringent 

limitations on manuscript length and extent of 

graphic presentations.

• TECHNICAL MEMORANDUM. Scientifi c 

and technical fi ndings that are preliminary or 

of specialized interest, e.g., quick release 

reports, working papers, and bibliographies that 

contain minimal annotation. Does not contain 

extensive analysis.

• CONTRACTOR REPORT. Scientifi c and 

technical fi ndings by NASA-sponsored 

contractors and grantees.

• CONFERENCE PUBLICATION. Collected 

papers from scientifi c and technical 

conferences, symposia, seminars, or other 

meetings sponsored or cosponsored by NASA.

• SPECIAL PUBLICATION. Scientifi c, 

technical, or historical information from 

NASA programs, projects, and missions, often 

concerned with subjects having substantial 

public interest.

• TECHNICAL TRANSLATION. Englishlanguage translations of foreign scientifi c and 

technical material pertinent to NASA’s mission.

Specialized services also include creating custom 

thesauri, building customized databases, organizing 

and publishing research results.

For more information about the NASA STI 

program, see the following:

• Access the NASA STI program home page at 

http://www.sti.nasa.gov

• E-mail your question via the Internet to help@

sti.nasa.gov

• Fax your question to the NASA STI Help Desk 

at 443–757–5803

• Telephone the NASA STI Help Desk at

 443–757–5802

• Write to:

 NASA Center for AeroSpace Information (CASI)

 7115 Standard Drive

 Hanover, MD 21076–1320
George Kopasakis

Glenn Research Center, Cleveland, Ohio

Atmospheric Turbulence Modeling for Aero 

Vehicles: Fractional Order Fits

NASA/TM—2010-216961

December 2010

National Aeronautics and

Space Administration

Glenn Research Center

Cleveland, Ohio 44135
Acknowledgments

The author would like to acknowledge the support of the NASA Supersonics Project for the research conducted in this paper.

Available from

NASA Center for Aerospace Information

7115 Standard Drive

Hanover, MD 21076–1320

National Technical Information Service

5301 Shawnee Road

Alexandria, VA 22312

Available electronically at http://gltrs.grc.nasa.gov

This work was sponsored by the Fundamental Aeronautics Program 

at the NASA Glenn Research Center.

Level of Review: This material has been technically reviewed by technical management. 
NASA/TM—2010-216961 1

Atmospheric Turbulence Modeling for Aero Vehicles: 

Fractional Order Fits

George Kopasakis

National Aeronautics and Space Administration

Glenn Research Center

Cleveland, Ohio 44135

Abstract

Atmospheric turbulence models are necessary for the design 

of both inlet/engine and flight controls, as well as for studying

coupling between the propulsion and the vehicle structural 

dynamics for supersonic vehicles. Models based on the 

Kolmogorov spectrum have been previously utilized to model 

atmospheric turbulence. In this paper, a more accurate model 

is developed in its representative fractional order form, typical 

of atmospheric disturbances. This is accomplished by first

scaling the Kolmogorov spectral to convert them into finite 

energy von Karman forms and then by deriving an explicit 

fractional circuit-filter type analog for this model. This circuit 

model is utilized to develop a generalized formulation in 

frequency domain to approximate the fractional order with the 

products of first order transfer functions, which enables 

accurate time domain simulations. The objective of this work 

is as follows. Given the parameters describing the conditions 

of atmospheric disturbances, and utilizing the derived 

formulations, directly compute the transfer function poles and 

zeros describing these disturbances for acoustic velocity, 

temperature, pressure, and density. Time domain simulations 

of representative atmospheric turbulence can then be 

developed by utilizing these computed transfer functions 

together with the disturbance frequencies of interest. 

Introduction

This paper addresses the need for a model that simulates 

atmospheric disturbance, in both the time and frequency 

domain, over a wide range of altitudes and variations in 

atmospheric turbulence conditions, that is relatively easy to 

implement and representative of the actual fractional order 

nature of atmospheric turbulence. This is applicable to both 

propulsion system flow field type disturbances as well as 

vehicle gust loads.

Atmospheric turbulence have been studied for some time, 

especially in the field of Atmospheric Sciences, Nastrom 

(1985), and Fairall (1991), and various models have been 

developed. These models are primarily based on the so called 

Kolmogorov spectrum, originally developed by Tatarski (1961), 

based partly on the studies of turbulence by the Russian 

mathematician Andrei Kolmogorov (Kolmogorov (1941) (a) 

and (b)). The Kolmogorov spectrum has an energy level that 

approaches infinity as frequency approaches zero. This 

characteristic makes it difficult to implement these types of 

models in the time domain. A suitable approximation to the 

Kolmogorov model that is commonly used with a finite energy 

spectrum is the von Karman type model, originally referred to 

as the isotropic-turbulence spectrum, (Houbolt (1964)).

However, simulating the von Karman type models in the time 

domain is still problematic because of their fractional order. 

The Kolmogorov model has also been extended (Tank (1994)) 

to develop a baseline of atmospheric turbulence for the High 

Speed Civil Transport (HSCT). The Tank model also covers 

atmospheric acoustic wave disturbance modeling utilizing the 

von Karman spectral. Hoblit (1988) introduced the Dryden 

model approximation to the fractional order von Karman 

atmospheric model. But this model is second order compared 

to the 5/3 fractional order of the acoustic velocity atmospheric 

turbulence spectral. Thus, the Dryden model underestimates 

the atmospheric disturbances, increasingly with frequency.

To alleviate some of these difficulties, the Tank model for the 

von Karman approximation is utilized in this paper to derive an 

explicit electrical circuit analog of atmospheric disturbances. 

This circuit analog will act as a low pass filter. The circuit 

elements are explicitly computed as functions of atmospheric 

parameters, such as eddy dissipation rate and integral length 

scale. However, like the actual atmospheric disturbances, 

Nastrom (1985), the circuit order also turns out to be fractional

order, which makes it difficult to simulate. Thus, the circuit 

model is used as the basis in this development to derive integer

order transfer function (TF) approximations to the fractional 

order model. These TF approximations are a product of first 

order poles and zeros, which are determined as a function of the 

parameters describing an atmospheric disturbance. This

approach alleviates the manual process of hand fitting the

approximation every time an atmospheric parameter is changed.

Fairall (1991), Tank ((1994) and (1996)) and others approach 

atmospheric disturbance modeling probabilistically, as an 

exceedance for controls design purposes. That means that the 

probability of a time to failure is computed for controls design. 

This time to failure is associated with a controls design that can 

tolerate a maximum specified disturbance, and computes the 

probability, in terms of flight miles or hours, that an atmospheric 

disturbance will exceed this threshold. In this paper the approach 

will be to compute the worst case disturbance expected, and 

assume that the control system will be designed to handle this 

disturbance. Thus, this paper will not cover exceedance rates, 

which can be computed separately, based on the worst

disturbance expected and the specifics of the controls design.
NASA/TM—2010-216961 2

Kopasakis (2010) extends the work discussed in this paper 

to cover considerations for atmospheric turbulence 

specifications for supersonic propulsion systems. This work 

also provides an example for a supersonic inlet shock position 

controls design in the presence of atmospheric turbulence. 

This paper is organized as follows. First, the Kolmogorov 

form of the atmospheric disturbance spectrum is presented 

followed by the Tank model for the von Karman spectrum 

forms of the acoustic disturbances. This is followed with 

formulations of the equivalent fractional order TF 

approximations of atmospheric turbulence. Finally, time 

domain atmospheric disturbances are discussed for the derived 

TF approximations, followed by concluding remarks. The 

formulations covered in the body of the paper are based on an 

electrical circuit analog developed for different types of 

atmospheric disturbances, which is covered in Appendix B. 

While the detailed derivations of these formulations are 

covered in Appendix C.

Kolmogorov Form of the Atmospheric 

Disturbance Model

Tank ((1994) and (1996)) utilizes a Kolmogorov onedimensional locally isotropic atmospheric turbulence

spectrum, mathematically developed in Tatarski (1961), which 

represents the spectral density of a structured random field of 

atmospheric turbulence as

2 / 3 5/ 3 ( ) − St k = αtε k (1)

In this equation the quantity ε signifies the eddy dissipation 

rate, in units of (energy/(mass × time)) = (m2

/sec3), and k

signifies the wavenumber, in units of (rad/m) or (cycles/m). 

The subscript t in the atmospheric turbulence spectral density, 

St(k), signifies the type of disturbance. Based on these 

definitions, the units of St(k) are (m/sec)

2

/(rad/m) (Tatarski

(1961)), where the units of rad or Hz are not affected by 

raising them to a power. Or St(k) will have units of 

(m/sec)2

/(cycle/m) (Tank (1994)), as rad = cycle/2π, with the 

2π factor that multiplies the spectral density for this 

conversion absorbed into the constants αt terms. In this 

treatment, more convenient units in terms of Hz will be 

utilized as ((m/sec)

3

/Hz), by substituting cycles with Hz*sec,

which allows displaying results in terms of the acoustic wave 

velocity in (m/sec) /Hz.

Tatarski (1961) includes a more detailed treatment of ε

along with a detailed derivation of Equation (1). In brief, 

'2 2 * v l l ε = ν is called the eddy dissipation rate and it 

represents the energy dissipated as heat per unit mass per unit 

time due to a velocity fluctuation that occurs in an 

atmospheric region of size (length) l, due to a flow instability 

that takes place when a certain critical Reynolds number value

is exceeded (Tatarski (1961)). This critical Reynolds number 

cannot be determined precisely, but for worst case 

atmospheric turbulence calculations it would be correct to 

assume that this number is exceeded. The quantity ν is the 

kinematic viscosity, and ' vl is a velocity fluctuation in an

atmospheric region of size l.

At lower altitudes, the eddy dissipation rate can vary 

significantly with altitude and severity of turbulence (Fairall

(1991)) as well as globally (Nastrom (1985)). Tank (1994)

used a worst value of 8.6e–5 based on data collected at 

altitudes ranging from 25,000 to 40,000 ft (~ 7.5 to 12 km), 

which is about four times the appropriate value of ε for North 

Atlantic cruise altitudes. For supersonic cruise altitude, around 

60,000 ft (~18 km), the same value of ε can be used, as ε is 

relatively constant at altitudes above 20,000 ft (~6 km). 

According to Tatarski (1961), Equation (1) holds for 

turbulence lengths much greater than what would be 

considered a microstructure (10’s of meters) for eddy 

dissipation and much less than the long length scale (that is 

considered to be about 400 km). Some sparse measurements 

performed by Obukhov (1949) in the lower troposphere agree 

with Equation (1), as well as other measurements by Nastrom

(1985). The constant αt is constant for each type of 

disturbance, given by Tank 1994, to fit observed data:

αl = 0.15 (longitudinal wind velocity gust, 

dimensionless)

αv = 0.2 (vertical or horizontal wind velocity 

gust, dimensionless)

αT = 0.39 (temperature disturbance, K2

s

2

m–2)

αP = 0.0005(Po/To)

2 (pressure disturbance, Pa2

s

2

m–2), 

where the subscript “o” denotes total 

atmospheric quantities. 

For a flight vehicle encountering an atmospheric 

disturbance, the disturbance frequency is defined as

kMa Ma f = λ = (2)

Where M is the vehicle Mach number, a is the local speed of 

sound (m/sec), λ is the wave length of the atmospheric 

disturbance (m/cycle), and k is 1/λ (cycles/m). The 

disturbance frequency is based on an aero vehicle traveling at

a certain altitude and Mach number that encounters an 

atmospheric disturbance with a certain wavelength. There 

would be a corresponding frequency based on Equation (2) 

around which atmospheric disturbances will be exhibited. 

However, a control design would need to be able to 

sufficiently attenuate these disturbances under a wide range of 

flight and atmospheric conditions. Thus, accounting for the 

full operating envelope, a controls design would need to 

attenuate atmospheric disturbances in a wide range of 

frequencies. Normally, lower frequency atmospheric 

disturbances will be easier for the control design to handle.

The speed of sound is related with temperature as 
NASA/TM—2010-216961 3

a = γRTs (3)

where Ts is the static temperature in K.

According to the International Standard Atmosphere (ISA)

or (Anderson (2000)), the temperature decreases at a rate of 

approximately 6.5 K per km, from sea level up to 36,000 ft 

(~11 km). From 36,000 to 65,000 ft (~11 to 20 km) which 

would be the approximate cruise altitude of a supersonic 

transport, the temperature remains constant at about 216 K. 

Thus, temperature (K) as a function of altitude, h (km), can be 

expressed as

( ) 6.5( ) T h = Toh − h − hoh (4)

where Toh is the temperature at an arbitrary reference height 

hoh, which is accurate up to 36,000 ft and then stays constant 

above that. Based on that, at supersonic cruise conditions the 

speed of sound would be approximately 295 m/sec.

Figure 1 shows the zonal (east-west) winds, meridional 

(north-south) wind, and potential temperature power spectral 

as reported in Nastrom (1985). The units of the wind power 

spectral are in (m3

/sec2)/rad and temperature is in (K2m/rad). 

As shown in this figure, at the long length scale (i.e., beyond 

400 km) the power of turbulence decreases as the –3 power 

with decreasing wavelength (increasing frequency), thereafter 

decreasing as the –5/3 power which agrees with Equation (1).

Therefore, Equation (1) shows a good representation of the 

atmospheric spectral density with a wavelength smaller than 

400 km.

Figure 1.—Wind and Potential Temperature Spectra as 

reported by Nastrom (1985). Note: for clarity, the meridional 

wind and potential temperature spectra have been shifted 

one and two decades to the right, respectively.

The Tank Model of the Atmospheric 

Disturbance

As an alternative to the Kolmogorov spectrum, Tank 

(Soreide and Tank (1996) and (1997)) scaled the von Karman 

spectrum to fit the Kolmogorov model in the limit (i.e., for 

large k), which compares well with other data. For 

longitudinal disturbances this spectrum is

( ) [ ( ( ) ) ]5 / 6 2

2 / 3 5 / 3 ,

1 1.339 2

2 2.7

Lk

Sl VK k L

+ π

= ε (5)

For transverse disturbances this spectrum is

( ) ( ( ) )

[ ( ( ) ) ]

11/ 6 2

2

2/ 3 5/ 3 ,

1 1.339 2

1.339 2

3

8 1

2.7

Lk

Lk

Sv VK k L

+ π

+ π

= ε (6)

The main difference of this von Karman form (compared to 

Kolmogorov model) is that the disturbance spectral levels off at 

low frequencies. Another difference is that in the von Karman 

form of the Tank model, the integral length scale, L, is explicitly 

employed, Equations (5) to (6). The integral length scale is 

related to the outer scale length, Lo (the length of the 

atmospheric turbulence patch) by L=Lo/14.7 for the longitudinal 

disturbance, and L=2Lo/19 for the transverse disturbance 

(Soreide and Tank (1996)). Figure 2 shows a plot of the 

Kolmogorov acoustic wave velocity spectral in (m/sec/Hz)

compared to the corresponding von Karman spectral for the

longitudinal and transverse acoustic wave velocity disturbances. 

Different values of L will produce the same curves

Figure 2.—Acoustic wave velocity spectral comparisons for the 

Kolmogorov and von Karman spectral.
NASA/TM—2010-216961 4

Figure 3.—Longitudinal acoustic wave velocity spectral 

comparisons between the Kolmogorov and the von Karman 

forms of four different integral scale lengths.

Figure 4.—Longitudinal acoustic wave velocity spectral 

comparisons for different integral scale lengths.

with either a higher or lower, low-frequency asymptote as 

shown in Figure 3 for the longitudinal disturbance. For 

comparison, a value of L=762 m pertains to an atmospheric 

turbulence patch of approximately 11 km for the longitudinal 

disturbance, while a value of L=10 km pertains to a turbulence

patch of approximately 147 km long. Tank used a value of 

L=762 m, which he called as standard in the airplane industry. 

Since a typical control design can sufficiently attenuate 

disturbance frequencies well beyond 1 Hz, based on Figure 3

differences in the value of L and thereby, differences in the 

lower frequency asymptote, will have negligible effect on the 

control design. For completeness, the Kolmogorov spectrum 

for the longitudinal acoustic wave velocity is shown in 

Figure 4 for different values of eddy dissipation rates. As can 

be deduced from Figure 4, orders of magnitude difference in 

eddy dissipation rate doesn’t produce a corresponding

appreciable change in the Kolmogorov spectral density. For 

instance, approximately four orders difference in ε only makes 

approximately a factor of 8 difference in velocity as

( ) ( ) 1

2 / 9

1

2

2 S f S f 















ε

ε = , ((m/sec)/Hz) (7)

Where S1 is a known spectral density with an eddy dissipation 

rate ε1, and S2 is the calculated spectra density for a different 

eddy dissipation rate ε2.

Fractional Order Fit of Atmospheric 

Turbulence Model

Atmospheric turbulence, as shown in Figure 1, and as 

described by Kolmogorov and the Tank models in 

Equation (1) and Equations (5) and (6), is fractional order. In 

Appendix B a circuit analog of the Tank model von Karman 

form is utilized, which serves as the basis for deriving integer 

order pole-zero product TF approximations to the fractional 

order atmospheric disturbances. The reason for the need to 

derive approximations to the fractional order equations is 

because of the difficulty of explicitly or numerically solving 

fractional order differential equations. The reason is fractional 

order derivatives, unlike integer order derivatives, do not obey 

the locality law (i.e., the limit theorem), as further explained in 

Appendix B.

The idea behind the integer order TF approximation for a 

fractional order TF is as follows (see Fig. C.1): Starting at a 

frequency near the beginning of the equivalent 3 dB point of 

the fractional order TF, an integer order TF approximation can 

be developed symmetrically centered about the fractional 

order TF, like a descending staircase shape, by interleafing 

integer order poles and zeros. As the number of steps of this 

staircase TF approximation increase, the steps become shorter, 

eventually collapsing to the straight line of the fractional order 

TF as the number of poles and zeros of this approximation is 

increased to infinity. For this approach to work, the frequency 

of the poles and zeros need to be related to the atmospheric 

disturbance parameters and also be derived in a way such that 

the staircase TF approximation is symmetrically centered 

about the fractional order TF.

Based on the fractional order circuit analog of Equations (5) 

and (6), see Appendix B for detailed derivations, the values

for the equivalent capacitance and resistance of atmospheric 

turbulences are determined as 

(a ) ( Ma)

C x

t

t

ε π

=

2

1

1/ 2 / 3 (8)
NASA/TM—2010-216961 5

R ( )(a ) L

x

t t

1/ 2 / 3 = 1.3390 2π ε (9)

and the corresponding natural frequency is computed as

t t

n n R C

Kω ω = (10)

where the adjustment factor Kωn is 1, but represented in this 

form in case any adjustments need to be made to the 

atmospheric disturbance natural frequency. 

Utilizing these circuit analog parameters, a time domain 

disturbance can be derived (see Appendix C for detailed

derivations) based on an integer order TF approximation for 

atmospheric turbulence, given the atmospheric disturbance 

parameters ε, L, q, t, and r for the units conversion factor (with 

q=xr; x=5/3 and r=1/3 for acoustic disturbances and 1⁄2 for 

temperature and pressure) as

( )

( ) m t

pi

m

zi

t o t fit W

s

s

W K p

z

∏

∏

ω +

ω +

≅

1

1 , ,

1

1 )

(11)

where Wt is series of input sinusoids with unit amplitude 

frequency components distributed over the frequency range of 

interest, which is discussed more in Kopasakis (2010) for 

propulsion system type disturbances. The frequencies of the 

poles can be computed as

( )

( )

p i

j

Hpi zi j

i q

i

j

pi Hpi Hpi pi j

pi i m

K

2,3, ,

10 1 1

1

1

1

(2 1)

1

1 = 

ω ω + −

ω ω ω +

ω =

∏

∏

−

=

− η −

−

=

ω −

(12)

where the first pole is computed as follows

( ) q

q

q

p n

−

η ω = ω −

1

1 10 1 (13) 

and the frequencies of the zeros can be computed as

( )

( ) z i

j

Hzi pi iq

i

j

zi Hzi Hzi zi j

zi i m

K

1,2, ,

10 1 1

1

1

2

1

1 = 

ω ω + −

ω ω ω +

ω =

∏

∏

=

− η

−

=

ω −

(14)

For n number of frequency decades desired to be estimated, the 

number of poles and zeros in the TF approximation can be 

determined as

mp= (n – 1)/η (15)

mz = mp – 1 (16)

For a desirable pole-zero density pair, ρpz, per decade to be 

used to approximate the fractional order disturbance

ρpz

η = 2

1 (17)

and the terms ωHpi and ωHzi in Equations (11) and (12) can be 

computed as

( ) p

q q i Hpi n 10 1 , i 2, ,m

1/

ω = ω η (2 −1) − =  (18)

( ) z

q qi Hzi n 10 1 , i 1,2, ,m

1/

ω = ω 2η − =  (19)

The utility of ωHpi and ωHzi are to maintain symmetry, so that 

the staircase pole-zero approximation is symmetrically located 

on top of the fractional order TF asymptote. The proportional 

gains Kωpi and Kωzi in Equations (12) and (14) are reserved for 

final adjustments that may be needed to these frequencies for 

more closely approximating the fractional order TFs 

representing atmospheric disturbance.

Longitudinal and Transverse Acoustic Wave 

Turbulence

The longitudinal and transverse acoustic wave atmospheric 

disturbances are in the form of pure wind gusts in the axial 

direction of vehicle motion and in the vertical direction of 

motion respectively. Based on the fractional order fits 

determined in Appendix C, using Equations (8) to (19), with 

q=xr=5/9 (r=1/3), then Kt,fit in Equation (11), based on 

inspection of Equations (5) and (6), is

( )

1/ 3 2 / 3 5/ 3 Kl,fit = 5.4ε L (20)

for the longitudinal disturbance, and

( )

1/ 3 2 / 3 5/ 3 Kv,fit = 2.7ε L (21)

for the transverse disturbance. For n=3 (i.e., the TF fit over the 

span of 3 decades in frequency), the proportionality factor 

adjustments to improve these fits have been found to be

(Appendix C),

[ ]

[2.4; 1 1 1/ 2.4 1/1.5; 1 1 1]

; ; ,

=

Kl ω = Kωn Kωpi Kωzi (22)

for the longitudinal acoustic wave, and

[ ]

[4.27; 1 1 1/ 2.4 1/1.5; 1 1 1]

; ; ,

=

Kv ω = Kωn Kωpi Kωzi (23)
NASA/TM—2010-216961 6

for the transverse. These TF fits are compared against the 

Tank von Karman spectral of Equations (5) and (6), with the 

units conversion exponent r, as shown in Figures 5 and 6 for

ε=8.6e–5 (m2

/sec3) and L=762 m. As shown in these figures, 

the fits do a good job in approximating these disturbances. 

More accuracy can be achieved around the knee of these 

spectral by increasing the density of pole-zero pairs per 

frequency decade in this region. A certain inaccuracy at the 

lower frequencies would be acceptable for a feedback control 

design, since a typical control design should have no problem 

attenuating disturbances at this low frequency range. 

Figure 5.—Longitudinal von Karman spectral, final adjusted TF 

fit (ε = 8.6e-5 m2

/sec3, L = 762 m).

Figure 6.—Transverse von Karman spectral, final adjusted TF 

fit (ε = 8.6e-5 m2/sec3, L = 762 m).

Temperature Turbulence

Atmospheric temperature turbulence causes both temperature 

as well as acoustic velocity disturbances due to the change in 

the speed of sound. For the most part, temperature disturbances 

result in vertical displacement of air parcels (so called gravity 

waves). Therefore, acoustic velocity disturbances caused by 

temperature will generate vertical wind gusts that add with any 

transverse acoustic velocity gusts. However, for the propulsion 

system, according to Ashun (2004), the wing forebody will turn 

a vertical gust into a longitudinal gust, by multiplying the 

vertical gust by the conversion factor ( 1) 1 2 M∞ − M∞ − . This 

factor amounts to 0.63 for M∞ = 2.35, and can be ignored for 

worst case purposes. However, this conversion factor 

decreases with speed and can be taken into account, especially 

at lower supersonic speeds, for longitudinal acoustic gusts due 

to temperature turbulence. In Appendix B, Section B.2, first the 

Kolmogorov spectrum of the longitudinal acoustic wave and 

that of the temperature, based on Equation (1) are plotted

(Fig. B.6). This results in parallel spectral lines with frequency, 

similar to the Kolmogorov spectra shown in Figure 2. Then the 

von Karman temperature spectral is constructed by scaling the 

horizontal von Karman type acoustic wave spectral (i.e., scaling 

it by the difference in magnitudes of these Kolmogorov spectra

as shown in Fig. B.6).

Temperature turbulence spectral density, as can be seen in 

Figure 1, also follows the 5/3 law, but its units are in terms of 

Kelvin squared. Thus, in order to convert the units to Kelvin, 

the exponent r becomes 1⁄2, which makes the fractional order 

q=5/6. Based on this, the TF fit performed for temperature is

done the same way as with acoustic wave in the previous 

section (i.e., utilizing Equations (8) to (19), with the additional 

relations.

2 / 3 5/ 3 KT ,fit(temp) = 14.0ε L (24)

[ ]

[1.5; 1 1 1/1.1 1/1.2; 1 1 1]

; ; ,

=

KT ω = Kωn Kωpi Kωzi (25)

Figure 7 shows this TF fit and the spectral of Equation (5) 

(i.e., Eq. (5) substituting the scaled Eq. (24) for the numerator

and by also applying the units conversion factor in the 

denominator, for ε=8.6e-5 an L=762 m.

The proportionality factors, Kt, for all the fits can also be 

considered part of the formulations. Their accuracy at higher 

frequencies starts deviating somewhat with higher values of L

(integral length scales of few thousand meters and above),

values that are considered to be outside the typical range. As 

discussed before, L=762 m is considered standard in airplane 

industry. 

Temperature turbulence also causes acoustic wave gusts due 

to change of the local speed of sound. By perturbing the 

relation of the speed of sound and temperature in Equation (3), 
NASA/TM—2010-216961 7

Figure 7.—Temperature von Karman spectral and 

its TF fit (ε = 8.6e-5 m2/sec3, L = 762 m).

and substituting the resulting expression for ∆a into the Mach 

number equation of M=v/a, the relation between the change in 

temperature and acoustic velocity can be obtained as 

T

a

M R

v o

∆ γ ∆ = 2 (26) 

This relation can be applied to Equation (24) to compute

KT,fit(acoustic) for the acoustic wave velocity disturbance due to 

an atmospheric temperature fluctuation as

2 / 3 5 / 3 ,fit(acoustic) 14.0

2

L

a

M R K o

T ε γ = (27)

Equation (27) is applicable to transverse or vertical wind gust 

disturbances generated by temperature turbulence, and 

accurate in the worst case sense for longitudinal wind gusts, 

especially at higher supersonic Mach numbers. Generally, 

longitudinal wind gusts due to temperature variations can be 

expressed as follows by utilizing the conversion factor 

( 1) 1 2 M∞ − M∞ −

( ) 2 / 3 5 / 3

2 ,fit(longitudinal) 14.0

2 1

1

L

a M

M M R

K

o

T ε

−

− γ =

∞

∞ ∞ (28)

This TF fit for the acoustic velocity disturbance at supersonic 

cruise altitude, at Mach 2.35, for L=762 m and two different 

values of ε is shown in Figure 8. The TF fit for ε=8.6×10–5

(m2

/sec3) and for two different values of integral scale lengths 

is shown in Figure 9. As can be seen from Figure 8 or Figure 9 

compared to Figure 5, the acoustic velocity gusts produced by 

fluctuations in temperature are much higher in amplitude than 

those produced by pure acoustic velocity gusts.

Figure 8.—Von Karman acoustic wave velocity spectral due to 

temperature gust and its TF fits for different eddy dissipation 

rates (L = 762 m).

Figure 9.—Von Karman acoustic wave velocity spectral due to 

temperature gust and its TF fits for different integral scale 

lengths (ε = 8.6e-5 m2/sec3).

As shown in Appendix B, Figure B.8, for the combined 

acoustic wave, the acoustic wave velocity due to temperature 

gust dominates in the low frequency range and at higher 

frequencies both acoustic waves affect the total. At worst case, 

it can be assumed that both affects of the acoustic wave 

velocity can combine. It is also assumed that both the 

longitudinal and transverse acoustic wave velocities will not 

combine to produce worst case affects. For engines under the 

wing, the transverse wave is assumed to be converted to a 

longitudinal wave via the wing forebody (Ahsun (2004)), see 

Appendix B.
NASA/TM—2010-216961 8

Pressure Turbulence

The longitudinal velocity wave spectral and the pressure 

spectral of Equation (1) was utilized the same way as 

described in the previous section for the temperature to scale 

the von Karman type acoustic wave spectral to come up with 

the pressure spectral of atmospheric disturbances, see 

Section B.3. The units conversion exponent for pressure is the 

same as that for temperature (i.e., r=1/2), which makes the 

fractional order the same (i.e., q=5/6). As such, Kp,fit and Kp,ω

were computed as follows

2 / 3 5 / 3 KP,fit = 11.6ε L (29)

[ ]

[1.5; 1 1 1/1.1 1/1.2; 1 1 1]

; ; ,

=

KP ω = Kωn Kωpi Kωzi (30)

Figure 10 shows this TF fit for ε=8.6e–5 (m2

/sec3) and 

L=762 m compared to the scaled von Karman type spectral—

i.e., Equation (29) substituted for the numerator and 

proportional factor of Equation (5) and the denominator of 

Equation (5), raised to the power r. 

Density Turbulence

Density disturbances are contained within temperature and 

pressure fluctuations and its inclusion in a simulation, along 

with pressure and temperature will produce unnecessary 

additional disturbances. However, if need be, the relation of 

density disturbance with temperature and pressure can be 

derived as (by perturbing the state equation)

o

o

o

o

RT

P

R T T

P P − + ∆+ ∆ ∆ρ = ( ) (31)

A plot of density disturbance utilizing Equation (31) is 

shown in Figure 11.

Simplified Atmospheric Turbulence Model for 

Controls Design

The attempt in this section will be to further simplify the 

fractional order TF fit developed in the previous sections, by 

identifying certain characteristics of this model that are 

important for controls design to handle atmospheric 

turbulence.

Based on observations, utilizing Figure 8, differences in 

eddy dissipation rates, ε, produce spectral densities with the 

same frequencies, offset by a fixed magnitude. This offset is 

accounted by differences in TF gain. Calculating these 

frequencies using the formulations developed so far also 

Figure 10.—Pressure von Karman spectral and its TF fit 

(ε = 8.6e-5 m2

/sec3, L = 762 m).

Figure 11.—Density von Karman spectral and its TF fit 

(ε = 8.6e-5 m2

/sec3, L = 762 m).

supports this observation. Differences in the fractional order of 

the disturbances, such as longitudinal acoustic wave velocity

vs. temperature, obviously, will produce different spectral 

frequencies. Unlike differences in eddy dissipation rates that 

produce the same frequencies with different TF magnitudes,

inspection of Figure 9 shows that spectral with different 

integral length scales, L, will be represented by different 

frequencies. However, different integral length scales only 

affect the low frequency spectral, where a typical feedback 

controls design should have no difficulty attenuating 

disturbances at such a low frequency range.
NASA/TM—2010-216961 9

Figure 12 shows spectral density TF fits, first, with a fixed 

integral length scale and two different values of ε. This shows 

that the differences in ε produce a fixed offset in magnitude, 

but with the same pole-zero frequencies. Then Figure 12 also 

shows comparison between two spectral densities with the 

same value of ε and two different values of L. This later 

comparison, with different values of L, is represented with 

different pole-zero frequencies. However, their high frequency 

asymptotes match, which means that for the purpose of 

controls design, either of these TF fits can be utilized, despite 

their differences in pole-zero frequencies. 

Therefore, the frequency domain atmospheric turbulence 

model can be simplified by utilizing the model developed in 

the previous section to calculate the fixed frequencies due to a 

fixed integral length scale and variable eddy dissipation rates. 

Based on that, the simplified model developed for atmospheric 

turbulence, which ignores differences in the pole-zero 

frequencies due to different integral length scales can be 

represented as 

Figure 12.—von Karman longitudinal acoustic wave velocity 

spectral and its TF Fits due to different values of eddy 

dissipation rates and Integral Length Scales. 

( /1.46 1)( / 30.1 1)( / 85.7 1)( /1593.1 1)

( / 9.2 1)( / 55.0 1)( / 335.5 1) ( ) 70 2 / 9

+ + + +

+ + + = ε

s s s s

s s s GLA s (32)

( /1.46 1)( / 30.1 1)( / 85.7 1)( /1593.1 1)

( / 9.2 1)( / 55.0 1)( / 335.5 1) ( ) 56 2 / 9

+ + + +

+ + + = ε

s s s s

s s s GVA s (33)

( /1.1 1)( / 25.1 1)( /109.8 1)( /816.3 1)

( / 33.0 1)( / 45.6 1)( / 602.4 1) ( ) 943 2/ 6

+ + + +

+ + + = ε

s s s s

s s s GT s (34)

( /1.1 1)( / 25.1 1)( /109.8 1)( / 816.3 1)

( / 33.0 1)( / 45.6 1)( / 602.4 1)

1

( 1) ( ) 472

2

2 / 6

+ + + +

γ + + +

−

− = ε

s s s s

s s s

a

M R

M

M G s o

TLA (35)

( /1.1 1)( / 25.1 1)( /109.8 1)( / 816.3 1)

( / 33.0 1)( / 45.6 1)( / 602.4 1) ( ) 472 2 / 6

+ + + +

γ + + + = ε s s s ss s s

a

M R G s o

TVA (36)

( /1.1 1)( / 25.1 1)( /109.8 1)( / 816.3 1)

( / 33.0 1)( / 45.6 1)( / 602.4 1) ( ) 859 2 / 6

+ + + +

+ + + = ε

s s s s

s s s GP s (37)

Where GLA, GVA, GT, GTLA, GTVA and GP are the simplified 

atmospheric disturbance TF for longitudinal and transverse 

acoustic wave velocity, temperature, longitudinal and vertical 

acoustic velocities due to temperature, and pressure 

respectively. The factor ( 1) 1 2 M∞ − M∞ − in Equation (35) 

can be set to one for worst case turbulence conditions, 

especially for relatively high Mach numbers, in which case it 

makes Equation (35) and Equation (36) the same.

The TF fit described in the previous section was carried up 

to 200 Hz, which should make it accurate even beyond this 

frequency range. However, the model simplification 

represented by Equations (30) to (34) introduces additional 

error at high frequencies and its validity should be limited to 

about 200 Hz. The frequencies of this model are based on an 

integral length scale of L=762 m. This means that this 

simplified model is not unique, and other equivalent models 

can be calculated using different values of L with the 

equations developed in the previous sections. Thus, for a 

controls design, which can adequately attenuate 1 Hz 

disturbances, this model is accurate for values of L as low as 

100 m, which is well below the expected range of L. A value 

of L=762 m is considered standard in the airplane industry as 

mentioned before.

Time Domain Disturbances

The TF fits in the previous section can be thought as filters, 

which take input disturbances at desirable frequencies and 
NASA/TM—2010-216961 10

convert them to representative free stream atmospheric 

disturbances. The TF fits contain the representative 

magnitudes of these disturbances. Thus, the discrete time 

domain frequencies, Wt, only need to have unity RMS values. 

For instance, a time domain input can be constructed as a sum 

of unit amplitude sinusoids, starting with a low frequency,

somewhere at the low frequency asymptote of the spectral, 

and continuing up with sinusoids that are dispersed along the 

high frequency asymptote up to frequencies that covers the 

controller bandwidth. Also, the control design can be checked 

with single frequency sinusoids, sequentially simulated. Of 

course, the total energy or the maximum amplitude of the 

disturbance will need to be reasonable. For instance, it may be 

considered that total wind gust will not exceed an amplitude of 

180 mph, or some other reasonable number that depends on 

altitude. For more information on propulsion type disturbance 

considerations, see Kopasakis (2010).

For most of the spectral density simulations covered so far, 

the worst case eddy dissipation rate for low to moderate 

turbulence was used for altitudes of about 36,000 ft (~11 km) 

and above. However, the highest atmospheric turbulence ever 

recorded had an eddy dissipation rate of 1.7e-3 (m2

/sec3)

(Tank (1994)). Theoretically, it is possible that even more 

severe turbulence can be encountered according to McMinn 

(1997). McMinn provides tables for altitudes, turbulence 

severity, and associated eddy dissipation rates. Kopasakis

(2010) covers more in this area of simulating atmospheric 

disturbances, including considerations for atmospheric 

turbulence specifications. Fairall (1991) also covers

implementation of time domain sinusoidal input disturbances 

by utilizing Fourier sinusoids.

As mentioned before, previous works cover atmospheric 

turbulence probabilistically, as an exceedance for controls 

design purposes. In terms of probabilistic atmospheric 

turbulence models, there are no comparisons drawn here 

between those models and the non-probabilistic approach 

described in this paper, even though exceedance can also be 

incorporated with this approach. The modeling approach 

developed in this paper is geared towards implementing 

controls design by considering worst case atmospheric 

turbulence. This also allows implementation of control designs 

to cover worst case turbulence encounters, without penalizing 

performance under normal atmospheric conditions, by 

incorporating control scheduling in the design.

Figure 13 shows a time domain longitudinal acoustic wave 

disturbance utilizing the acoustic wave TF fit discussed in the 

previous section. The input to this TF is a combination of 

unity amplitude sinusoids distributed along the frequency 

range shown in Figure 5. For a single unity amplitude sinusoid 

at low frequency, the amplitude of the acoustic disturbance 

would be approximately 9 m/sec according to Figure 5. For 

this summation of unity amplitude input sinusoids, their

amplitudes combine to give approximately 15 m/sec 

Figure 13.—Longitudinal acoustic wave velocity disturbance 

created using the TF fit of the longitudinal acoustic wave 

with unity input sinusoids.

maximum disturbance. The same approach can be used to 

construct other atmospheric disturbances for time domain 

simulations.

In Kopasakis (2010), considerations for atmospheric 

turbulence specifications are provided, with more discussion 

on the frequencies of these turbulences, also involving the 

vehicle Mach number. In this reference a control design 

example is also provided for the supersonic inlet shock control 

problem, in the presence of atmospheric turbulence. 

Conclusion

In this paper derivations were carried out to approximate the 

fractional order atmospheric turbulence with integer order 

pole-zero transfer function fits for more accurate time domain 

simulations. An existing von Karman type model form of the 

Kolmogorov spectral is utilized initially for the acoustic 

velocity gusts disturbances. This model is scaled in this paper 

in order to also develop the von Karman model forms for other 

type of atmospheric disturbances like temperature, pressure, 

and density. Because of the fractional order of these models, a 

circuit equivalent is developed that is used as the basis to 

derive the integer order pole-zero approximations. Utilizing

the formulations presented in this paper, the transfer function 

poles and zeros of the approximations to the atmospheric 

disturbances can be directly computed for different parameters 

describing atmospheric turbulence. The model is derived to 

duplicate the fractional order form of atmospheric turbulence 

and is more accurate than previous models developed. These 

new models could be used to design controls for aerospace

vehicle propulsion or flight control systems.
NASA/TM—2010-216961 11

Appendix A.—Nomenclature

a speed of sound, (m/sec)

ao ambient speed of sound, (m/sec)

Ct equivalent electric circuit capacitance for t-type of 

disturbance, (farads)

f frequency, (Hz)

fn natural frequency, (Hz)

f c frequency of computed correction factor, (Hz)

h altitude, (km)

H vector of intersection points to match estimated and 

fractional order TFs 

k wavenumber, (rad/m) or (cycles/m)

K proportional gains

Kt von Karman horizontal asymptote for t-type of 

disturbance

Kωn adjustment factor to disturbance natural frequency

l length, (m)

L integral length scale, (m)

Lo outer scale length, (m)

mp number of poles in the TF approximation

mz number of zeros in the TF approximation

M mach number

n number of frequency decades desired to estimate the 

fractional order TF

P pressure, (Pa)

Po standard atmospheric pressure, (Pa)

q fractional order of equivalent electrical circuit

R universal gas constant, 287 (N*m)/(kg*K)

r units conversion exponent of atmospheric disturbance 

Rt equivalent electric circuit resistance for t-type of 

disturbance, (ohms)

S atmospheric turbulence spectral density 

Sl-or-v atmospheric velocity turbulence spectral density, 

(m3

/sec2)/(rad) or (m/sec)3/(Hz) also converted to 

(m/sec)/(Hz) by taking 1/3 root 

s Laplace operator

T temperature, (K)

To standard atmospheric temperature, (K)

Toh temperature at an arbitrary reference height, (K)

v flow velocity, (m/sec)

' vl velocity fluctuation in a region of size l of the basic 

laminar flow, (m/sec)

W unity disturbance time domain signal

Wt disturbance time domain signal for t-type of 

disturbance

x fractional order of atmospheric disturbance spectral

A.1 Greek

α constant associated with atmospheric turbulence 

spectral density

αl longitudinal, (unit less),

αv transverse or vertical, (unit less)

αT temperature, (°K2

sec2m-2)

αP pressure, (Pa

2

sec2m-2) 

γ ratio of specific heats, (γ = 1.41)

ε eddy dissipation rate, (m2

/sec3)

η the ratio of each decade interval where a pole or a 

zero will be used to estimate the fractional order TF

λ wave length of atmospheric disturbance (m/cycle)

ν kinematic viscosity of a laminar flow, (m2

/sec)

ρ weight density, (kg/m3

)

ρpz density of pole-zero pairs per decade for estimated 

TF

ωn natural frequency, (rad/sec)

ωpi frequency of TF pole i of integer order approx.

(rad/sec)

ωzi frequency of TF zero i of integer order approx.

(rad/sec)

A.2 Subscripts

a constant associated with circuit capacitance 

correctional factor

A variable associated with an adjustment

C constant associated with circuit capacitance 

correctional factor

H associated with symmetry frequencies for the approx. 

fit variable associated with TF fit or approximation

l variable associated with longitudinal atmospheric 

disturbance

LA variable associated with longitudinal acoustic 

disturbance

VA variable associated with vertical or transverse 

acoustic disturbance

P variable associated with pressure disturbance

t variable associated with type of disturbance

T variable associated with temperature disturbance

TLA variable associated with longitudinal acoustic 

disturbance due to temperature
NASA/TM—2010-216961 12

TVA variable associated with vertical or transverse 

acoustic disturbance due to temperature

v variable associated with vertical or transverse 

atmospheric disturbance

VK variable associated with the von Karman spectral 

density

VKA variable associated with the von Karman spectral 

density circuit approximation

o variable associated with an output

p variable associated with TF poles

s variable associated with a static quantity

z variable associated with TF zeros

ωpi variable associated with pole frequencies

ωzi variable associated with zero frequencies

∞ variable associated with freestream conditions
NASA/TM—2010-216961 13

Appendix B.—Circuit Analog of the von Karman Model Form

In this appendix the von Karman approximations to the 

different atmospheric disturbances based on the Kolmogorov 

spectrum given by (Soreide and Tank (a) and (b)) are

approximated by utilizing an electrical circuit analog. The 

purpose of the circuit analog is to serve as basis for deriving 

integer order pole-zero product TF approximations to the 

fractional order atmospheric disturbances. The reason for the 

need to derive approximations to the fractional order equations 

is because of the difficulty of numerically solving fractional 

order differential equations. For instance, unlike integer order, 

fractional order derivatives obey the law of non-locality. This 

property of fractional order derivatives makes the state 

transition matrix a convolution integral, with each integration 

step necessitating the computation of the convolution of the 

state all the way back to time zero where the state had an 

initial value of zero.

B.1 Longitudinal and Transverse Disturbances

The Kolmogorov atmospheric turbulence spectral model as 

described in the main body, Equation (1), of this paper is 

2 / 3 5 / 3 ( ) − St k = αtε k (B.1)

And the von Karman longitudinal and transverse models 

respectively, reproduced from Equations (5) and (6) from the 

main body of this report are:

( ) [ ( ( ) ) ]

5 / 6 2

2 / 3 5 / 3 ,

1 1.339 2

2 2.7

Lk

Sl VK k L

+ π

= ε (B.2)

( ) ( ( ) )

[ ( ( ) ) ]

11/ 6 2

2

2 / 3 5 / 3 ,

1 1.339 2

1.339 2

3

8 1

2.7

Lk

Lk

Sv VK k L

+ π

+ π

= ε (B.3)

The shape of the von Karman spectrum approximations of 

Equations (B.2) and (B.3) is that of a parallel network 

consisting of a resistive element corresponding to the horizontal 

von Karman asymptote and capacitive element corresponding to 

the high frequency Kolmogorov asymptote as 

( ) [ ( )]

(K ) S (k )

K S k S k r

t r t

r t t t VKA

+ , = , (m/sec)/Hz (B.4)

Where r represents the units conversion exponent from 

(m/sec)3

/Hz to (m/sec)/Hz (i.e., r=1/3 in this case—for the 

longitudinal and transverse acoustic wave disturbances).The 

variable St(k) signifies the Kolmogorov spectrum of 

Equation (B.1) or the high-frequency asymptote of the von 

Karman spectrum, and Kt, signifies the von Karman lowfrequency asymptote or TF gain for the type of disturbance of 

Equations (B.2) and (B.3) as 

2 / 3 5 / 3 Kl = 5.4ε L (B.5)

2 / 3 5 / 3 Kv = 2.7ε L (B.6)

Equation (B.4), for a flow, such as that of the acoustic velocity 

spectral, is equivalent to a circuit as shown in Figure B.1. The Kt

element is a gain for the type of disturbance and is a function of 

the parameters describing the atmospheric disturbance. The 

parameter Wt represents an input flow and Wt,o represents the 

resulting output flow, where the exponent q represents the 

fractional order of these elements. The circuit is an analog of 

this representation, acting as a disturbance source connected to a 

filter, meant to be connected to a sink or a system simulation in 

order to establish an output flow. Applying circuit theory, the 

TF of this circuit can be computed as

( ) t q

t t

r

t

t o W

R C s

K W

1 , + = (B.7)

Where Wt in the circuit stands for a unity flow signal 

representing the atmospheric disturbance, containing the 

disturbance frequencies of interest. More on that will be 

discussed later. The output signal, Wt,o, will be the time 

domain signal representing the von Karman spectral 

approximation. The subscript t stands for the type of 

disturbance as before, the exponent q represents the fractional 

order of the circuit, and the factor Kt is given by Equation (20) 

or (21). Where 1/(Cts)

q and Rt

q are the capacitive and resistive 

impedances of the circuit, with s representing the Laplace

operator.

The values of resistance Rt, capacitance Ct, and the 

fractional exponent q remain to be determined. The natural 

frequency of this circuit, ωn, can be directly obtained by

inspection of the TF of Equation (B.7). This frequency will 

coincide with the frequency at which the low-frequency

asymptote (Figs. 2 and 3) intersects the high-frequency

asymptote associated with the Kolmogorov spectral or the 

capacitive impedance of the circuit as

Figure B.1.—Equivalent circuit TF of atmospheric 

disturbance.
NASA/TM—2010-216961 14

t t

n n R C

Kω ω = (B.8)

For now Kωn is used as a placeholder in case any adjustments 

are needed to this natural frequency and its value here is set to 1.

In Equation (B.1), k stands for the wavenumber and it is

related to the frequency of the disturbance as

kMa Ma f = λ = (B.9)

For the circuit analogy, the natural frequency of the circuit 

should be set to the same as that for the von Karman spectral 

of Equations (B.2 and B.3). By inspection, the natural 

frequencies of both the von Karman longitudinal and 

transverse disturbances are the same and can be calculated by 

substituting Equation (B.9) for f into Equation (B.2) or (B.3) 

to obtain the relationship

1.339(2 ) . Ma

f L s = π

ω

(B.10)

Then by substituting the magnitude s for 2πf on the right hand 

side of Equation (B.10) and also setting ωn = 2πfn on the left 

hand side, we can solve for the natural frequency as

( )L

Ma fn π =

1.339 2 (B.11)

Additionally, the TF of Equation (B.7) needs to represent the 

von Karman approximation of the Kolomogorov model (i.e.,

the Kolmogorov asymptote). However, the Kolomogorov 

model, as a straight line decreasing with frequency in a log 

scale, can be represented in Figure B.1 with a capacitive 

impedance. Thus, Equation (B.1) must be this capacitive 

impedance. This capacitive impedance is raised to the r=1/3

power to account for the units change from (m3

/(sec2cycle)) to 

(m/sec)/Hz. Therefore, taking Equation (B.1) and raising it to 

the power r, replacing k with f/(Ma) using Equation (B.9), then 

multiplying both numerator and denominator by (2π)

(5/3)r

, and 

by setting the resulting equation to 1/(Cts)

q

, the following 

expression can be obtained.

( )

( )

( )

r

t

q t f

Ma

C s 















π

α ε π = 5/3

2/3 5/3

2

1 2 (B.12)

By inspection of Equation (B.12), since |s| = 2πf, then 

q= (5/3)r or in general q=xr (i.e., with x=5/3 in this case), then 

the capacitance, Ct, can be computed as 

(a ) ( Ma)

C x

t

t

ε π

=

2

1

1/ 2 / 3

(B.13)

Plugging Equations (B.11) and (B.13) into (B.8), with

ω = 2πf, Rt can be solved as 

R ( )(a ) L

x

t t

1/ 2 / 3 = 1.339 2π ε (B.14)

Neither the capacitance Ct nor the resistance Rt depends on the 

exponent r, utilized for units conversion, which could be 

indicative of the fundamental nature of these parameters. 

With this analysis, the circuit elements for the longitudinal 

and transverse von Karman model approximations, pertaining 

to Equation (B.7) have been computed. Figures B.2 and B.3

show the von Karman longitudinal and transverse spectral of

Equations (B.2) and (B.3) respectively and their circuit 

approximation of Equation (B.7). These approximations can 

be improved by increasing the magnitude of the circuit 

approximation (Kt in Equation (B.7) for the longitudinal wave 

in order to reduce the high frequency asymptote error at the 

lower frequency. 

To do this, an offset dB level is chosen (∆dB) so that the 

resultant horizontal asymptote of the circuit approximation is 

not too high above the corresponding horizontal asymptote of 

the von Karman spectral. This can be done by multiplying 

Equation (B.7) by a proportional gain Ka, corresponding to

this ∆dB level, as Ka=10∆dB /20. However, doing that will also 

raise the high frequency asymptote of the circuit 

approximation by the same dB level. To counter this affect,

the high frequency asymptote can be lowered by the same 

amount. Lowering the high frequency asymptote can be done 

by lowering the capacitive impedance, since, the capacitive 

impedance is inversely proportional to the capacitance, see the 

circuit in Figure B.1. Thus, this will also require multiplying 

the capacitance by the same proportional gain, Ka. 

For the transverse approximation, illustrated in Figure B.3, 

it is seen that the high frequency asymptotes of the circuit 

approximation runs in parallel with the von Karman Spectral. 

Therefore, the asymptote corresponding to the transverse

Figure B.2.—von Karman longitudinal spectral of Equation (B.2) 

and their respective circuit TF approximation of Equation (B.7), 

(ε = 8.6e-5 m2

/sec3, L = 762 m).
NASA/TM—2010-216961 15

Figure B.3.—von Karman transverse spectral of Equation (B.3) 

and their respective circuit TF approximation of Equation (B.7) 

(ε = 8.6e-5 m2

/sec3, L = 762 m).

approximation needs to first be raised by doing the reverse of 

what was done to the capacitance in the previous step (i.e., 

multiplying the capacitance by the reciprocal of a certain 

correction factor, Kc. Secondly, the magnitude of the 

transverse approximation can be adjusted upwards by a 

desired level, Ka, by performing the two steps described for 

the adjustment of the longitudinal approximation. These 

correction factors can be computed as follows.

(desired dBlevelmay range from 0 to 6 dB)

10 dB/ 20 ,

∆

= ∆ Kt a (B.15)

( ) ( ) l c v c v VK c v o c K K S f W f , , , , = 1, = − , (B.16)

where fc is the frequency where this correction factor is to be 

computed. With these factors computed, Equation (B.13) for 

the capacitance can be modified utilizing Equations (B.15) and

(B.16). But the approach would be to keep the values of the 

capacitance and resistance generic and instead incorporate 

these adjustment factors directly into Equation (B.7) as

t q

t t

q

t c

t a

t A t oA W

R C s

K

K

K

W

1

1/

,

,

, ,

+ 



































= (B.17)

where 

r Kt,A = Kt,aKt (B.18)

and r Kt (for the type of disturbance) comes from Equations

(B.5) and (B.6) raised to the power r.

Figure B.4.—Von Karman longitudinal acoustic wave spectral of 

Equation (B.2) and their circuit TF approximations with different 

magnitudes of adjustment (ε = 8.6e-5 m2/sec3, L = 762 m).

Figure B.5.—Von Karman transverse acoustic wave spectral of 

Equation (B.3) and their circuit TF approximations with different 

magnitudes of adjustment (ε=8.6e-5 m2

/sec3, L=762 m).

Figures B.4 and B.5 show the longitudinal and transverse 

acoustic wave spectral for the von Karman forms of Equations 

(B.2) and (B.3) and the equivalent circuit adjustments based on 

Equation (B.17) for different magnitude adjustments. The Kt,a

values for these adjustments are based on Equation (B.15). The 

Kv,c value used is 1.424. As discussed before, the important 

frequencies to match are the higher frequencies, above 1 Hz, 

because a typical control system design should have no problem 

attenuating disturbances at the lower frequencies. Therefore, the 

6 dB or even the 3 dB, adjustments would normally provide 

sufficient approximations to the von Karman spectral.
NASA/TM—2010-216961 16

B.2 Temperature Disturbance

Unlike the von Karman model forms of Equations (B.2) and 

(B.3) for the acoustic wave given by Soreide and Tank (1996), 

no comparable models exist for temperature, temperature 

generated acoustic wave, pressure, or density disturbances. 

Thus the approach here is to utilize the Kolmogorov spectral 

of Equation (B.1) for the acoustic wave disturbance and 

compute the differences from that and those of the 

Kolmogorov Spectral of Equation (B.1) for temperature. Then 

scale the respective von Karman models of Equations (B.2) 

and (B.3) by these differences to generate the von Karman 

models for temperature, pressure and density. 

Figure B.6 shows this scaling, with the resulting von 

Karman spectral for the temperature disturbance. This scaling 

can be done by graphically finding or by computing 

the difference in dB magnitude (at any frequency using 

Equation (B.1)) between the Kolmogorov temperature and the 

Kolmogorov longitudinal spectral, and then multiplying the 

von Karman longitudinal spectral by this scale factor 

(corresponding to the difference in dB magnitude) to come up 

with the von Karman temperature spectral. As a result, the von 

Karman temperature spectral is computed as 

( )

[ ( ( ) ) ] ((K *m sec) Hz)

,

1 1.339 2

2 7.0

2

5 / 6 2

2 / 3 5 / 3

,

Lk

L

ST VK k

+ π

= ε (B.19)

The horizontal asymptote for this spectral, based on 

Equation (B.19), is 

2/3 5/3 KT =14.0ε L (B.20)

The circuit approximation for this model is similar to that of 

Equation (B.17), except for some differences in the 

proportional gains and the fractional exponent as

T q

T T

q

T c

T a

T A

T oA W

R C s

K

K

K

W

1

1/

,

,

, ,

+ 



































= (B.21)

Where KT,a is the same as Equation (B.15) for the desired 

adjustment and 

r KT ,A KT ,aKT = , KT ,c = 1 (B.22)

As before, the fractional exponent is q=xr. However, this time 

for the temperature disturbance, the units conversion exponent 

is r=1/2 (see Fig. 1 for the temperature units). As such, the 

fractional exponent for temperature is q=5/6, instead of q=5/9 

for the acoustic wave disturbances. Figure B.7 shows a plot of 

the circuit approximation for temperature based on Equation

(B.21), with KT,a = 1.4125 (i.e., raising the circuit 

approximation magnitude by 3 dB) , compared to its scaled 

von Karman spectral of Equation (B.19).

Figure B.6.—Kolmogorov spectral of the longitudinal wave and 

temperature and scaling of the longitudinal von Karman 

spectral to come up with the von Karman temperature model 

form (ε = 8.6e-5 m2

/sec3, L = 762 m).

Figure B.7.—Kolmogorov and von Karman spectral for 

temperature and its von Karman circuit approximation 

(ε = 8.6e-5 m2

/sec3, L = 762 m).
NASA/TM—2010-216961 17

By perturbing the relation of the speed of sound and 

temperature, a = γRT , and substituting the resulting 

expression for ∆a into the Mach number equation of M=v/a, 

the relation between the change in temperature and acoustic 

velocity can be obtained as 

T

a

M R

v o∆ γ ∆ = 2 (B.23)

Where γ is the ratio of specific heats (~1.41), R is the universal 

gas constant (~287 – N*m/(kg*K)), and ao is the standard 

speed of sound at that altitude (~ 295 m/sec at supersonic 

cruise altitude). Figures B.8 and B.9 show the von Karman 

spectral for the longitudinal and transverse acoustic waves 

respectively, the temperature gust based on Equation (B.19)

turned to either longitudinal or vertical acoustic wave velocity 

gust, and the combination of the two. As shown in these 

figures, for the combined acoustic wave the acoustic wave 

velocity due to temperature gust dominates in the low 

frequency range. At higher frequencies, both affect the total 

and the resultant slope is more of a combination of the two 

disturbances. The question in encountering atmospheric 

turbulence is if it’s possible that the two effects (i.e., the 

acoustic wave due to pure wind gusts and acoustic wave 

disturbance due to a temperature fluctuation) can combine at 

the same time to produce worst case conditions. As for 

longitudinal and transverse disturbances combining, in terms 

of the propulsion system, the assumption is that the 

longitudinal disturbance would be the worst case (i.e., the 

longitudinal and the transverse or vertical would not combine 

at the same time). This is especially true for propulsion 

systems located under the wing, since, the vertical acoustic 

wave loses some velocity when it’s turned to a longitudinal 

acoustic gust via the vehicle wing forebody (Ahsun (2004)).

B.3 Pressure Disturbance

For an atmospheric pressure disturbance, as discussed in the 

previous section, the Kolmogorov spectral of Equation (B.1) 

will be used to come up with the von Karman model form by 

scaling the respective models of the longitudinal acoustic 

wave disturbance. Figure B.10 shows this scaling, with the 

resulting von Karman spectral for the pressure disturbance. As 

a result, the von Karman pressure spectral is computed as 

((Pa *m s ec) Hz)

, [1 (1.339(2 ) ) ]

2 5.8

( )

2

2 5/ 6

2 / 3 5/ 3

,

Lk

L

SP VK k

+ π

ε

=

(B.24)

The horizontal asymptote for this spectral, based on 

Equation (B.24), is 

2 / 3 5/ 3 KP =11.6ε L (B.25)

Similar to the Temperature spectral, the circuit 

approximation for this model is the same as that of Equation 

(B.17), except for some differences in the proportional gains 

as

P q

P P

q

P c

P a

P A

P oA W

R C s

K

K

K

W

1

1/

,

,

, ,

+ 



































= (B.26)

Figure B.8.—Von Karman longitudinal acoustic wave spectral, 

temperature disturbance turned to acoustic velocity, and the 

two combined.

Figure B.9.—Von Karman transverse acoustic wave spectral, 

temperature disturbance turned to acoustic velocity, and the 

two combined.
NASA/TM—2010-216961 18

Figure B.10.—Kolmogorov spectral of the longitudinal wave 

and pressure and scaling of the longitudinal von Karman 

spectral to come up with the von Karman pressure model 

form (ε = 8.6e-5 m2

/sec3, L = 762 m).

Figure B.11.—Kolmogorov and von Karman spectral for 

pressure and its von Karman circuit approximation, (ε = 8.6e-5 

m2/sec3, L = 762 m).

Where KP,a is the same as Equation (B.15) for the desired 

adjustment and 

r KP,A KP,aKP = , KP,c = 1 (B.27)

As with the temperature in the previous section, the fractional 

exponent is q=xr, with r=1/2. As such, the fractional exponent 

for pressure is the same as that for temperature (i.e., q=5/6,). 

With the value of KP,a = 1.4125. Figure B.11 shows a plot 

of the circuit approximation for pressure based on 

Equation (B.26), compared to its scaled von Karman spectral 

of Equation (B.24).

Figure B.12.—Kolmogorov and von Karman spectral for 

density and its von Karman circuit approximation (ε = 8.6e-5 

m2/sec3, L = 762 m).

B.4 Density Disturbance

The density disturbance won’t be needed when acoustic 

velocity, temperature and pressure disturbances are included 

in a simulation. However for completeness, density 

disturbance will be discussed here in brief. The magnitude of 

the density disturbance can be computed from the equation of 

state, P=ρRT, by perturbing this equation, which can be 

solved as

o

o

o

o

RT

P

R T T

P P − + ∆

+ ∆ ∆ρ = ( ) (B.28)

where Po and To are the standard atmospheric pressure and 

temperature in Pa and Kelvin respectively, which would be 

around 5500 Pa and 216 K at supersonic cruise. Figure B.12

shows a plot of the Kolmogorov, the Von Karman, and the 

circuit approximation spectral for density, all based on 

Equation (B.28) from values of temperature and pressure 

obtained in the previous sections. It can be seen from this 

figure that the Kolmogorov spectral obtained from Equation

(B.28), actually looks like a finite von Karman spectral at low 

frequencies. Also, its low frequency asymptote would be 

approximately –20 dB, which equates to a density of 

0.1 kg/m3

. From standard atmospheric tables, the atmospheric 

density at about 62,000 ft (i.e., at approximate cruise altitude) 

is also approximately 0.1 kg/m3

. Even though, at the limit as 

frequency approaches zero, the Kolmogorov spectral 

approximates the mean atmospheric density at cruise, the 

appropriate density disturbance would be that shown for the

von Karman or the circuit approximation.
NASA/TM—2010-216961 19

Appendix C.—Fractional Order TF Fits of Atmospheric Disturbances

Normally, with the derivation of circuit analogs and their 

respective TFs to describe atmospheric disturbances (as they 

were presented in the Appendix B), time domain simulations 

would be straightforward. However, these circuits and their 

TFs are fractional, which complicates performing time-domain 

simulations. The complication comes from a certain property 

of fractional order differential equations, which deals with 

non-locality. There are several works that deal with time 

domain solutions of fractional order differential equations

Schmidt (2006) and Lorenzo (2008), but even to date, solving 

these types of differential equations is still challenging. Unlike 

integer order differential equations, the solution space of 

fractional order derivatives is not local and it depends on the 

whole prior history of the state, resulting in a state transition 

matrix that is in the form of a convolution integral. Therefore, 

solving the convolution at each time instant, taking into 

account all the prior history of the state, can be both 

complicated and taxing on the computing resources. 

In this section the attempt will be to formulate an integer 

order approximation of the fractional order TF. Given the 

parameters describing an atmospheric disturbance, like eddy 

dissipation rate and integral length scale, the formulations will 

explicitly solve for the integer order TF poles and zeros that 

approximate the fractional order disturbance. This approach 

avoids the laborious process of hand fitting such an 

approximation, for every type of disturbance, and for every 

time an atmospheric disturbance parameter changes. After 

forming this approximation to the fractional order TF, and by 

selecting the frequencies of interest for the disturbance as an 

input to this TF approximation, the time domain atmospheric 

disturbance can be formed. 

If the desire is to estimate the fractional order TF at each 

frequency decade by a single first order pole-zero pair, then 

the estimation process for each frequency decade starts with a 

pole followed by a zero (i.e., ωpi < ωzi), with one last pole 

placed after the last decade in order to keep attenuating the TF

magnitude as

W i n s s

s W K t

pn

n

i pi

zi t o , 1,2, , 1

1

1

1

1 1

, =   ω +















ω +

ω +

≅ + =

∏ (C.1)

Where n is the number of pole-zero pairs of the estimated TF 

and K is the DC gain. So if the desire is to approximate the TF 

for three consecutive decades, with one pair of pole/zeros for 

each decade, starting at some desired frequency, then n would 

be three in this case. A sketch of this TF approximation as a 

staircase symmetrically located on top of the fractional order 

TF is shown in Figure C.1, with more description to follow 

later.

Thus, if one first order pole-zero pair is used to estimate 

each decade of the fractional order TF, the first pole should be 

placed at a certain frequency such that the pole gain dropping 

at a rate of –20 dB/decade, intersects the gain of the fractional

Figure C.1.—Pictorial diagram of fractional order TF fit.

order TF midway through, in the first decade. A fractional 

order TF gain would be decreasing at a rate q*20 dB/decade, 

where q is the fractional order. Therefore, midway through the 

first decade the initial gain is expected to decrease by a factor 

of (1⁄2)*q*20 dB as

K

s

K q

p s

(20)/ 20

2

1

1

10

1 1/ 2decade

−

= ω + = (C.2)

Divide by 20 is because of the 20log10(x) base scale. Solving 

for the frequency we get 

10 1 (1/ 2) 1 1/ 2decade − ω = = q p

s s (C.3)

Also, midway in the first decade, the fractional TF, based on 

its circuit approximation in Appendix B, would have the same 

gain as that of Equation (C.2) as 

K

R C s

K q

q t t s 2

1

10

( ) 1 1/ 2decade

−

= + = (C.4)

As was discussed in Appendix B, the natural frequency of the 

fractional order TF in Equation (C.4) is

t t

n R C

1

ω = (C.5)

Knowing ωn, Equation (C.4) can be solved for s to find out at 

what frequency to place the first pole in order for the gains of 

these two TFs to intersect midway through the first decade as

( ) q q n s

s 1/ (1/ 2) 10 1 1/ 2decade= ω − =

(C.6)

Substituting Equation (C.6) into Equation (C.3) for s, the 

frequency of the first pole can be solved as

1st order pole-zero TF 

approximation

freq.

amplitude, dB

ωp1ωz1 ωp2ωz2 ωp3ωz3 ωp4

Fractional order TF

d1

d2

ωHz1

ωHp2

ωHz2

ωHp3 ωHp4

ωHz3
NASA/TM—2010-216961 20

( ) q

q

q

p n

−

ω = ω −

1

(1/ 2) 1 10 1 (C.7)

So far in this development the assumption was that each 

decade of the fractional order TF will be approximated with 

one pole-zero pair as shown in Equation (C.1). Based on that, 

the estimated product of first order poles and zeros of 

Equation (C.1) will match the fractional order TF at half 

decade intervals. For that, the only distinction in 

Equation (C.7) is the 1/2 exponent in the right side of the 

equation. This successive ratio of desired decay matching, 

expressed here as H, can be generalized as

( ) ( ) q

q

H q

p n p

−

ω = ω −

1

1 10 1 1 (C.8)

where Hp1 is the decay ratio for the first pole, where the two 

TFs would match or the magnitudes will intercept. 

Following the placement of the first pole, the need is to find 

out at what frequency to place the first zero of this pair. At the 

end of the first decade, the gain of the fractional TF would 

have dropped by 10–q K as

( ) K

R C s

K

q s q t t− = +

10

1 at first decade (C.9)

By substituting 1/ωf, using Equation (C.5), for RtCt, the 

frequency s where the gain has dropped by the above 

magnitude can be calculated similarly as Equation (C.6)

( ) q q n s 1/ 10 1 at 1decade

= ω − (C.10)

Or in general terms, as a function of the second portion of the 

approximation, Hz1, desired to be matched (i.e., the first 

portion at 1⁄2 decade was matched with a pole, the second 

portion at 1 decade is matched with a zero (Hz1=1), and so on), 

this frequency can be expressed as

( ) q H q Hz n z 1/ ( ) ω 1 = ω 10 1 −1 (C.11)

where 

(decade) 1 atHz1 Hz s ω =

The purpose of the ωH frequencies is to ensure symmetry (i.e., 

to symmetrically locate the TF approximation on top of the 

fractional order TF. Thus, in order to accomplish symmetry, 

the desire is to also have the estimating TF with the first polezero pair to equate to the same gain at this frequency, ωHz1,

(i.e., the same gain as the right hand side of Equation (C.9), 

see also Figure C.1, as

( ) ( ) K

s

K s H q

p

z z s 1 at 1decade 10

1

1

1

1 − = ω +

ω + (C.12)

From this equation, the frequency of the first zero can be 

calculated as

( ) 10 1 ( 1 1 1) 1

1

1

ω ω + − −

ω

ω =

Hz p

q z

Hz z H (C.13)

Thus, choosing the ratio of decades at which it’s desirable to

match the TF estimate, Equations (C.2), (C.8), (C.11), and 

(C.13) can be used to compute the frequency of the first zero.

Following the same procedure to calculate the frequency of 

the second pole, by using Equation (C.9), but with the gain 

10–(3/2)q

K, (i.e., Hp2=1.5 decades), 

( ) ( ) q H q Hp n p 1/

ω 2 = ω 10 2 −1 (C.14)

Then using the TF with two poles and one zero equated to this 

gain, similar to Equation (C.12), the frequency of the second 

pole can be calculated as

( )

( ) 10 ( 1) 1

1

2 1

2 2 1

2 2 ω ω + −

ω ω ω +

ω =

Hp z H q

Hp Hp p

p p

(C.15)

The same procedure can be followed to calculate the rest of 

the frequencies of the poles and zeros of the estimating TF. In 

order to generalize these formulations (for approximating the 

fractional order TF), let’s define ρpz as the desired density of 

pole-zero pairs per frequency decade. In addition, let’s also 

define the number of equal decade subintervals (in log 

frequency scale) as η, where a pole or a zero will be used to 

closely match this TF as

ρpz

η = 2

1 (C.16)

Then, pole and zero vectors of decade intervals can be formed, 

where these TF will be closely matched as

Ηp = η [2(1) – 1 2(2) – 1...2(mp) – 1] (C.17)

Ηz = η[2(1) 2(2)...2(mz)] (C.18)

For n number of decades over which to estimate the fractional 

order TF, mp and mz can be computed as

mp = (n – 1)/η (C.19)

mz = mp – 1 (C.20)
NASA/TM—2010-216961 21

For instance, for a density of one pair of pole-zeros per decade 

(i.e., ρpz =1) and estimating the TF for three decades (i.e., 

n=3), these vectors would have the values Hpi = [1/2 3/2 5/2 

7/2], Hzi = [1 2 3].

Inspecting Equation (C.17), it can be seen that the first 

element of the vector, corresponding to the calculation of the 

first pole frequency, will always be equal to the value of η. 

Therefore, for convenience, Equation (C.8) can also be

expressed in terms of η.

( ) q

q

q

p n

−

η ω = ω −

1

1 10 1 (C.21)

Also, inspection of the equations that describe the 

frequencies of the zeros and poles (i.e., Eqs. (C.11), (C.13),

(C.14) and (C.15) respectively, it can be seen that a 

generalized relationship for these equations can be written as 

follows

( ) ( ) q H q Hpi n pi

1/

ω = ω 10 −1 (C.22)

A more generalized expression of Equation (C.22) can be 

arrived by substituting into it Hp from Equation (C.17) in 

terms of the pole number, i and η as 

p

i q q ωHpi = ωn (10η(2 −1) −1)

1/ i = 2,,m (C.23)

Then a generalized expression for the calculation of the pole

placement frequencies can be formulated as 

( )( ) ( )

[( )( ) ( )]

( )

( ) p i

j

Hpi zi j i q

i

j

pi Hpi Hpi pi j

Hpi zi Hpi zi Hpi z H q

pi Hpi Hpi pi Hpi pi Hpi pi

pi

i m

K

ω ω

K

pi

, 2, ,

10 1 1

1

10 1 1 .... 1 1

1 1 .... 1

1

1

(2 1)

1

1

1 2 1 ( )

1 2 1

= 

ω ω + −

ω ω ω +

=

+ ω ω + ω ω + −

ω ω ω + ω ω + ω ω + =

ω =

∏

∏

−

=

− η −

−

=

ω −

− −

ω − −

(C.24)

The proportionality factors Kωpi have been inserted in case any 

final adjustments are needed to these frequencies to improve 

the TF fit.

Similarly, for the frequencies for placing the zeros

( ) ( ) z

q H q Hzi n 10 zi 1 i 1,2, ,m 1/

ω = ω − =  (C.25)

Or by substituting in Equation (C.21) the Hz vector from 

Equation (C.18) in terms of the zero number, i, and η as

( ) z

q iq Hzi n 10 1 i 1,2, ,m 1/

ω = ω 2η − =  (C.26)

Then a generalized expression for the calculation of the zero

placement frequencies can be formulated as 

( )( ) ( )

[( )( ) ( )]

( )

( ) z i

j

Hzi pi iq

i

j

zi Hzi Hzi zi j

Hzi pi Hzi pi Hzi p H q

zi Hzi Hzi zi Hzi zi Hzi zi

zi

i m

K

K

zi

for 1,2,

10 1 1

1

1 1 .... 1 1 ( 10

1 1 .... 1

1

2

1

1

1 1 )

1 2 1

=

ω ω + −

ω ω ω +

=

ω ω + ω ω + ω ω + − −

ω ω ω + ω ω + ω ω + =

ω =

∏

∏

=

− η

−

=

ω −

−

ω − −

(C.27)

These formulations complete the computation of the pole-zero 

frequencies for the TF approximation of fractional order 

atmospheric disturbances. In terms of the indexes defined 

here, Equation (C.1) can be rewritten as

( )

( ) t m

pi

m

zi

t o t fit W

s

s

W K p

z

∏

∏

ω +

ω +

≅

1

1 , ,

1

1

(C.28)

As described in Appendix B, the von Karman spectral,

Equations (B.2) and (B.3), can be approximated with a 

fractional order circuit with a natural frequency that can be 

determined using Equation (B.8). Equation (B.8) is 

conveniently reproduced below in Equation (C.29).

t t

n n R C

Kω ω = (C.29)

with an equivalent capacitance and resistance for the type of 

disturbance t, derived as Equation (B.13) and (B.14)

(a ) ( Ma)

C x

t

t

ε π

=

2

1

1/ 2 / 3 (C.30)

R ( )(a ) L

x

t t

1 2 3 = 1.339 2π ε (C.31)

Based on these derivations, TF approximations of fractional 

order atmospheric turbulences can be developed (for the von 

Karman spectral, for any fractional order 0 < q < 1) by using 

Equation (C.28), together with supporting Equations (C.16), 

(C.19) to (C.21), (C.25) to (C.27), (C.23) to (C.24) and (C.26) 

to (C.27), (C.29) to (C.31). 

In a pictorial sense, by referring back to Figure C.1, the first 

order pole-zero pair TF approximation of a fractional order 

atmospheric turbulence can be viewed as a staircase 

symmetrically located on top of the fractional order TF. The 

utility of ωHpi and ωHzi are to maintain this symmetry. This is 

accomplished in the derivations above by formulating the 

symmetry frequencies, ωH’s, such that the length segments 

between intercepts (like d1 and d2) or alternatively the 

horizontal distances on either side of the intercepts are equal.

In the subsections that follow, more detail will be discussed 

about developing the approximations for each type of 
NASA/TM—2010-216961 22

disturbance. In addition, some final adjustments that need to 

be made to the pole-zero frequencies and the proportionality 

constants K will be presented. Without these final adjustments, 

the fits will resemble the example fit shown in Figure C.2, 

with the fractional order fit winding its way around the 

fractional order TF 

C.1 Longitudinal and Transverse Acoustic 

Wave Disturbances—TF Fits

By applying the formulations derived in Section C.1, TF fits 

for the longitudinal and transverse atmospheric disturbances 

can be constructed based on Equation (C.28). Figure C.3 

shows a plot of the longitudinal von Karman disturbance of 

Equation (B.2), it’s circuit TF approximation reproduced from 

Figure B.2, and the TF fit computed using the equations

derived previously in this Appendix for q=xr=5/9 (i.e., x=5/3, 

r=1/3), η=1/2, n= 3 (decades), and with atmospheric 

turbulence parameters ε=8.6e-5 (m2

/sec3) and L=762 m. 

The maximum error in this TF fit for the circuit 

approximation, using 4 poles and 3 zeros, over a span of 31⁄2

decades, amounts to approximately 1.5 dB. The error for the 

actual von Karman spectral is larger; about 4 dB at 0.1 Hz and 

decreases as the frequency increases. Assuming, that the 

control system design can sufficiently attenuate frequencies at 

this low frequency range (as it should) then this TF fit 

approximation can be acceptable for the actual von Karman 

spectral. For the TF fit shown in Figure C.3, which is 

an approximation to the longitudinal von Karman spectral 

of Equation (B.2) the pole and zero frequencies in 

Equation (C.28) are

[ ]

[3.82 22.92 312.38]

0.6 12.54 85.71 2062.0 ,

ω =

ω =

zi

pi (rad/sec) (C.32)

with Kl,fit in Equation (C.28) given by the numerator in 

Equation (5) raised to 1/3 power to account for the units 

conversion as

( )3

1

2 / 3 5 / 3 Kl,fit = 5.4ε L (C.33)

The preceding development to derive the equations to 

approximate a fractional order TF was based on a single 

fractional order TF, and a circuit approximation to the

atmospheric disturbance was used for this development. 

However, as indicated by Equation (B.2) and (B.3), the von 

Karman spectral is not that of a single order fractional TF,

especially in the lower frequency spectral. Therefore, for 

better accuracy, the derived equations can be adjusted to better 

fit the von Karman spectral. Examining Figure C.3, it can be 

seen that the TF fit can better approximate the von Karman

Figure C.2.—An example Fit without final adjustments to 

proportionality factors and frequencies.

Figure C.3.—Longitudinal von Karman spectral, circuit 

approximation, and TF fit (ε = 8.6e-5 m2/sec3, L = 762 m).

spectral if the natural frequency ωn, Equation (C.5), is 

increased, about proportionally to the difference in magnitudes

between the fractional order TF and the approximation TF, at 

that frequency. However, doing that will also increase the 

pole-zero frequencies of the TF fit based on the equations 

derived. Therefore, some equivalent adjustment to some of 

these frequencies will need to be made. Figures C.4 and C.5

show the TF fit for the longitudinal and transverse von 

Karman disturbances respectively, with the following 

adjustments to Equations (C.24), (C.29) and (C.27)
NASA/TM—2010-216961 23

[ ]

[2.4; 1 1 1/ 2.4 1/1.5; 1 1 1]

; ; / ,

=

Kl v ω = Kωn Kωpi Kωzi (C.34)

With Kl/v,ω symbolizing the proportionality constants for the 

frequency adjustments for both the longitudinal and transverse 

disturbances, which gives the frequency values of ( for ε=8.6e5 m2

/sec3 and L=762 m)

[ ]

[9.18 55.02 335.48]

1.46 30.10 85.71 1593.1 ,

ω =

ω =

zi

pi (C.35)

The proportionality factor, Kl,fit, in Equation (C.28) is given 

again by Equation (C.33). Similarly, Kv,fit is developed using 

Equation (B.6) as follows, but also adjusting by 3 dB which

Figure C.4.—Longitudinal von Karman spectral, adjusted TF fit 

(ε = 8.6e-5 m2

/sec3, L = 762 m).

Figure C.5.—Transverse von Karman spectral, adjusted TF fit 

(ε = 8.6e-5 m2

/sec3, L = 762 m).

will appropriately raise its magnitude in order to match the 

high frequency asymptote as

( )3

1

2 / 3 5/ 3 Kv, fit =1.4 2.7ε L (C.36)

As can be seen in Figure C.5, the TF fit of the transverse 

disturbance is not very accurate at the low frequency, which 

would normally be fine for control system design purposes as 

discussed before. But it may not be acceptable, for let’s say, 

wing loading. Thus, further adjustments could be made if 

desirable. For instance, instead of multiplying Kv,fit with the 

3 dB gain of 1.4, the natural frequency ωn of Equation (C.29) 

could be multiplied by this factor to raise its frequency. Then 

for final adjustments, the natural frequency can be also 

multiplied by the ratio of the difference of frequencies between 

the fit and the actual von Karman spectral obtained at some 

fixed amplitude, at some high frequency. This can be done 

either graphically or analytically. Performing these 

manipulations, a better fit of the transverse disturbance is shown 

in Figure C.6, with a final value of Kωn = 4.27. This will result 

in the following adjustments for the transverse disturbance

[ ]

[4.27; 1 1 1/ 2.4 1/1.5; 1 1 1]

; ; ,

=

Kv ω = Kωn Kωpi Kωzi (C.37)

which gives the frequency values of ( for ε=8.6e-5 m2

/sec3 and 

L=762 m)

[ ]

[16.33 97.92 597.07]

2.60 53.56 152.55 2835.3 ,

,

,

ω =

ω =

v zi

v pi (C.38)

with a Kv,fit expression as

( )3

1

2 / 3 5 / 3 Kv,fit = 2.7ε L (C.39)

Figure C.6.—Transverse von Karman spectral, final adjusted 

TF fit (ε = 8.6e-5 m2

/sec3, L = 762 m).
NASA/TM—2010-216961 24

C.2 Temperature Disturbance—TF Fit

As discussed before, the units conversion factor for 

temperature is r=1/2, which makes the fractional exponent 

value q=5/6 (x=5/3). Substituting this value for the fractional 

exponent, with the same values of η, n, and with the same 

atmospheric parameter values as those used for the 

longitudinal and transverse acoustic disturbances, and again 

using the formulations previously developed in this Appendix, 

the poles and zeros for the TF fit can be directly computed as 

before. The resulting TF fit with Kt,fit in Equation (C.28) given 

by Equation (B.20) raised to the 1⁄2 power as

( )2

1

2 / 3 5 / 3 KT = 14.0ε L (C.40)

This fit (not shown here) will end-up with the same low and 

high frequency asymptotes as with the von Karman form, but

at somewhat lower frequency. Similarly to the longitudinal 

disturbance, this can be corrected by multiplying Kωn in 

Equation (C.29) by 3/2, which turns out to be the ratio of the 

frequencies of these two spectra at some high frequency, for 

select fixed amplitude. Finally, by some minor adjustments to 

the high frequency poles, these spectra can be matched fairly 

close, as shown in Figure C.7. This will result in the following 

adjustments for the temperature disturbance

[ ]

[1.5; 1 1 1/1.1 1/1.2; 1 1 1]

; ; ,

=

KT ω = Kωn Kωpi Kωzi (C.41)

which gives the frequency values of (for ε=8.6e-5 m2

/sec3 and 

L=762 m)

[ ]

[33.04 45.64 602.36]

1.10 25.11 109.77 816.35 ,

,

,

ω =

ω =

T zi

T pi (C.42)

with a KT,fit, using Equation (B.20) as

( )

1/ 2 2 / 3 5 / 3 KT ,fit(temp) = 14.0ε L (C.43)

If the desire is to simulate the acoustic disturbance 

generated by a temperature gust, this conversion can be done 

by utilizing the perturbation relation of the speed of sound 

with temperature described in Appendix B (Eq. (B.23)) as 

T

a

M R

v o∆ γ ∆ = 2 (C.44)

Together with Equation (C.36) to adjust KT,fit as

2 / 3 5 / 3 ,fit(acoustic) 14.0

2

L

a

M R K o

T ε γ = (C.45)

Figure C.7.—Temperature von Karman spectral and its TF fit 

(ε=8.6e-5 m2/sec3, L=762 m).

Figure C.8.—Von Karman acoustic wave velocity spectral due to 

temperature gust and its TF fit (ε = 8.6e-5 m2/sec3, L = 762 m).

The acoustic velocity disturbance due to a temperature gust, at 

approximate supersonic cruise altitude is shown in Figure C.8.

Based on all the derivations carried out so far, the TF fits 

should work for different values of atmospheric parameters 

(like those of Eq. (C.45)). Figures C.9 and C.10 as compared 

to Figures 3 and 4 for their shape, show the TF fits with 

different values of eddy dissipation rates and integral scale 

lengths, by employing the same TF fit equations and without 

changing any other parameters, like the KT,ω values, 

Equation (C.41).
NASA/TM—2010-216961 25

Figure C.9.—Von Karman acoustic wave velocity spectral due 

to temperature gust and its TF fits for different integral scale 

lengths (ε = 8.6e-5 m2

/sec3).

Figure C.10.—Von Karman acoustic wave velocity spectral due 

to temperature gust and its TF fits for different eddy 

dissipation rates (L=762 m).

C.3 Pressure Disturbance—TF Fit

The units conversion factor for pressure is the same as that 

for temperature (i.e., r=1/2), which also makes the fractional 

exponent value q=5/6. Substituting this value for the fractional 

exponent, with the same values of η, n, and H as those used 

for the longitudinal, transverse and temperature disturbances, 

and again using Equations (C.16), (C.19) to (C.21), (C.25) to 

Figure C.11.—Pressure von Karman spectral and its TF fit 

(ε = 8.6e-5 m2

/sec3, L = 762 m).

(C.27), (C.23) to (C.24) and (C.26) to (C.27), (C.29) to (C.31), 

the poles and zeros for the TF fit for pressure can be computed 

as before. The resulting TF fit, with Kt,fit in Equation (C.28) 

given by Equation (B.25) raised to the 1⁄2 power, is shown in 

Figure C.11. It turns out, that adjustments to KP,ω are exactly

the same as those for KT,ω, Equation (C.41). For completeness,

these frequency adjustments and the KP,fit equation for 

pressure disturbances are listed below. 

[ ]

[1.5; 1 1 1/1.1 1/1.2; 1 1 1]

; ; ,

=

KP ω = Kωn Kωpi Kωzi (C.46)

which gives the frequency values of (for ε=8.6e-5 m2

/sec3 and 

L=762 m)

[ ]

[33.04 45.64 602.36]

1.10 25.11 109.77 816.35 ,

,

,

ω =

ω =

P zi

P pi (C.47)

with a KP,fit , using Equation (B.25) as

2 / 3 5 / 3 KP,fit = 11.6ε L (C.48)

C.4 Density Disturbance—TF Fit

As discussed in Section B.4, once the temperature and 

pressure TF fits have been calculated based on the 

formulations in the previous two sections, the density 

disturbance can be directly obtained by utilizing 

Equation (B.28). A plot of this density spectral is shown in

Figure C.12.
NASA/TM—2010-216961 26

Figure C.12.—Density von Karman spectral and its TF fit 

(ε = 8.6e-5 m2/sec3, L = 762 m).

References

Anderson, J.D., Jr., (2000): “Introduction to Flight,” Fourth Edition, 

McGraw-Hill Series in Aeronautics and Aerospace Engineering.

Ashun, U., (2004): “Dynamic Characterization and Active Control of 

Unstarts in a Near-Isentropic Supersonic Inlet,” Master’s Thesis, 

Massachusetts Institute of Technology, Feb. 2004.

Fairall, C.W., White, A.B. and Thompson, D.W., (1991): “A 

Stochastic Model of Gravity-Wave-Induced Clear-Air 

Turbulence,” Journal of Atmospheric Sciences, Vol. 48, No. 15, 

Aug. 1991. 

Hoblit, F.M., (1988): “Gust Loads on Aircraft: Concepts and 

Applications,” America Institute of Aeronautics & Astronautics, 

AIAA Education Series.

Houbolt, J.C., Steiner, R., and Pratt, K.G., (1964): “Dynamic 

Response of Airplanes to Atmospheric Turbulence Including 

Flight Data on Input Response,” NASA TR R–199, Jun. 1964. 

Kolmogorov, A.N., (1941(a)): “Dissipation of Energy in the Locally 

Isotropic Turbulence,” Comptes rendus (Doklady) de 

l'Académie des Sciences de l'U.R.S.S., 32, 16–18, 1941.

Kolmogorov, A.N., (1941(b)): “The Local Structure of 

Turbulence in Incompressible Viscous Fluid for Very Large 

Reynolds’ Numbers,” Comptes rendus (Doklady) de l'Académie 

des Sciences de l'U.R.S.S., 30, 301–305.

Kopasakis, G., (2010): “Modeling of Atmospheric Turbulence as

Disturbance for Control Design and Evaluation of High Speed 

Propulsion Systems,” GT2010–22851, ASME TurboExpo, Jun.

2010.

Lorenzo, C.F. and Hartley, T.T. (2008): “Initialization of FractionalOrder Operators and Fractional Differential Equations,” ASME 

Journal of Computational and Nonlinear Dynamics, Vol. 3, 

Issue 2, Apr. 2008.

McMinn, J.D., (1997): “Extension of a Kolmogorov Atmospheric 

Turbulence Model for Time-Based Simulation Implementation,” 

AIAA Guidance, Navigation and Control Conference, New 

Orleans, LA, Aug. 11–13, 1997, AIAA–97–3532.

Nastrom, G.D. and Gage, K.S. (1985): “A Climatology of 

Atmospheric Wavenumber Spectra of Wind and Temperature 

Observed by Commercial Aircraft,” Journal of Atmospheric 

Sciences, Vol. 42, No. 9, May 1985, pp. 950–960.

Obukhov, A. M., (1949): “Structure of the Temperature Field in 

Turbulent Flow,” Translation Izvestiin Akademii, Nauk 

S.S.S.R., Ser. Geogr. i Geofiz., Vol. 13, No. 1, pp. 58–69.

Schmidt, A. and Lothar, G., (2006): “On the numerical Evaluation of 

Fractional Derivatives in Multi-Degree-of-Freedom Systems,”

Elsevier North-Holland, Inc., Fractional Calculus Applications 

in Signals and Systems, Vol. 86, Issue 10, pp. 2592–2601, Oct.

2006.

Soreide, D.C. and Tank, W.G., (1996): “Proposed Model of the 

Atmosphere for the High Speed Civil Transport Program,” 

Program Report, March 1996.

Soreide, D.C. and Tank, W.G., (1997): “Validation of a statistical 

Atmospheric Model Using 747 Vertical Gust Exceedance Data,” 

AIAA–975579.

Tank, W.G. (1994): “Atmospheric Disturbance Environmental 

Definition,” NASA CR–195315, Feb. 1994.

Tank, W.G. and Gillis, J., (1996): “Atmospheric Disturbance Models 

for Linear and Nonlinear System Response Analysis,” AIAA 

34th Aerospace Sciences Meeting and Exhibit, Reno, NV, 

AIAA–1996–394, Jan. 1996.

Tatarski, V.I., (1961): “Wave Propagation in a Turbulent Medium,” 

McGraw-Hill.
REPORT DOCUMENTATION PAGE Form Approved 

OMB No. 0704-0188 

The public reporting 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 information, including suggestions for reducing this 

burden, to Department of Defense, Washington Headquarters Services, Directorate for Information Operations and Reports (0704-0188), 1215 Jefferson Davis Highway, Suite 1204, Arlington, VA 22202-4302. 

Respondents should be aware that notwithstanding any other provision of law, no person shall be subject to any penalty for failing to comply with a collection of information if it does not display a currently valid OMB 

control number. 

PLEASE DO NOT RETURN YOUR FORM TO THE ABOVE ADDRESS. 

1. REPORT DATE (DD-MM-YYYY)

01-12-2010 

2. REPORT TYPE

Technical Memorandum 

3. DATES COVERED (From - To)

4. TITLE AND SUBTITLE

Atmospheric Turbulence Modeling for Aero Vehicles: Fractional Order Fits 

5a. CONTRACT NUMBER

5b. GRANT NUMBER

5c. PROGRAM ELEMENT NUMBER

6. AUTHOR(S)

Kopasakis, George 

5d. PROJECT NUMBER

5e. TASK NUMBER

5f. WORK UNIT NUMBER

WBS 984754.02.07.03.20.02 

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

National Aeronautics and Space Administration 

John H. Glenn Research Center at Lewis Field 

Cleveland, Ohio 44135-3191 

8. PERFORMING ORGANIZATION

 REPORT NUMBER

E-17566 

9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESS(ES)

National Aeronautics and Space Administration 

Washington, DC 20546-0001 

10. SPONSORING/MONITOR'S

 ACRONYM(S)

NASA 

11. SPONSORING/MONITORING

 REPORT NUMBER

NASA/TM-2010-216961 

12. DISTRIBUTION/AVAILABILITY STATEMENT

Unclassified-Unlimited 

Subject Categories: 01, 02, 03, 05, 07, 08 and 15 

Available electronically at http://gltrs.grc.nasa.gov 

This publication is available from the NASA Center for AeroSpace Information, 443-757-5802

13. SUPPLEMENTARY NOTES

14. ABSTRACT

Atmospheric turbulence models are necessary for the design of both inlet/engine and flight controls, as well as for studying coupling 

between the propulsion and the vehicle structural dynamics for supersonic vehicles. Models based on the Kolmogorov spectrum have been 

previously utilized to model atmospheric turbulence. In this paper, a more accurate model is developed in its representative fractional order 

form, typical of atmospheric disturbances. This is accomplished by first scaling the Kolmogorov spectral to convert them into finite energy 

von Karman forms and then by deriving an explicit fractional circuit-filter type analog for this model. This circuit model is utilized to 

develop a generalized formulation in frequency domain to approximate the fractional order with the products of first order transfer functions, 

which enables accurate time domain simulations. The objective of this work is as follows. Given the parameters describing the conditions of 

atmospheric disturbances, and utilizing the derived formulations, directly compute the transfer function poles and zeros describing these 

disturbances for acoustic velocity, temperature, pressure, and density. Time domain simulations of representative atmospheric turbulence 

can then be developed by utilizing these computed transfer functions together with the disturbance frequencies of interest. 

15. SUBJECT TERMS

Atmospheric turbulence; Atmospheric disturbances; Flight controls; Propulsion controls; Flight loads; Propulsion disturbances; 

Aeroservoelastic disturbances; AeroPropulsoServoElasticity; Gust loads; Vehicle gust loads 

16. SECURITY CLASSIFICATION OF: 17. LIMITATION OF

 ABSTRACT

UU 

18. NUMBER

 OF 

 PAGES

32 

19a. NAME OF RESPONSIBLE PERSON

STI Help Desk (email:help@sti.nasa.gov) 

a. REPORT

U 

b. ABSTRACT

U 

c. THIS 

PAGE

U 

19b. TELEPHONE NUMBER (include area code)

443-757-5802 

Standard Form 298 (Rev. 8-98)

Prescribed by ANSI Std. Z39-18