George Kopasakis

Glenn Research Center, Cleveland, Ohio

Atmospheric Turbulence Modeling for Aero 

Vehicles: Fractional Order Fits

NASA/TM—2010-216961/REV1

July 2015

This Revised Copy, numbered as NASA/TM—2010-216961/REV1, July 2015, supersedes the previous 

version, NASA/TM—2010-216961, December 2010, in its entirety. 
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 Technical Report Server—Registered 

(NTRS Reg) and NASA Technical Report Server—

Public (NTRS) 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 counter-part 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 co-sponsored 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.

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 to help@sti.nasa.gov

• Fax your question to the NASA STI 

Information Desk at 757-864-6500

• Telephone the NASA STI Information Desk at

 757-864-9658

• Write to:

NASA STI Program

 Mail Stop 148

 NASA Langley Research Center

 Hampton, VA 23681-2199
George Kopasakis

Glenn Research Center, Cleveland, Ohio

Atmospheric Turbulence Modeling for Aero 

Vehicles: Fractional Order Fits

NASA/TM—2010-216961/REV1

July 2015

This Revised Copy, numbered as NASA/TM—2010-216961/REV1, July 2015, supersedes the previous 

version, NASA/TM—2010-216961, December 2010, in its entirety. 

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

Level of Review: This material has been technically reviewed by technical management. 

Revised Copy

This Revised Copy, numbered as NASA/TM—2010-216961/REV1, July 2015, supersedes the previous version,

NASA/TM—2010-216961, December 2010, in its entirety. 

Extensive corrections and revisions have been made to the content and additional text has been inserted into the report.

NASA STI Program

Mail Stop 148

NASA Langley Research Center

Hampton, VA 23681-2199

National Technical Information Service

5285 Port Royal Road

Springfi eld, VA 22161

703-605-6000

This report is available in electronic form at http://www.sti.nasa.gov/ and http://ntrs.nasa.gov/

This work was sponsored by the Fundamental Aeronautics Program 

at the NASA Glenn Research Center.
NASA/TM—2010-216961/REV1 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 itdifficult 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/REV1 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

s2m–2) 

P = 0.0005(Po/To)

2

 (pressure disturbance, hPa2s2m–2 – 

Tank 1994, multiply P by 100 to 

change units to Pa2

s2m–2) 

Note Tank 1994 does not define the symbols Po and To for 

pressure and temperature in the expression for P. McMinn 

1997, who is referencing Tank 1994, defines these symbols as 

standard atmospheric quantities. However, the assumption here 

is that Tank 1994 meant to define these quantities as freestream 

atmospheric conditions. The reasons for this assumption are as 

follows: a) because if these quantities were referring to standard 

atmospheric conditions, P would turn out to be a constant, just 

like the rest of the t constants (no need to display it in a variable 

form), and b) because Tank 1994 describes the atmospheric 

turbulence measurements taken for pressure as being freestream 

measurements, and as such, it would have not made much sense 

to convert these measurements into stagnation quantities, either. 

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 
NASA/TM—2010-216961/REV1 3 

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 

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( ) s oh oh T h  T  h  h (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 (northsouth) wind, and potential temperature power spectral as reported 

in Nastrom (1985). The units of the wind power spectral are in 

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. 

(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. 

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 

Figure 2.—Acoustic wave velocity spectral comparisons for the 

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

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 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 

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. 

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 does not 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. 
NASA/TM—2010-216961/REV1 5 

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) 

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. Note that the natural frequency, 

n, does not depend on the different type of disturbances or the 

eddy dissipation rate, as these factors cancel out when 

resistance and capacitance are multiplied together in 

Equation (10). 

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) 
NASA/TM—2010-216961/REV1 6 

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) 

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 

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). 

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. 

Potential Temperature Turbulence 

For the sake of simplicity, the terminology for temperature 

turbulence/disturbances, T, is utilized in this development to 

indicate potential temperature. Maintaining potential temperature 

in the formulations and plots versus actual temperature, also 

allows for direct comparisons between Kolmogorov (see Fig. 1), 

von Karman, and fractional order fit models. Potential 

temperature is defined as the temperature that a parcel of fluid at 

pressure, P, would acquire if adiabatically brought to a standard 

reference pressure Pst, usually 100,000 Pa. To convert to actual 

static temperature, Ts, the levels displayed in the equations and 

the graphs should be adjusted by the following equation: 

 ܶ௦ ൌܶቀ ௉

௉ೞ೟

ቁ

ோ/௖೛

 (24) 

 

where R is the universal gas constant and ܿ௣ is the specific heat 

capacity at constant pressure. R/cp = 0.286 for air. Notice, that at 

altitudes near sea level, potential and actual temperatures are 

about the same, while they increasingly deviate with altitude (eg. 

Ts ~ 0.5T at 50,000 feet). 

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 
NASA/TM—2010-216961/REV1 7 

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 (25) 

 

  1.5; 1 1 1/1.1 1/1.2; 1 1 1

; ; ,



KT   Kn Kpi Kzi (26) 

Figure 7 shows this TF fit and the spectral of Equation (5) (i.e., 

Equation (5) substituting the scaled Equation (25) for the 

numerator and by also applying the units conversion factor in 

the denominator, for =8.6e-5 an L=762 m. Note that the units 

of temperature include (m/sec)1/2, which is preserved in order to 

maintain consistency with the derivations. However, the actual 

units would not include the 1⁄2 power factor. Also, as discussed 

previously, temperature in this development refers to potential 

temperature. 

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), 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 (27) 

Figure 7.—Temperature von Karman spectral and 

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

This relation can be applied to Equation (25) 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    (28) 

Equation (28) 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 



  



  (29) 

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. 

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 effects 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. 

Figure 8.—Von Karman acoustic wave velocity spectral 

due to temperature gust and its TF fits for different 

eddy dissipation rates (L = 762 m). 
NASA/TM—2010-216961/REV1 8 

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). 

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). Utilizing Equation (B.25), with 

this unit conversion factor (same fractional order as 

temperature, q=5/6), Kp,fit can be expressed as 

2/ 3 5/ 3

,fit KP  5.4(P /l) L (30) 

and Kp, was computed as follows: 

 

  1.5; 1 1 1/1.1 1/1.2; 1 1 1

; ; ,



KP   Kn Kpi Kzi (31) 

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 (30) 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 

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). 

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       ( ) (32) 

A plot of density disturbance utilizing Equation (32) 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. 
NASA/TM—2010-216961/REV1 9 

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 

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. 

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 polezero frequencies. 

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. 

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 

( /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 (33) 

( /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 (34) 

( /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 (35) 

( /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 (36) 

( /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 s

s s s

a

M R G s o

TVA (37) 

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

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

   

    s s s s

s s s G s P P l    (38) 

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 temperature related TFs of Equations (35) to (37) are 

expressed in terms of potential temperature, in order to allow 

for direct comparisons for the simplified TF fits. However, in 

order to convert them to actual temperature and associated wind 
NASA/TM—2010-216961/REV1 10 

velocities, these TFs should be multiplied by the corresponding 

conversion factor of Equation (24) or ቀ ௉

௉ೞ೟

ቁ

ோ/௖೛

. 

The factor   1 1 2 M  M  in Equation (36) can be set to 

one for worst case turbulence conditions, especially for 

relatively high Mach numbers, in which case it makes Equation 

(36) and Equation (37) 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 (33) to (38) 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. 

The frequencies calculated in this simplified model are also 

dependent on the vehicle Mach number, see Equations (8) and 

(10). For this simplified model, the Mach number used is 2.35, 

which makes this model sufficiently accurate down to a Mach 

number of about 1.5. However, a simplified model can be 

calculated at any atmospheric Mach number, using the 

equations developed in this paper and the assumptions utilized. 

Time Domain Disturbances

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

which take input disturbances at desirable frequencies and 

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 maximum 

disturbance. The same approach can be used to construct other 

atmospheric disturbances for time domain simulations. 

Figure 13.—Longitudinal acoustic wave velocity disturbance 

created using the TF fit of the longitudinal acoustic wave 

with unity input sinusoids. 

 
NASA/TM—2010-216961/REV1 11 

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 polezero 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/REV1 13 

Appendix A.—Nomenclature

a speed of sound, (m/sec) 

ao ambient speed of sound, (m/sec) 

cp specific heat capacity at constant pressure, (J/(kg⋅K)) 

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 freestream static pressure, (Pa) 

Pst standard reference pressure, (100,000 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 Potential temperature, (K) 

To freestream static temperature, (K) 

Ts static 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, (Pa2

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 
NASA/TM—2010-216961/REV1 14 

TLA variable associated with longitudinal acoustic 

disturbance due to temperature 

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 

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/REV1 15 

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/REV1 16 

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/REV1 17 

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 from0 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/REV1 18 

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/REV1 19 

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. In this case, however, the constant P in 

Equation (1) or Equation (B.1) has a dependency on freestream 

(static) atmospheric conditions. Therefore, the scaling factor 

P/l will be explicitly used to scale Equation (B.2) to come up 

with the von Karma spectral, recognizing that the difference in 

magnitude between the Kolmogorov spectral of the acoustic 

longitudinal velocity and pressure is exactly equal to this factor. 

As such, the von Karman model form for pressure can be 

expressed as 

    Pa * m sec Hz

, [1 (1.339(2 ) ) ]

2 2.7( / )

( )

2

2 5 / 6

2 / 3 5 / 3

,

Lk

L

S k

P l

P VK



  





 (B.24) 

The horizontal asymptote for this spectral, based on 

Equation (B.24), is 

2 / 3 5 / 3 KP  5.4( P / l ) 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. 
NASA/TM—2010-216961/REV1 20 

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

temperature disturbance turned to acoustic velocity, and the 

two combined. 

Figure B.10.—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 ,, aPAP KKK P  , K ,cP  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.10 shows pressure 

comparisons of the Kolmogorov spectral of Equation (B.1), the 

respective scaled von Karman spectral of Equation (B.24) (both 

converted from Pa2

 to Pa), and the respective circuit 

approximation of Equation (B.26). 

Figure B.11.—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

TTR

PP    )( (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.11 

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

circuit approximation spectral for density, all based on Equation 

(B.29) 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/REV1 21 

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 nonlocality. 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 order TF 

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

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

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/REV1 22 

  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 polezero 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 pole-zero 

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/REV1 23 

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/REV1 24 

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 spectral if the 

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). 

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/REV1 25 

 

  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/REV1 26 

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/REV1 27 

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/REV1 28 

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-07-2015

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-1

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-REV1

12. DISTRIBUTION/AVAILABILITY STATEMENT

Unclassified-Unlimited

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

Available electronically at http://www.sti.nasa.gov

This publication is available from the NASA STI Program, 757-864-9658

13. SUPPLEMENTARY NOTES

This Revised Copy, NASA/TM-2010-216961-REV1, July 2015, supersedes the previous version, NASA/TM-2010-216961, December 

2010, in its entirety.

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

36

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)

757-864-9658

Standard Form 298 (Rev. 8-98)

Prescribed by ANSI Std. Z39-18