EPA-600/4-76-007
January 1976
Environmental Monitoring Series
SPECTRAL MODELING OF
ATMOSPHERIC FLOWS AND TURBULENT DIFFUSION
Environmental Sciences Research Laboratory
Office of Research and Development
U.S. Environmental Protection Agency
Research Triangle Park, North Carolina 27711
-------
RESEARCH REPORTING SERIES
Research reports of the Office of Research and Development,
U.S. Environmental Protection Agency, have been grouped into
five series. These five broad categories were established to
facilitate further development and application of environmental
technology. Elimination of traditional grouping was consciously
planned to foster technology transfer and a maximum interface in
related fields. The five series are:
1. Environmental Health Effects Research
2. Environmental Protection Technology
3. Ecological Research
4. Environmental Monitoring
5. Socioeconomic Environmental Studies
This report has been assigned to the ENVIRONMENTAL MONITORING series.
This series describes research conducted to develop new or improved
methods and instrumentation for the identification .and quantification
of environmental pollutants at the lowest conceivably significant
concentrations. It also includes studies to determine the ambient
concentrations of pollutants in the environment and/or the variance
of pollutants as a function of time or meteorological factors.
This document is available to the public through the National
Technical Information Service, Springfield, Virginia 22161.
-------
EPA-600/4-76-007
January 1976
SPECTRAL MODELING OF ATMOSPHERIC FLOWS
AND TURBULENT DIFFUSION
by
Arthur Bass, Steven A. Orszag
Flow Research, Inc. (N.E. Division)
1 Broadway
Cambridge, MA 02142
Contract No. 68-02-1297
Project Officer
Kenneth L. Calder
Meteorology and Assessment Division
Environmental Sciences Research Laboratory
Research Triangle Park, North Carolina 27711
U.S. ENVIRONMENTAL PROTECTION AGENCY
OFFICE OF RESEARCH AND DEVELOPMENT
ENVIRONMENTAL SCIENCES RESEARCH LABORATORY
RESEARCH TRIANGLE PARK, NORTH CAROLINA 27711
-------
DISCLAIMER
This report has been reviewed by the Environmental Sciences Research
Laboratory, U.S. Environmental Protection Agency, and approved for
publication. Approval does not signify that the contents necessarily
reflect the views and policies of the U.S. Environmental Protection
Agency, nor does mention of trade names or commercial products constitute
endorsement or recommendation for use.
11
-------
Preface
In order to provide a strong quantitative basis for air
quality management through selective emission controls, there
have been increasing demands in recent years for air quality
simulation models of considerable sophistication and complexity.
In these models a major component is that describing the atmos-
pheric transport and turbulent dispersion of the pollutants, and
many efforts are in progress to improve the specification of these
effects. For the turbulent dispersion this is normally done through
the use of diffusion coefficients (eddy diffusivities) or by more
sophisticated turbulence parameterizations that are developed
through various "closure" schemes for the governing equations of
flow. These closure approximations are, of course, necessary
because of the well-known mathematical intractability of the
general non-linear time-dependent Navier-stokes equations that
govern the turbulent flow of a viscous fluid. Closure approxi-
mations all involve the use, at some stage in their development,
of empirically based constants of parameters, that in some cases
may require ad hoc field experiments or observations.
During the past few years, and largely stemming from the
pioneering efforts of Professor Steven A. Orszag and his associates,
a major breakthrough has been made in developing powerful new mathe-
matical-computational approaches for the accurate numerical simu-
lation of turbulent flow. It has been shown that many aspects of
-iii-
-------
turbulent motion may be numerically simulated with impressive
accuracy, directly from the equations of motion, and without
the somewhat arbitrary closure approximations involved in normal
parameterizations. It was these remarkable results that stimulated
the present contract. It was hoped that it might provide the ini-
tial step towards providing a basic technique for atmospheric tur-
bulence and dispersion estimates by accurate non-empirical methods.
In due course such results might provide a set of diffusion data
against which various approximations might be tested.
The recent methods of Professor S. A. Orszag involve so-called
spectral techniques, along with sophisticated computer programming.
These techniques require expansion of the dependent variable into
a Fourier (or other orthogonal) series of smooth functions. The
calculation proceeds by performing certain of the mathematical
operations in Fourier (or other) space and others in physical con-
figuration space, according to the mathematical form of the opera-
tion and its relative efficiency. Transformations between the two
spaces are made by Fast Fourier Transform (FFT) methods and the
use of special computational techniques. In addition to having
very powerful applications to the problems of turbulent flow and
dispersion, the spectral technique has a wide potential for numer-
ical simulation of many atmospheric flows. The present report
attempts to provide an overview of such applications. It also com-
pares the accuracy and advantages of spectral techniques, with the
more conventional finite-difference numerical techniques for the
solution of partial differential equations. In some places a tu-
torial-type presentation results, from the attempt to provide ade-
-iv-
-------
quate references and to make the potential of spectral methods
more familiar to the developers of meteorological air quality
simulation models. In addition, it provides an account of an
initial research attempt to develop direct solutions, without
parameterization approximations, for turbulent dispersion of
an inert pollutant in a thermally stratified atmospheric layer
that is governed by the Boussinesq equation. Unfortunately,
this research failed to reach its desired goal during the rela-
tively short period of the contract, although in anticipation of
renewed efforts along similar lines in the future, a fairly
detailed account is desirable.
As an example of the power of spectral techniques, some new
results in turbulence dynamics are presented, concerning the rate of
relaxation of anisotropic flows to asymptotic isotropy.
Finally, general recommendations are made concerning the
potential utility of spectral methods in a number of specific
problem areas, and specific suggestions are given concerning the
prospective improvement of several existing air quality simulation
models, using these spectral methods.
Environmental Portection Agency Kenneth L. Calder
Research Triangle Park Project Officer
North Carolina
October, 1975
-v-
-------
CONTENTS
Page
Title Page i
Disclaimer Page ii
Preface iii
Contents vii
List of Figures ix
List of Tables x
Acknowledgements xi
1.0 Introduction 1
2.0 Overview of Numerical Techniques for Meteorological
Flow Simulations 3
2.1 Comparative Accuracy of Finite Difference and Finite
Spectral Methods H
3.0 Spectral Methods 25
3.1 Boundary Conditions and Geometries 31
3.2 Treatment of Boundary Layers 33
3.3 Pseudospectral Methods 39
4.0 Applications of Numerical Spectral Methods 45
5.0 Comparison of Methods of Spectral and Finite
Difference Methods for Two Examples 48
5.1 Comparison of Methods for the Diffusive Color
Problem 48
5.2 Comparison of Methods for Neuringer's Problem 65
6.0 Spectral Modeling of the Turbulent Diffusion of a
Passive Scalar 77
6.1 The Numerical Model 79
6.2 Physical Parameters 86
6.3 An Unsuccessful Experiment 88
-vi*-
-------
CONTENTS (CONT'D)
Page
7.0 Relaxation of Anisotropic Turbulence to Isotropy 107
8.0 Recommendations for Future Research 113
9.0 References 115
Appendix 1 The Data Content of Discrete-Grid
Representations 118
Appendix 2 Discrete Fourier Transform Representations 121
Appendix 3 Elimination of Aliasing in Discrete Real-
Space Multiplications 124
Appendix 4 The Fast Fourier Transform 132
Appendix 5 Construction of Initial Perturbation
Flow Field 138
-vii-
-------
LIST OF FIGURES
Number
5.1
5.2
5.3
5.4
5.5
5.6
5.7
5.8
5.9
5.10
5.11
6.1A
6. IB
6. 1C
6.2A
6.2B
6.2C
6. 3A
6.3B
6.3C
Initial Point-Source Distribution (t = 0) . . .
Finite-Difference Scheme, 32 x 32 point resolu-
Spectral Scheme, 32 x 32 point resolution,
Finite-Difference Scheme, 32 x 32 point resolu-
tion v = 05, t = 1/2 revolution
Spectral Scheme 32 x 32 point resolution,
v = .05, t ~ 1/2 revolution
Finite-Difference Scheme, 32 x 32 point resolu-
Spectral Scheme, 32 x 32 point resolution,
Finite-Difference Model, 32 x 32 point resolu-
Spectral Model, 32 x 32 point resolution,
v = .005, t = 1/2 revolution
Correct Neuringer FCN for y = 1 > 't = 0 . 5 . .
Neuringer Problem (y = 1.0, t = 0.5)
Analytic Value of F (x,y) /F (0 , 0)
U Field (y— z projection) ,t=0
U Field (y-z projection) , t = 10 sec ....
U Field (y-z projection) , t = 2 x 10 sec . .
U Field (x-z projection) , t = 0
U Field (x-z projection) , t = 10 sec ....
U Field (x-z projection) , t = 2 x 10 sec . .
U Field (x— y projection) ,t=0
U Field (x-y projection) , t = 10 sec . . - .
3
U Field (x-y projection) , t = 2 x 10 sec . .
Page
52
53
54
55
56
57
58
59
60
67
68
90
91
92
93
94
95
96
97
98
-viii-
-------
LIST OF FIGURES (CONT)
Number Page
6.4A V Field (y-z projection), t = 0 99
6.4B V Field (y-z projection), t = 10 sec ...... 100
6.4C V Field (y-z projection), t = 2 x 103 sec .... 101
6.5A W Field (y-z projection), t=0......... 102
6.5B W Field (y-z projection) e t = 10 sec 103
6.5C W Field (y-z projection), t = 2 x 10 sec .... 104
7.1 Time History of R (t) for E,(0) = 0 and En (0) =
1/2E2(0) .. .a. ............... 110
-ix-
-------
LIST OF TABLES
Number Page
2.1 Phase Errors for a = 0.1 13
5A Analytic Evaluation of F(x,y)/F(0,0) for y = 1.0,
t = 0.5 69
5B Analytic Value of F (x,y)/F (0,0), y = 4.0, t = 0.5 71
5C Absolute Errors in Spectral Run 72
5D Absolute Errors in 48*48 Point Finite Difference
Run 73
5E Absolute Errors in 96*96 Point Finite Difference
Run 74
-------
Acknowledgements
We acknowledge with thanks the constructive suggestions of
K. L. Calder, R. Eskridge, and K. L. Dermerjian, who have contri
buted greatly to the effectiveness of this report. We are
especially greatful to K. L. Calder for suggesting to us the
spectral simulation of Neuringer's problem.
-xi-
-------
1.0 Introduction
Advances in computer technology, together with important
achievements in the design of efficient numerical algorithms,
have begun to make practical the use of direct spectral methods
for numerical simulations of many fluid flow problems. In a
spectral method, the dependence of the flow variables upon one
or more of the spatial coordinates is represented by expansions
in a discrete series of smooth orthogonal functions (Fourier
series expansions being the most familiar). The governing
dynamic equations are then reformulated as equations relating
the spectral coefficients of the fields so analyzed. To permit
numerical computation, the expansions are truncated to some
finite upper mode. The resolution and accuracy properties of
spectral methods are very different from those of the more
familiar finite-difference methods, although both classes of
methods are related to discrete spatial grids.
Although rapidly growing in popularity, these spectral
methods are as yet relatively unfamiliar to large segments of
the meteorological communities. We have attempted, therefore,
to provide a brief review of some important finite spectral
methods for a variety of flow simulations relevant to air quality
modeling. The report discusses at some length the relative ad-
vantages and disadvantages of finite-difference and finite spec-
tral techniques, especially in regard to ease of implementation,
accuracy, and physical utility. Implicit in this description
is a clear understanding of the considerations which influence
the choice of a particular spectral approach, and determine its
range of utility. We have included in Appendix 1, a discussion
-------
of the resolution implications of discrete-grid methods, from
the point of view of both finite-difference and finite spectral
representations.
Next, a review is made of the current state of spectral
modeling in a variety of geophysical fluid flow problem areas
for which spectral methods have proved fruitful and/or hold
much prospective promise: (a) numerical weather prediction,
(b) turbulence modeling, (c) methods for advective-diffusive
simulations.
In the course of the current contract effort, we have
developed a spectral model for the direct simulation of fully
three-dimensional, non-linear, time-dependent, turbulent,
incompressible flows, and a corresponding computer program
for the solution of the Navier-Stokes equations using high-
resolution Fourier spectral methods. The model has been de-
signed to treat, within the Boussinesq approximation, the influ-
ence of differing stratification profiles upon the spectral evolu-
tion of turbulent flow fields. We have attempted, using this pro-
gram, to study the evolution of a passive scalar field in a
sheared, turbulent velocity field, but unforeseen limitations
in computer core and running time precluded the attainment of
physically realistic results. More successfully, the program
was used to look at an important problem in turbulence transport
theory, that is the rate at which an initially anistropic flow
will relax to isotropy; the relaxation time constant is an impor-
tant term in some turbulence closure models.
A two-dimensional version of these finite spectral methods
has been used for the accurate simulation of two advective-diffu-
-2-
-------
sive point-source problems, which may have useful .implications
for air quality modeling efforts in which discrete point sources
are to be represented. It is shown that, in the two cases,
finite spectral methods offer significant advantages, both with
regard to accuracy and computer storage requirements, compared
with finite-difference representations of the same problems.
Finally, recommendations are made concerning the practioc'il
uses to which such finite-spectral methods can be put, in the
present state of air quality modeling, and in the light of more
realistic assessment of available computational resources.
2.0 Overview of Numerical Techniques for Meteorological Flow
Simulations
The historical development of numerical methods for the
simulation of meteorological and other fluid flows can be traced
back to the early part of the century. Especially noteworthy is
the pioneering manual forecast effort of L.F. Richardson, during
the first World War. Although unsuccessful, this effort provided
later generations of numericists with much insight into the limi-
tations imposed by small-scale phenomena which are meteorologically
irrelevant, but whose implicit presence in a numerical model leads
to computational catastrophe. The impetus of technological pro-
gress during the second World War spurred greatly the development
of numerical methods for Shockwave and other supersonic model
calculations, and resulted in the first systematic examinations
of finite-difference numerical models for such calculations.
With the advent of high-speed calculational devices by the early
fifties, it became practical to attempt the first numerical weather
-3-
-------
predictions (p-jortof t , Von Neumann, Charney), using simple models
which isolated, by suitable scaling arguments, the meteorologi-
cally-relevant large scales of motion and eliminated the small-
scale structural behavior which had frustrated Richardson's
efforts thirty years before. The first such numerical forecasts
were performed by constructing finite-difference approximations
for the evolution of the two-dimensional inviscid barotropic
vorticity equation
(2.1) ( + v-V) (c + f) = 0
where V (x,t) •= u(x,t)l + v(x,t)j is the horizontal velocity
field, f is the Coriolis parameter, and r, = 1Z _ — is the
3x 9y
geostrophic vorticity. Considerably more ambitious simulations
were soon attempted in a variety of fluid-mechanical contexts,
beginning from the viscous , incompressible Navier-Stokes equa-
tions :
(2.2) (3- + V-V)V = - ¥1 + WV2V
3 1 p —
(2. 3a) V-V = 0
where P(x,t) is the pressure, p the (constant) density, and
v the kinematic viscosity. Using the incompressibility condi-
tion (2.3a) one obtains a Poisson equation for the pressure field,
(2.3b) yV2P = - V- (V-V)V .
Extensions of such models to slightly compressible flows (in the
Boussinesq sense) provide the needed coupling of flow kinematics
and thermodynamics.
It was soon recognized that the proper formulation of boundary
-4-
-------
conditions is among the most crucial difficulties encountered
in numerical flow simulations. Appropriate boundary conditions
may include (a) the specification of the velocity V at no-slip
<-\
boundaries; (b) the specification of V and — V_ at no-stress
., n 3 n T
boundaries (n and T are respectively, the normal and tangential
directions to the boundary); (c) the specification of boundary-
layer momentum, flux when such fluxes cannot be computed explicitly;
(d) thermal boundary conditions, (e) the choice of inflow and out-
flow boundary conditions for systems whose natural boundaries lie
outside the designated computational boundaries, (f) the incorpora-
tion of topographic influences and differences in land-sea surface
properties. The problems associated with imposition of boundary
conditions remain, by and large, the principal impediments to
truly accurate numerical flow evolution, and continue to be areas
of intensive investigation.
The most widely applicable and popular numerical methods
for flow simulations are based upon finite-difference approxi-
mations, in one or more geometries, to various versions or exten-
sions of (2.1 - 2.3); comprehensive reviews of finite-difference
methods have been provided by Fromm (1970), Kreiss and Oliger
(1973), Roache (1972), to name but a few. Difference methods
>
fall into several general classes. We may distinguish, first,
between those difference schemes developed for Eulerian, and
those developed for Lagrangian, versions of the equations of
motion. In the Lagrangian formulation, the coordinate system
is embedded in the moving fluid. For flows which are relatively
uncomplicated, and not subject to extreme distortions, such
methods can lead to accurate results, while accurately preserving
— 5—
-------
advection properties and real discontinuities. However, for suffi-
ciently disorganized flows, the Lagrangian mesh becomes badly dis-
torted, and the calculational accuracy quickly deteriorates.
Lagrangian methods have been especially popular in the simulation
of strongly compressible flows.
Eulerian systems, in which the coordinate grid is fixed,
are much less susceptible to violent computational breakdown as
the fluid undergoes large distortion. The penalty paid is the
corresponding loss of accuracy in the representation of advection,
diffusion, and flow discontinuities. We shall discuss some of
these difficulties shortly.
Difference methods may be also distinguished with respect
to the choice of variables which are computed prognostically,
and those obtained diagnostically- A variable is said to be
computed prognosfcically if the equation used to govern its
temporal evolution explicitly includes its time-dependence.
For example, the velocity field in (2.2) is a prognostic variable.
A variable is said to be diagnostic when the equation used to
govern its temporal evolution makes no explicit reference to its
time dependence. For example, the pressure field in (2.3) is a
diagnostic variable, since its change with time is entirely
implicit in the time variation of the velocity field. These
choices are far from trivial to make and govern the physical
realism which will inhere to the computation. For example, in
two-dimensional versions of (2.1 - 2.3) one may solve directly
for the variables V and P, or alternatively, convert to the
vorticity-stream function formulation
-6-
-------
(2.4) + J(4»,C) = W2?
where i|> is the stream function defined by
(.2.5b) c = V
-------
However, I ho vorticity-streamfunction idea is still very useful
for the determination of boundary conditions in primitive vari-
able simulations (Orszag and Israeli 1974).
The ability of an inviscid numerical sgheme to conserve
energy and enstrophy (square vorticity) is an important consider-
ation from the point of view of physical realism, and long-term
computational stability. Let us consider briefly how well such
properties are preserved by finite-difference schemes. If the
Jacobian in (2.4) is evaluated using second-order centered differ-
ences on a uniform grid the resulting scheme violates energy and
enstrophy conservation. Therefore, Arakawa (1966, 1970) proposed
that the Jacobian be approximated by second-order centered differ-
ences when it is first written in the equivalent form
(2.6) J = ^[(CT|I ) - (Sty ) + \\) l, - 4> 5 + (\K ) - (iK ) ]
3 Yx y y x rx y y x y x x y
If (2.6) is approximated by second-order centered differences,
the discrete analogs of energy and enstrophy conservation
(2. 7) //*J dx dy = 0; {/ C J dx dy = 0
respectively, are preserved. Arakawa1s energy-conserving methods
have been very popular in recent years, apparently because it has
been believed that energy conservation is necessary to avoid non-
linear computational instability (Phillips 1959), usually called
"aliasing" instability. However, energy-conserving finite-differ-
ence schemes may not be as useful as expected. The Arakawa schemes
"semi-conserve" energy; that is, they conserve energy only in the
absence of time-differencing errors, but Kreiss & Oliger (1972)
have given examples of schemes that semi-conserve energy but
that are unstable when leapfrog time-differencing is used. in
-------
general, though, Arakawa-type schemes are less susceptible to
numerical instability than straightforward centered schemes
(GrammeItvedt 1969). They are, however, considerably less effi-
cient and yet not noticeably more accurate. It is now understood
that aliasing instability will appear only in simulations with
inadequate grid resolution) (Orszag 1971a, 1972, Fornberg 1972,
Kreiss & Oliger 1973), and the deliberate use of an energy-con-
serving scheme to avoid aliasing instability may lead to a grossly
inaccurate simulation.
Among the most attractive features of the spectral methods
to be discussed is the ease with which in the inviscid case
energy and enstrophy are preserved in the absense of time-differencing
errors. This is mainly due to the accuracy with which phase and
amplitude information is retained using a spectral form of the
Jacobian term.
An essential part of the vorticity-streamfunction method is
the inversion of the Poisson equation (2.5b). In a finite differ-
ence scheme, the inversion is performed by a standard relaxation
method. In a spectral scheme, the inversion is trivial, and
direct spectral methods, rather than relaxation methods, should
be used whenever possible (Hockney 1970, Dorr 1970, Swartztrauber
1972).
Primitive-variable methods are increasingly widespread,
especially in three-dimensional NWP simulations. In the staggered-
mesh second-order scheme (Fromm 1963, Lilly 1965, Harlow and Welch
1965, Orszag 1969, Williams 1969, Deardorff 1970), velocities are
defined at cell boundaries and pressures at cell centers. It has
the attractive features of efficiency, small truncation error and
-9-
-------
energy semiconservation. The Poisson equation for the pressure
is most conveniently solved by direct spectral methods. 1-Iarlow
and Welch (1965) and Piacsek and Williams (1970) have developed
techniques to ensure compliance with the incompressibility con-
dition; these are especially useful when relaxation methods are
used to solve (2.2, 2.3).
-10-
-------
2 . 1 Comparative Accuracy of Finite Difference and Finite
Spectral Methods
In this section we compare finite difference and finite
spectral methods with regard to several important classes of
numerical error. These may be identified as follows:
(2.1.1) First-Differencing Phase Errors
Consider the usual approximation to the first derivative
3A/9x' where A is a scalar function of the spatial variable
x :
(2.8
2Ax
Suppose we are attempting to represent in this manner the first
i KX
derivative of the sinusoid function A(x) = e . Then the
ratio of the approximation (2.8) to the exact result iK A(x)
is given by sin (KAx) /KAx- This ratio approaches unity
only in the limit KAx->-0. Note, that for large K and finite
Ax, the approximation (2.8) becomes progressivly poorer. In
particular for K-HT/AX , i.e., as we approach the limit of spec-
tral resolution of the finite mesh, the ratio approaches zero,
which is hardly acceptable if there is significant information
contained at wavenumbers near the truncation limit N/2. Con-
sider now the effect of such first-differencing errors upon the
solution of the one-dimensional advection equation (N is the
number of mesh intervals in the grid)
(2.9) -- A(x,t) + U A(x,t) = 0
-11-
-------
where (J is a constant, uniform advecting velocity. The solu-
tions to (2.9) can be synthesized by linear combinations of
waves of the form
(2.10) A (x,t) = eiKn(x-Ut).
n
In fact, if the values of A(x,t) are given on N discrete
points equally spaced by Ax, an arbitrary excitation on these
N points may be approximately represented as a linear combina-
tion of the N solutions (2.10) with KAx - 2rrm/N (-^ f_ m< - ) ,
m integral.
The leapfrog centered "2Ax difference" approximation to
(2.9) is
(2.11) _
where A. = A(jAx, nAt) , At is the time increment, and
a= UAt/Ax. The finite-difference equation (2.11) is satisfied
by
(2.12) An = ei(KjAX - n9)
D
where sin0 =otsin (KAx). If |a| < 1 then 6 is real for all
K and hence the scheme is neutrally stable. In each time step,
the exact change of phase of A(x,t) given by (2.10) is -KUAt
= -KaAx while the change of phase of (2.12) = -0 = - sin~
[asin KAx]. Since |9/(KaAx)|
-------
Let us compare the ratio of
-------
The values of 6 listed in Table 2.1 show that the fourth-
order scheme has smaller phase errors than the second-order
scheme, and that the truncated Fourier scheme has phase errors
of less than 1% for all waves satisfying | KAx | <-y, when
It should be noticed that the second-order scheme has lagging
phase error (6<1) , while the fourth-order scheme has mostly lagging
errors. [There is a region of leading phase error 6>1 for very
small K Ax for the fourth-order scheme] . The lagging phase
errors provide an explanation for the wakes of bad numbers given
by finite-difference approximations (see also the discussion in
Chapter 6) .
The truncated Fourier scheme has only leading phase errors
for finite At, and even these disappear in the limit At^-0 ,
if KAx0
(if the exact solution that is being approximated is infinitely
dif ferentiable) . Finite-difference approximations have non-zero
phase errors in the limit At+0.
2.1.2 Second-differencing phase errors
These errors are introduced in evaluating the Laplacian
2
operator V by the standard finite-difference operation
(2.13) (v2A)x ,. A(x£+i'V + A(*^ym) + A(x£--ym-i} + ^
X£'Ym ~ ----- ' - - - ' - 2
Ax
Consider the two-dimensional sinusoid A(x,y) = e ix 2
Then the ratio of the approximate expression (2.13) to the exact
-14-
-------
result -(kj + k22 ) is just
• 2 klA* . . 2
(2.14) 4 __fL: : ^ '
Ax2 2 , . 2
kl + k2
which approaches unity only in the limit k Ax-»0 and k Ax->0.
Note again that for large k] or k,, the second-order finite
difference scheme seriously misrepresents the Laplacian operator
2
as, for example, in the evaluation of the viscous term vV V.
Reduction of these errors is possible by use of a fourth-order
0
finite-difference approximation to V V, due to the higher order
of approximation, but little benefit is gained.
There are no second-differencing errors in the truncated-
2
Fourier expansion method, as the transform of vV V is exactly
p O
-vk tJ(fi). The algebraic form -vk U(k) makes it a simple matter
to treat viscous dissipation implicitly in time in the Fourier
method. However, this advantage is slight, because in fully
turbulent flows, it is the advective terms which are most impor-
tant, and the numerical simulations are thereby usually limited
by advective rather than diffusive stability criteria.
It should be noted here that, for systems where the pressure
field (or the streamfunction field) is recovered diagnostically
by solution of a Poisson equation (cf. (2.3b) or (2.5b), for
example), the second-differencing phase errors associated with
the representations of the Laplacian of these fields are consider-
ably more detrimental to the evolution of the system than those
introduced into the viscous terms.
-15-
-------
2.1.3 I _n compressibility Errors
The imposition of the imcoinpressi bill ty condition is ,1
trouble-some problem in many f inite-di f feronce schemes and com-
putationally expensive iterative procedures are often required
to maintain incompressibility to within a specified precision.
Incompressibility can be imposed exactly in a truncated Fourier
scheme by the following simple procedure.
Suppose that the three-dimensional velocity field is repre-
sented in Fourier space by U.(k_) , where U(k) is the three-dimen-
sional Fourier transform of the real-space velocity field V(x) .
The divergence of U(k) is given in Fourier space by ik-U(k_).
Let .us subtract out vectorially the divergent component of U(k)
by replacing U(k) with
(2.5) U' (k) = U(k) - -(-'-(k)) .
k-k
Note that if we now take the divergence of this entire expression
(that is, by performing the operation ik- upon both terms, we
get the identity
ik.u-(k) = ik-u(k) - JQi-kXk'lHiO) = o.
(k-k)
Hence, incompressibility is maintained exactly if, at the con-
clusion of each time step, we replace U(k_) by U' (k) . We have
thus removed the divergent component of the velocity field. (Note
how trivially such a constraint can be imposed in the Fourier
spectral scheme). What is truly remarkable about this method
is that, if the Navier-Stokes equations are written in rotational
-16-
-------
form, i.e.
1 i~\
(2.16) -— V(x,t) - V(x,t)* u>(x,t)- V(P + y V-V) + vV V(x,t)
dt — — — — — ^ —
where u = VxV is the vorticity, then this system may be evolved
in time without any reference to the pressure head term, V (P+^- V-V)
so long as incompressibility is imposed at each time step in the
manner indicated above. This is because, in an incompressible
system, the pressure field serves only to enforce incompressibility,
by dynamic readjustment of the velocity fields.
2.1.4 Coupling errors
Among the crucial difficulties with finite-difference approaches
are the coupling errors introduced when non-linear differential opera-
tors are represented. Consider again for example the simple two-
dimensional barotropic vorticity equation;
(2.17) ^L v2(/;= _v-v (V2^)
^\
V = k x Vip
where if; is the geostrophic streamfunction. If we represent the
Fourier transform of ^ (x) by 0 (k) as defined in Appendix 2,
then it can be shown (Lilly 1965) that the wave-space version
of (2.17) is just
(2.18) -^ {k.k<|>(k)} = II*1 x *_" (k"-k" - k1 -k1 KMk" )$ (kj1) dk
k'+k"=k .
For the continuous system, the interaction between components
with wavenumbers k' and k" must vanish identically if either
k' is parallel to k or k1 | = k"|. However, the Arakawa J
_L
-17-
-------
algorithm leads instead to a non- vanishing, spurious coupling
in this case (Lilly, op.cit.)- (Other finite-difference schemes
lead to different spurious coupling relationships.) Moreover,
the misrepresentation of spectral interactions occurs every-
where throughout the spectral range, and cannot be surmounted by
ignoring a particular subrange, as discussed below.
2.1.5 Aliasing
Another fundamental difficulty with finite-difference tech-
niques is "aliasing" associated with a finite grid interval.
The phenomenon of aliasing is intrinsic to the representation
of non-linear interactions defined over a discrete physical-
space grid. It can be illustrated in the following simple manner.
Consider the spatial interval of length L, over which is defined
an N-ippirit discrete grid • x. , x.=jAx, j = l? ... N (L=NAX)
Let us now represent the discrete analogues of the fundamental
2fTX
sinusoid f (x) = sin — — , and the first two harmonics
LI
4 71 }£ 6 TT 5C
g(x) = sin — — and h(x) = sin — - — . The discrete version of
LI LI
the fundamental sinusoid, defined over this grid, is just
f,(x_;) = sin — ^- and in like manner g , (x . ) = sin
,_; , .
d j N a ] N
h , (x . ) = sin 7r-' . Consider now the discrete-grid representation
d j N
of the two product fields
P(x) = f (x)g(x) = i(cos lj£ - cos
Q(x) = g(x)h(x) = |(cos - cos
We might then suppose that
P (x ) = -Tcos ^Li - cos
^ ; ^
N N
-18-
-------
and that
To make matters concrete, choose N=8. If we were to tabulate
P (x) and Q, (x) , we would find them to be numerically identical
over the discrete-grid j = l, ... 8, whereas the continuous-space
functions certainly differ with regard to the terms cos -j— ,
cos 107TX . what has occurred? Trivially, of course, we have im-
L
posed an identity, since, for N=8, cos -^- = cos ~~- • But
the physical consequence of this little example is that, over
a discrete grid of length eight points, the sixth harmonic of
the fundamental spatial waveform is indistinguishable from the
second spatial harmonic, and that, in fact, the information
contained in the sixth harmonic, if it occurs as a result of
non-linear interaction among, say, two lower modes, would be
misrepresented by, or aliased down upon, the second harmonic.
In like manner, the seventh harmonic would appear aliased as
the first harmonic, and so on. We have illustrated the notion
of a "folding" frequency (the Nyquist frequency N/2) associated
with a discrete spatial grid of finite length N points.
Let us examine the aliasing problem in more formal mathe-
matical manner. The values of the arbitrary functions f(x), g (x)
at N equally spaced grid points xn = -— - (n=0 , 1, ... N-l) may
be expanded in the discrete Fourier series
f(x ) = f - I U(k)elkXn
(2.19) n n |-k.|
-------
where k is an integer satisfying -K^k
-------
approximation to the Navier-Stokes equations may lead to difficulty,
because of aliasing errors in the discrete-grid approximation of
the (scalar) product v-V v- Phillips (1959) gives a simple example
where aliasing terms in the representation of a product induce an
instability not present without aliasing. This aliasing property
of finite-difference representations has several other consequences
in the numerical simulation of non-linear fluid flows. We must
recognize,- first, that the physical process of non-linear turbulent
spectral cascade of energy (or vorticity) to small scales is falsi-
fied, inasmuch as the folding frequency acts, inthe absense of
dissipation, as a reflective boundary in wavespace. That is,
since no modes above the Nyquist frequency exist to receive energy,
the energy transferred by non-linear interaction must be redistri-
buted among the available modes which lie below the Nyquist cutoff,
in order that energy remain a conservative quantity. [This is
true, as well, of vorticity, or enstrophy]. Accordingly, it might
suggest itself to ignore the small-scale structure, in the expec-
tation that the largest scales are at least adequately represented.
But, in fact, all the scales of the system are contaminated by
downward-aliased spectral information.
One may ask, therefore, why conventional finite-difference
methods are reasonably successful for certain classes of geophysica]
fluid flow simulations, most notably, for example, in large-scale
operational numerical weather prediction. The answer is rather
subtle, and concerns the spectral structure of the flow fields
being simulated. If the spectrum falls off sufficiently rapidly
with increasing wavenumber, then the aliasing-induced contamination
of the lowest modes, while present, may be so small as to be ignorabl
-21-
-------
It is precisely this property of large, geostrophic-scale atmos-
pheric flow systems (which seem to follow minus-third or minus-
fourth spectral power laws) that permits physically credible
numerical weather forecasts on these scales. In general, how-
ever, the spectral power law of the simulated physical system does
not fall off with sufficient rapidity, so as to permit the ignoring
of aliasing contaminations. For such problems, therefore, it is
especially important that the numerical technique employed can
eliminate aliasing errors or strongly limit their extent, if
such errors are deemed important. Appendix 3 presents the method
developed by Orszag for efficient alias-free multiplication of
fields defined over discrete spatial grids.
2.1.6 Treatment of Singularities and Local Discontinuities
Among the most troublesome problems encountered in numerical
flow simulation is the treatment of flow structures with singular-
ities and local discontinuities embedded in such flows, as for
example, point sources of pollution. Eulerian finite-difference
methods permit the numerical evolution of strongly discontinuous
flows, but only at the cost of smearing out, or smoothing over,
these discontinuities, often in an entirely unrealistic manner.
These effects have been demonstrated by Orszag (1971b),who compared
second and fourth order finite-difference methods with truncated
spectral methods, for the problem of the two-dimensional advection
of a passive scalar field, whose initial distribution is circular
(hence referred to as a scalar "cone")/ by a uniform rotation
velocity field. (An extension of this study to viscous flows
is included in Chapter 5 of the present report.) One of the most
-22-
-------
remarkable features of the truncated spectral method that is
demonstrated by this problem is the result that the cone is
better localized in space, using these truncated spectral repre-
sentations, than by finite-difference schemes formulated directly
in physical space. That is, the cone location is preserved more
accurately by expanding the initial excitation in discrete Fourier
series, permitting these Fourier coefficients to evolve in time
in accordance with the spectral form of the governing dynamical
equations, and finally, reconstructing the discrete real-space
field from the Fourier coefficients, than by a direct evolution
of the real-space field using a finite-difference method. In
other words, flow structure with intermittencies in physical
space seem to be more accurately followed in terms of their Fourier
components than directly in physical space. To explain this rather
counter-intuitive result, we discuss in Appendix 1 the data content
of discrete real-space representations of continuous spatial func-
tions. The principal consequence, in the present context, is that
spectral representations (as by Fourier expansions), allows much
better interpolation between discrete real-space grid points than
is permitted by the structural data preserved by the local grid-
point values alone. That is, the global information preserved in
the spectral decomposition of a spatial function permits more pre-
cise reconstruction of local details than does the grid-point
decomposition of corresponding resolution. This may be better
understood by reference to the previous discussion of the phase
and aliasing errors which underlie finite-difference methods; in
particular the point made that finite-difference schemes have
mostly lagging phase errors, so that some of the data needed to
reconstruct a spatially local disturbance at a given instant in
-23-
-------
time is erroneously phase-lagged and so is as yet unavailable
at the appropriate local spatial region.
-24-
-------
3. 0 Spectral Methods
The utility of spectral transformations in the analytic
solution of many partial differential equations of classical
physics has been long recognized. It is only relatively recently,
however, that advances in available computational power and algor-
ithmic breakthroughs have made possible the wide use of numerical
spectral methods for the solution of PDE's, in particular, for
the solution of a variety of geophysical and other fluid flow
problems.
The essential principles which govern the use of spectral
methods in numerical simulations may be stated quite simply;
instead of solving multidimensional finite-difference equations,
one applies semianalytic methods in which the dependence of the
flow on one or more independent variables (generally, but not
exclusively, the spatial dimensions) is expanded into a series
of smooth, and usually orthogonal, functions. The potential
choice of orthogonal expansion functions is quite broad; e.g.,
Fourier, Chebyshev, Lengendre, and Hermite series, among others,
are possible candidates in various different flow problems. The
specific choice made depends upon the following kind of consider-
ations :
(a) the nature of the boundary conditions to be imposed,
and the suitability of the expansion functions chosen
in regard to the imposition of these boundary conditions;
(b) the ease with which the governing differential equation
can be formulated and numerically evolved with time, in
terms of the expansion functions;
-25-
-------
(c). the relative computational effort involved in the
numerical evaluation of the expansion coefficients;
(d) the accuracy achieved in the expansion, in particular,
the rate of convergence of the expansion coefficients;
(e) the relative generalizability which a given class of
orthogonal polynomials offers in a variety of flow
problems.
A few simple examples will suffice to illustrate these
points. For example, in flow in cylindrical geometries the
velocity field must be periodic with respect to the zonal coor-
dinate so that the dependence on cj> is expandible in a
Fourier series
oo
V (r,
-------
(3.1) + = 0 (
d t d X
with the following mixed initial-boundary value conditions:
(3.2a) U(x,0) = 0
(3.2b) U(-l,t) = sin irni t.
The exact solution is just
U(x,t) = 0; (tx+l).
(The index m designates the number of full wavelengths within
the interval |x|<
-------
the boundary conditions. Chebyshev polynomials (Orszag 1971d)
are particularly appropriate. The nth-degree Chebyshev polynomial
Tn(x) is defined by T (cos8)=cos n9 . If u(x) has p piecewise
continuous derivatives on the interval [-1, 1] , then u(x) may be
expanded as (Fox and Parker 1968, Orszag 1971e)
00
(3.3) u(x) = I AnT (x)
n=0 n n
(3.4)
In particular, if u(x) is infinitely dif ferentiable , then the
Chebyshev coefficients go to zero faster than any finite power
of 1/n ; we say that the approximation obtained by truncating
(3.4) at N terms is of infinite order, since the error goes
to zero faster than any finite power of 1/N. Also, the latter
Chebyshev series may be differentiated termwise any number of
times. Notice that the convergence rate of the expansions is
governed by the smoothness of the function u(x), irrespective
of boundary conditions.
At present, there are three principal classes of discrete-
spectral methods for the representation of continuous systems
such as (3.1). These are the Galerkin approach (Collatz, 1960),
the tau method (Lanczos, 1956), and the collocation method
(Collatz, 1960, Lanczos 1956). The three methods can be use-
fully -compared in regard to equations (3.1-3.2). The Galerkin
methods employ a new dependent variable v (x, t) =u (x, t) -sin mTrt,
so that v(x,t) satisfies the homogeneous boundary condition
v(-l,t)=0 and the inhomogeneous wave equations, and equations
-28-
-------
for the time-dependent coefficients are obtained by requiring
that an appropriate inner product of the dynamical equation
(3.1) with each of the retained Tn(x) be zero. When rewritten
in terms of u(x,t), the semidiscrete Chebyshev-Galerkin equa-
tions become (Orszag 1971d)
N
(3.5) u(x,t) = I a (t)Tn(x)
n=0
da N
(3.6) c -J-P- = -2 V pa + (-l)nb(t) (0l) . Here b(t) is a "boundary" term
that is determined by compliance with (3.7).
The tau method is based on using the finite expansion
(3.5) to solve exactly the approximate differential equation
u +u = T (t)T (x). The resulting equations for a (t) are
L x j i n ti
precisely (3.6) and (3.7) but with the boundary term b(t) set
to zero for n=0, ..., N-l and b(t) = TM(t) for n=N.
Finally, collocation uses the spectral representation (3.5)
to interpolate values of u(x,t) at (N+l) "selected points"
x (p=0, ..., N) . With Chebyshev polynomials, it is most natural
to choose x =cosTp/N, so that (3.5) is just a discrete Fourier
transform. The polynomial interpolation (3.5) is used merely to
evaluate 3u/9x at the points x by termwise differentiation.
In terms of a , collocation gives (Orszag 1972) equations (3.5-3.7)
with b(t) replaced by b(t)/C ; the only change is in (3.6)
-29-
-------
for n=N where (-l)Nb(t) is replaced by 1/2 (-l)Nb(t).
Collocation (called pseudospectral approximation in a later
section) is particularly convenient when there are no-neonstant
coefficients or nonlinear terms. In these latter cases, pseudo-
spectral approximation is quite different from, and much simpler
than, spectral (Galerkin or tau) approximation because of the
appearance of aliasing terms in pseudospectral approximation.
Notice that solution of 3.6, 3.7 requires only 0(N)
arithmetic operations per time step if care is taken to accumu-
late required sums, and that the transformation (3.5) to the
"physical space" u(x,t) can be accomplished in 0(N log N)
operations when N is a power of 2, because (3.5) is just a
Fourier cosine series for which fast Fourier transforms are
applicable.
-30-
-------
3 . 1 Boundary Conditions and Geometries;
3.1.1 Periodic Boundary Conditions
Physical intuition leads us to expect that for planar
geometries and/or periodic boundary conditions, the use of
discrete Fourier series is appropriate. In this case, our
intuition is indeed correct; finite Fourier series are the
most natural and efficient expansion method.
3.1.2 Free-slip Boundary Conditions
Free-slip--free-slip boundary conditions such as
(a) U =V =0 z=0, IT
z z
(b) W=0 Z=O,TT
are mathematically equivalent to symmetry boundary conditions,
inasmuch as conditions (a) are most naturally implemented in
cosine series and conditions (b) in sine series representations
of (U,V), and W respectively. Hence Fourier-series expansions
are the appropriate choice for free-slip--free-slip boundaries.
For mixed f ree-slip--np_-slip conditions, the use of Fourier
series would, in general, be ill-advised.
3.1.3 No-Slip Boundary Conditions
By contrast with free-slip conditions, no-slip boundary
conditions impose stringent constraints upon the dynamical fields,
which, in general do not lend themselves readily to computationally
efficient or accurate implementation in Fourier-series expansion
-31-
-------
due to Gibbs phenomena at the boundaries. Rather, as has been
demonstrated by Orszag (1972) Chebyshev expansion is the tech-
nique of choice, principally because the rate of convergence of
such .an expansion can be specified a priori, without reference
to the boundary conditions, and because Chebyshev methods are
efficiently implemented by skillful use of Fast Fourier (cosine)
transforms. Chebyshev series are not restricted to planar geo-
metries; they may profitably be applied, for example, to handle
the radial expansion in flows in cylindrical or spherical geo-
metries, or in any geometry that can be smoothly mapped, by con-
venient transformation, into either rectangular or circular
domains. Metcalfe (1973) has employed Chebyshev expansions in
the analysis of a variety of viscous incompressible flows in
various cylindrical and spherical geometries.
3.1.4 Inflow and Outflow Boundary Conditions
Spectral methods are capable of high-order accuracy at
inflow and outflow boundaries, for example, as in a limited
forecast area model (LFM), without the need for special methods
to impose the boundary conditions. By contrast, high-order
finite-difference methods deteriorate, or may break down com-
pletely, near boundaries because grid points outside the physical
domain must be invoked. The necessary modifications to maintain
accuracy near the boundary can get quite complicated (Kreiss and
Oliger 1973).
-32-
-------
3.2 Treatment of Boundary^Layers
Among the most valuable properties of Galerkin methods
using Chebyshev expansions, is the much greater fidelity with
which boundary layer structure can be resolved using Chebyshev
expansions, compared to a conventional finite-difference repre-
sentation with comparable number of grid intervals. This follows
naturally from the definition of the Chebyshev polynomials, namely:
(3.8) T (0) = cos n 8 , 0<8 0
i e = .
n
1, nN
In particular, for n=N-2, we have
(3.15)
_ = 4N(N-l)a
N
-35-
-------
Substituting (3.13) and (3.15) into (3.14) we may eliminate
the terms b , and obtain a linear tridiagonal matrix system
of the form
where a is the vector
L a
'NJ
rao
a is the vector
This system may be inverted by any of a number of standard
matrix inversion algorithms, the standard LUD algorithm being
especially useful.
Suppose Up(z) is of the form 1-
cos 26
(z=cos
Then the only non-zero modes in the Chebyshev representation
of this function are a =1, and a =-1/2.
u £
Let us now compare the Chebyshev solutions to (3.12)
(using (3.13) - (3.15)) as we vary parameters y and N,
the number of retained modes. As before we designate the Cheby-
shev modes of IF by a (n=0„ ..., N), and the modes of U
will be designated by a :
p+1
Case 1: Y =10 '
n=0
2
4
6
8
10
12
j
N=8
.9385
-.6093
-.0827
-.0451
-.2013
a
N=16
.927039
-.63308
-.10976
-.07966
-.05092
-.02863
-.01398
-36-
N=32
.927056
-. 63306
-.10976
-.07971
-.05102
-.02885
-.01443
n
1.0000
-.5000
0
-------
14
16
18
20
22
24
26
28
30
32
N=16
00538
00561
N=32
-.00641
-. 00253
-8.9 xlO
-2. 8 xlO
-8. 0 xlO
-2.0 xlO
-4.7 xlO
-9. 8 xlO
-1.8 xlO
-4.6 xlO
-4
-4
-5
-5
-6
-7
-7
-8
Note that for N=8 retained modes, no convergence (a +0 for n->-™)
is apparent. In the case N=16, there is just marginal conver-
gence, but for N-32 modes the convergence is extremely strong.
This indicates that there is just marginally enough resolution,
using 16 modes, to resolve the boundary layer whose characteristic
thickness is 6 ~— ~ .03, but that merely by doubling the number
/v
of modes, the resolution is about as perfect as we could ever
hope to achieve (another doubling would give no better resolution),
Case 2: Y =10
-2
n
n=0
2
4
6
8
N=8
. 857279
-.698663
-. 106009
-. 036378
-. 016238
N=16
. 857280
--698875
-.106911
--0391647
-.0101173
N=32
same as N=16
-37-
-------
N=8 N=16 N=32
10 -.00190785 -.00190785
12 -2.71303 xlO~4 -2.71333 xlO~4
14 -2.966 xlO~5 -2.991 xlO~5
16 -3.089 xlO~6 -2.619 xlO~6
18 -1.86 xlO~7
20 -1.09 xlO~8
22 -5 xlO~10
24 -2 xlO'11
26 -5 xlO~13
28 -2 xlO~14
30 -6 xlO~16
32 -2 xlO~17
In this case, for the characteristic boundary layer thickness
6 ~ .1, N=8 modes is just marginally adequate, while 16 modes
already displays rapid convergence, and an additional doubling
gives essentially nothing more.
Orszag has found an approximate'rule which relates the
boundary layer thickness implicit in the computation, to the
number of modes just necessary to adequately resolve such
a layer:
.. ^ ^-1/2
NCRIT 36
Note that for Case 1, NCRIT ~ 16.8>,. .and for Case2f NCRIT ~ 9.5,
which is entirely consistent with the results shown.
Stated in more familiar terms, the number of Chebyshev
modes required for an accurate spectral simulation of a laminar
1/4
flow with a viscous boundary layer increases only as R ' (where
-38-
-------
R is the Reynolds number of the flow), whereas the resolution
requirement of a finite difference model varies as R
For the high Reynolds numbers characteristic of realistic mete-
orological flow problems, the computer storage advantage of
Chebyshev methods is enormous. Such methods thus offer the
potential of extremely accurate numerical simulation of boundary
layer structure.
If, in conjunction with these Chebyshev methods, appropriate
coordinate-stretching techniques (Israeli 1971, Ph.D. Thesis) are
employed to assist in the resolution of a boundary layer of thick-
ness S , then the spectral errors are decreased faster than any
finite power of 6 as S+0, whereas the errors in finite differ-
ence methods, due to such coordinate transformations, decrease by
finite power of 6 .
3.3 Pseudospectral Methods
Our emphasis thus far has been focused upon the purely spec-
tral (Galerkin) methods. We have seen that, within limitations
imposed by geometry and boundary conditions, these methods cire verv
powerful. We have suggested also that Galerkin methods are rela-
tively easy to implement for certain important classes of non-linear
interactions, as for example the rotational term V x U) j_n the Wavier
Stokes equations. The accuracy with which phase information is pre-
served, and with which aliasing can be controlled, has been discussed
and the elimination of aliasing is described in detail in Appendix 3,
(1) in a finite-difference representation of vj , say
v zz -I/?1
—j~ (U. , - 2U . + u . _, ) i.t is evident that Az must vary wilh N
z 1/2
or eq;.ivalently with R
-39-
-------
However, there are many cases of interest for which aliasing
errors can be expected to be small, and/or tolerable, and we would
not wish to pay the computational price required to eliminate alias-
ing in such cases. Secondly, it is sometimes necessary to evaluate
highly non-linear terms within a spectral framework, and such terms
would be extremely difficult to formulate directly in wavespace.
Thirdly, if we wish to represent physical processes that are local
in real space (such as, for example, the onset of condensation in
a numerical model of cloud dynamics) such a local-space condition
would be extremely awkward to formulate in a Galerkin framework
because a single point in real space is represented in wavespaco
as a superposition of many modes, and a product term in real space
corresponds to a convolution sum in wavespace. Accordingly, Orszag
has introduced a modified Galerkin-like approach, which he has termed
"pseudospectral." This approach is computationally very economical,
retains much of the power of the purely spectral Galerkin method,
yet permits us to handle more general non-linear terms, complicated
boundary conditions, and processes which are local in
real space. The idea behind the pseudospectral method is to repre-
sent certain terms in physical space when it would be computationally
prohibitive to do so spectrally, but to switch back to wavespace
when it is very important, and computationally easy, to preserve
accurate phase and amplitude information.
Consider, for example, the temperature equation written in
an incompressible Boussinesq system
(3.16) ^~ f£ = ~ ^~ V-(V T) - N2(z)W + -3- KV2T
Too at Too ~ oo
-40-
-------
where T is the variable temperature field, TQO is a standard
reference temperature, V is the three-dimensional velocity vector,
W is the vertical velocity component, N(z) is the Brunt-VaLsala
frequency, which depends upon altitude z, g is the acceleration
of gravity, K is the (constant) thermal diffusion coefficient.
An appropriate pseudospectral procedure for evaluation of the
right hand side might be as follows:
(a ) Represent each component of V, and the temperature
field T, in real apace, and form the vector quantity
(VT) at each physical grid point, that is
VT = U(x,y,z)T(x,y,z)i + V(x,y,z)T(x,y,z)j
/-\
+ W (x , y , z ) T (x , y , z ) k ;
(a ) transform the vector quantity (VT) to a vector term
in wavespace, denoted by (VT).
(a.,) Evaluate the divergence of the quantity directly in
wavespace by multiplying by the wavespace divergence
operator ik- (where k_ is the wave vector) ; this has
the virtue of preserving the exact phase relationships
in the divergence expression, without the first-differ-
encing errors introduced by a finite-difference approxi-
mation to the divergence term in real space.
(b) Assuming an explicit evaluation of the thermal diffusion
term is adequate, transform T to wavespace, denoting
its transform by T, and evaluate the Laplacian term
directly in wavespace by multiplying T by the scalar
2
operator - k -k = -k . This procedure avoids the second-
differencing errors associated with the finite-difference
form of the Laplacian operator. [The diffusion term can
-------
be treated implicitly in a spectral formulation very
simply, if, for example, (3.16) is written in the
transformed form
(VT) - (N2(z)W),
T
00 00
where, as, before, the tilda (~ ) designates the wave-
space transform of the real -space field.]
(c) Represent the vertical velocity W in real space, and
at each point of the physical grid form the product
2 2
N (z)W; then transform the expression N (z)W to
wavespace. By first evaluating the product in real
space, we avoid a nasty convolution sum, but accept
the inevitable aliasing penalty, which may not be a
severe limitation, especially if the term is small
or if N is a slowly varying function of z.
(d) Since each term on the right hand side of (3.16) is
now available in its spectral form, we may evaluate
j\ m
each of the spectral components TTT. Finally, then,
d t
r>"T
we inverse transform v~- to real space to obtain
o t
the local rate of change of temperature.
Thus, the pseudospectral method makes effective use of finite
spectral operations when it is convenient and/or crucial to do so,
but returns to finite real space collocation operations when the
equivalent spectral operation is computatively prohibitive, or
aliasing errors are tolerable. A more elegant statement of the
method, which has perhaps been adequately motivated by the example
-42--
-------
above, has been provided by Orszag (1971a). That is, the
philosophy of the pseudospectra.l method is to calculate in
either spectral or physical space according t.o whichever repre-
sentation is more natural. In particular, one chooses a suitable
orthogonal expansion as an interpolatory tool, in order to be
able to differentiate without phase error. Derivatives are evalu-
ated entirely by explicit multiplications in spectral space (rather
than by finite differences in real space). Operations like differ-
entiations in finite-difference schemes and convolutions in spec-
tral schemes are avoided.
Let us apply the pseudospectral idea for efficient evaluation
of the non- linear terms of the equations of motion. The principal
error of finite-difference schemes is first-differencing or phase
error (cf. section 2.1). These errors are avoided so long as spec-
tral methods are used to compute derivatives. For example, 9v/9x
at x=x . , where x. = 2Trj/N (j = l, ..., N) , may be computed without
phase error by first Fourier analyzing v-=v(x.) into its discrete
Fourier components v(k) (|k| < N/2) and then setting
A term like v9v/3x in an equation of motion can then be evaluated
at the grid point x. as the local physical-space product of v(x-)
with Ov/9x) . If derivatives are computed in this manner,
.A. ^C «
the resulting scheme is not subject to phase error, but it is "fully
aliased" in the sense that spectral terms modulo +2k appear in
the spectral analysis of the scheme (see section 2.1.5).
-43-
-------
If a pseudospectral approximation is made to the Navier-
Stokes equations written in the form (2.2) by evaluating deriva-
tives as in (3.18), the resulting fully-aliased equations may be sub-
ject to unconditional instability (Orszag, 1971b). However, if
the Navier-Stokes equations are first expressed in rotation form
(3-19) j£ V(x,t) = V(x,t)x
-------
4.0 Applications of Numerical Spectral Methods
In this section we discuss briefly some of the principal
areas of fluid research, relevant to meteorological problems,
in which numerical spectral methods have been found especially
fruitful, namely:
(a) the simulation of atmospheric dynamics;
(b) the simulation of turbulence and tests of turbulence
theories;
(c) advective-diffusive modeling.
Recent years have witnessed an avalanche of papers concerning
spectral simulations of various problems in atmospheric flows and
numerical weather prediction. For example, Eliasen, et al. (1970),
have considered spectral harmonic methods for integration of the
primitive equations on the sphere; Merilees (1973, 1974) has de-
scribed pseudospectral methods for such problems; Bourke (1972),
Gordon and Stern (1974) have described experiments with multilevel
baroclinic global general circulation models constructed in spec-
tral manner. An especially important spectral model is the strato-
spheric dynamics model of the MIT group (Cunnold, et al. 1975) which
incorporates elements of both the purely Galerkin and the pseudo-
spectral approaches for the simulation of the dispersion of chemical
pollutants in the upper atmosphere. Many of the recenb advances
in the field of atmospheric spectral modeling are described in
the Copenhagen symposium volume (GARP, 1974), to which the reader's
attention is directed.
The employment of spectral methods for direct numerical simula
Lion of turbulent flows is of very recent origin. Spectral simula-
-45-
-------
tions of this sort, which require explicit retention of 1-3x10
modes, have become practical only with the introduction of effi-
cient transform techniques (Orszag 1969, Patterson and Orszag
1971). Simulations have been performed for homogeneous two-
dimensional turbulence (Herring,et al. 1974) and homogeneous
isotropic three-dimensional turbulence (Orszag and Patterson
1972), at moderate Reynolds numbers. Extensions of these studies
to inhomogeneous flows, turbulent wakes and jets, etc., are being
actively pursued (Herring 1974, Orszag and Pao 1974).
It would not be appropriate here to attempt to review the
salient conclusions regarding the structure of turbulence, and
the validation of turbulence theories, which have emerged from
such studies. One point, however, is certainly worth emphasizing
here. That is, throughout the course of such spectral simulations,
it has become empirically evident that the large-scale flow fea-
tures exhibit a remarkable degree of Reynolds-number independence.
As a consequence, it may not be necessary to simulate turbulence
at the huge Reynolds numbers appropriate to the atmosphere, if
our principal concern is with synoptic-scale features; that is,
much valuable information can be obtained from simulations whose
Reynolds numbers are such that all dynamically relevant scales
of motion can be included explicitly in the computation, without
recourse to closure arguments or subgrid-scalemodeling. It appears
that, once a certain critical Reynolds number is attained (and such
Re . are remarkably low!) the large-scale features are sensibly
independent of the detailed mechanisms by which energy (or ens trophy)
are cascaded to the smallest scales, nrovided that the initial flow
conditions are Reynolds-number independent.
-46-
-------
Of more immediate importance, with regard to the improvement
of air-quality simulation modeling efforts, are the insights which
are emerging concerning the utility of spectral methods in the
numerical simulation of advective-diffusive problems. In a pio-
neering effort, Orszag (1971b) considered the Crowley color pro-
blem (Crowley, 1968), namely the two-dimensional convection of
a passive scalar field by a uniform rotation velocity- This
problem has been considered as well by Molenkamp (1968) , from
the finite-difference point of view. Orszag compared solutions
to this problem using second-order and fourth-order Arakawa
finite-differencing schemes, and the Galerkin Fourier approach
discussed before. He demonstrated that results obtained using
the truncated Fourier-expansion scheme on a 32x32 space grid
(cutoff wavenumber K=16) are at least as good as those obtained
by the fourth-order scheme on a 64x64 grid, and significantly
better than those obtained by the fourth-order scheme on a 32x32
grid, or by the second-order scheme on 32x32 and 96x96 grids.
It appears that, for those problems for which the "color" problem
and the Neuringer problem discussed in Chapter 5 are typical exam-
ples, second-order schemes require at least four times as many
grid intervals in each spatial direction, and fourth-order schemes
require at least two times as many grid intervals, as the truncated
Fourier expansion scheme.
We have recently extended this study to :i ncorporate diffusive
terms. These new results are presented in Chapter 5 of this report
As will be seen, the conclusions which emerge support strongly the
argument for the use of spectral techniques where feasible, espe-
cially in regard to high Reynolds number flow simulations.
-47-
-------
5 - 0 A compg£i£gn__of Spectral and Finite Difference Methods
for Two Examples
To date, little effort has been directed toxvard exploiting
spectral techniques for numerical air-pollution simulations. This
section illustrates some of the computational advantages of the
spectral technique for the simulation of advective-dif f usive prob-
lems, by reference to (a) a diffusive generalization of the Crowley
"color " problem (Crowley 1968, Molenkamp 1968, Orszag 1971a),
that is, the advection-dif fusion of an instantaneous point-source
scalar in a two-dimensional field of solid rotation; and (b) a prob-
lem proposed by Neuringer (Neuringer, 1968) , concerning the advection
diffusion of a point source in a sheared wind field. We compare
numerical solutions using the spectral (truncated Fourier) method
with those using the conventional second-order finite-difference
method. It is shown that, for simulations of comparable accuracy,
the spectral method permits very substantial reductions in required
core storage, and less computational time as well.
We also present a rebuttal to some criticisms offered in a
private communication by Chan and Stuhmiller (denoted by C-S) to
a preliminary version of these results.
5 . 1 Comparison of Methods for the Diffusive Color Problem
We follow the evolution of the two-dimensional system
6
-------
rotation; \; is a diffusion coefficient; and the
-------
where k^'k2 represent the two horizontal wavenumbers ^
and the tilda (~) designates the two-dimensional discrete Fourier
transform of the indicated field.
In particular, the term
k k
K1'K2
designates the (k ,k ) component of the alias-free product of the
velocity component U with the-real-space field C at time step n.
The achievement of alias-free non-linear products is performed as
described in Appendix 3. The time-differencing used for the advec-
tive terms is the Adams-Bashforth scheme, while the viscous term is
represented with Crank-Nicholson time differencing. The expression
(5.3) makes use of the fact that the constant advecting velocity
field is non-divergent.
To emphasize the power of the spectral approach, the finite-
difference computations are performed with both 32x32 and 96x96 points
grid resolution, whereas the spectral calculations are performed only
on a 32x32 point grid. In each case we used a true point singularity
as the initial source term, of non-dimensional magnitude 1.0 units.
Thus the initial conditions are identical in each case, which satis-
fies an objection raised by (C-S) to the preliminary results. Pre-
vious experience with spectral calculations of the color equation
show that the spectral approach is remarkably accurate for a distrib-
uted scalar field (Orszag, 1971a) but should-be appreciably less
accurate for a true point source. In every case the numerical
time step was chosen sufficiently small that time-differencing errors
-50-
-------
were negligible.
For each system, the physical-space grid is square, and the
center of rotation coincides with the center of the grid. The
initial point-source distribution is located midway between the
center of rotation and the left boundary of the grid. The fre-
quency of rotation is such that one full rotation requires 100
time steps.
We present results for the following cases:
(a) Finite-Difference Scheme with 32x32 point resolution,
viscosity v = .05 w
(b) Spectral Scheme with 32x32 point resolution (i.e. wave-
number cutoff K=16), viscosity v = .05
(c) Finite-Difference Scheme with 32x32 resolution, v = .005
(d) Spectral Scherfte with 32x32 point resolution, v =.005.
Figure 5.1 displays the initial point-source distribution,
that has been enticofLy machine-constructed with ten equal contour
intervals, dashed lines representing negative-valued contour lines.
The contour intervals were generated automatically with contour
spacing
C - C
. _ max min
AL ~ 10
where C , C . are respectively, the maximum and minimum values
of the scalar field C encountered anywhere on the grid at the given
time step. Thus, no two figures show the same contour spacing, a
point which will become important shortly. However, the choice of
contour interval is now made in entirely objective manner, as was
not the case in the preliminary results seen by (C-S).
-51-
-------
Figure 5.1 Initial Point-Source Distribution (t=0)
'','."('
Contour from 0. to l.OOE+00
-52-
-------
Figure 5.2 Finite-Difference Scheme
(v = .05, t = 1/4 revolution)
32 x 32 point resolution
F.D. Scheme T = 25.0000
Contour from -1.42E-02 to 5.31E-02 Interval 6.73E-03
-53-
-------
Figure 5.3 Spectral Scheme
32 x 32 point resolution
v = .05
t = 1/4 revolution
Spectral Scheme T = 25.0000
Contour from -1.19E-03 to 6.21E-02 Interval 7.91E-03
-54-
-------
Figure 5.4 Finite-Difference Scheme
32 x 32 point resolution
v = .05
t = 1/2 revolution
F.D. Scheme T = 50.0000
Contour from -4.58E-03 to 2.84E-02 Interval 3.30E-03
-55-
-------
Figure 5.5 Spectral Scheme
32 x 32 point resolution
v = .05
t = 1/2 revolution
Spectral Scheme T = 50.0000
Contour from -5.92E-04 to 3.15E-02 Interval 4.01E-03
-56-
-------
Figure 5.6 Finite-Difference Scheme
32 x 32 point resolution
v = .005
t = 1/4 revolution
F.D. Scheme T = 25.0000
Contour from -1.09-01 to 1.13E-01 Interval 2.22E-02
-57-
-------
Figure 5.7 Spectral Scheme
32 x 32 point resolution
v = .005
t = 1/4 revolution
C>
o o
Spectral Scheme T = 25.0000
Contour from -6.38E-02 to 3.70E-01 Interval 4.34E-02
-58-
-------
Figure 5.8 Finite-Difference Model
32 x 32 point resolution
v = .005
t = 1/2 revolution
F.D. Scheme T = 50.0000
Contour from -5.02E-02 to 7.58E-02 Interval 1.26E-02
-59-
-------
Figure 5.9 Spectral Model
32 x 32 point resolution
v = .005
t = 1/2 revolution
o,
• o
o
Contour from -5.62E-2 to 2.55E-1 Interval 3.11E-1
-60-
-------
In figures 5.2, 5.3 we compare, respectively, cases (a) and
(b) after one-quarter resolution. Note that the contour interval
in the spectral case (5.3) is about 17% larger than that in the
finite-difference case. It is observed that case (b) more accu-
rately reproduces the circular symmetry to be expected, while case
(a) displays a small phase lag, relative to the expected center
of symmetry, axial distortion, and a spurious numerical noise
manifested as characteristic negative "trailing wake" (the explan-
ation for which we have considered in section 2.1). Proceeding
with the temporal evolution of both systems out to t=l/2 revolu-
tion, figures 5.4, 5.5 show that the spectral run continues to
yield a distribution which is quite circularly symmetric, with no
wake structure, and preserves the expected center of mass of the
distribution much more closely than does the finite-difference
run, wherein the axial distortion and phase lag are becoming quite
pronounced. The differences between these two cases (for v = .05)
are not yet striking, however. The situation changes greatly as
v->0.
Let us now compare the two methods for a ten-fold reduction in
"iscosity, v = .005, i.e., cases (c) and (d) as shown in figures
5.6, 5.7, for t=l/4 resolution. With a reduction in v, we
v;ould expect that the effect of the strong discontinuity in the
source distribution would become more apparent, since the smoothing
effect of the diffusion term is correspondingly reduced. Now the
spectral scheme shows to considerably better advantage with regard
to preservation of phase, symmetry of distribution, and suppression
of superfluous trailing wake structure. In fact, it would already
be difficult to argue that the finite-difference model corresponds
-61-
-------
to physical reality in any useful sense. Let us look as well at
i-he results for t=l/2 resolution (figures 5.8, 5.9). The contin-
ued high fidelity of the spectral model to the expected continuum
solution is truly impressive. Note that the phase information is
still preserved remarkably well, and the spurious numerical noise
is of small amplitude (< 10% of the peak).
The (32x32) spectral runs required about eight times longer
cpu execution time than the (32x32) finite-difference runs. How-
ever, it should be emphasized that these runs did not use optimally-
coded Fast Fourier Transform routines (which would save about 50%
execution time). Moreover, we employed for these runs a fully
alias-free spectral algorithm, and Orszag has shown that the partly
aliased pseudospectral method offers, in this case, a four-fold re-
duction in required computation time, without appreciable degrada-
tion of results. Together, therefore, these modifications would
effect at least an eight-fold savings in computation time. These
results suggest therefore that spectral methods are cost-effective for
simulation of the diffusive color problem and offer significant core
storage savings compared to finite-difference methods, for equiva-
l3nt accuracy.
Finally, it is necessary to correct a misconception entertained
by Chan and Stuhmiller, as evidenced by a comment of theirs concerning
the use of spectral techniques for the color problem:
"The example used in the comparison is a linear problem, i.e.,
the velocity field is given explicitly ..."In spectral method, there
are two general approaches to nonlinear problems in which the non-
linearity is of the product type. The first approach is to do the
entire calculation in the Fourier space; the nonlinear products
-62-
-------
lead to evaluation of 'convolution sum' which is in general very
time consuming. The second approach is to transform from Fourier
space into physical space and perform the multiplication there for
the nonlinear terms. When transforming back into Fourier space,
"aliasing" has to be controlled to avoid the buildup of errors.
Bass was probably referring to the second approach as the "spectral
method" in his memo. The first approach is not susceptible to phase
errors since the amplitude and phase of each Fourier mode are dealt
with directly as dependent variables. But the evaluation of con-
volution sums for nonlinear problems makes this approach rather
expensive. On the other hand, the second approach is not free of
phase error because of the calculation performed in the physical
space."
In response to these comments we should note that, first, the
color problem is "linear" in the sense that the velocity fields do
not enter quadratically, but inasmuch as it is necessary to evaluate
the advective product V-VC where V=V(x,y) is a spatially variable
field, the computation is nonlinear in the usual sense; and we require
an accurate efficient technique for the alias-free representation of
multiplicative fields.
Second, it is not true that the "spectral method" to which we
referred is implicitly aliased, or is phase-error contaminated. The
spectral method (a truncated Galerkin Fourier method) which we used
in these calculations is
(a) entirely free of phase error
(b) entirely free of aliasing misrepresentations, since we
used the alias-renewal technique described in Appendix 3.
Apparently Chan and Stuhmiller have confused spectral with
pseudospectr; 1 methods. The alias-free evaluation of nonlinear
-63-
-------
products does not require the explicit evaluation of convolution
sums. Perhaps the most single advantage of the methods introduced
by Orszag is that we can achieve the same results as would be ob-
tained by computationally expensive convolution, with much greater
efficiency via the "shifted-grid" algorithms (Orszag 1971a, 1971b).
-64-
-------
5.2 Comparison of Methods for Neuringer's Problem
In this section we compare spectral and conventional second-
order finite-difference numerical solutions to the problem solved
analytically by Neuringer (1968), namely, the advection-diffusion
of an instantaneous point source in a two-dimensional linearly
sheared wind field.
Neuringer presents an analytic solution for the evolution
of the two-dimensional system
(5.4)
CTT = d
3y
(S(X_X', y-y', t)
where the 6-function represents an instantaneous source and
a, c, d are constants. He shows that with the transformations
t = at
2
x = (a/c)(x-aty-act /2)
y = (a/c) (y+ct)
and with the 5-function source at the origin, x' = y' =0,
the analytic solution is
(5.5)
where
F(x,y,t) =-
exp
12
-y
(x/2-t.y/4)
t(l+t2/12)
2 2 . ,
Y = c /ad .
Figs. 2-5 in Neuringer's paper attempt to show the variation of
the dirnensionless distribution function F. It was observed that
these figures are entirely incorrect, except for the special para-
-65-
-------
etric choices x = 0 or y = 0. We have therefore reevaluated
numerically Neuringer's analytic expression for F(x,y~,t), equa-
tion (5.5) above. Fig. 5.10 exhibits some correct representative
curves of F vs. x for the case y = 1.0, t = 0.5. (Note, for
example, that by comparison with fig. 5.10, Neuringer's curve for
y = -5 is wrong by more than six orders of magnitude, while the
curve for y = 0 is correct.)
Superimposed on fig. 5.10 are some points sampled from a spec-
tral simulation of equation 5.4 with 32x32 point resolution, and
the same choices of parameters y, t, and (x,y) dependence. On
the scale of the figure the spectral results are indistinguishable
from the analytic results. Table 5A presents some sampled values
for the ratio F (x,y )/F(0, 0) spaced at unit intervals of x and y,
obtained from direct numerical evaluation of equation 5.4. In this
and subsequent tables the notation 2 (-8) is to be understood as
— P
2 x 10 , and so forth. In Fig. 5.11 we display a two-dimensional
contour plot of this ratio, constructed from more detailed numerical
tables. The corresponding contour plots constructed from the spec-
tral and finite difference runs cannot be distinguished by eye from
those shown in the figure and so are not included.
To permit more exact comparison of the spectral and finite-
difference computations for Neuringer's problem, we have performed
a second series of three runs, with parameter y now chosen for
convenience as 4.0, and the computation carried out to time
t = 0.5. For each run, the time increment was chosen small enough
so that time-differencing errors were negligible. The three runs
compared were:
(a) spectral computation on a 32x32 point space grid
(hereafter called run (a)) ;
-66-
-------
Figure 5.10
Legend
—• Y » o
Y • -2
V « -3
-•- y . .5
».... Y • -» 2
Y ' +3
Y » + 5 curve is
off -scole to lower
right ,
Figure 5.10 Correct Neuringer FCN for j =1, t = 0.5
A = Points Obtained from Spectral Run
for Same Combinations (X,Y)
-67-
-------
Fiqure 5.11
--2
Figure 5.11 Neuringer Problem (y = 1.0, t = 0.5)
Analytic Value of F(x,y)/F(0,0)
-68-
-------
TABLE 5A
Analytic Evaluation of F(x,y)/F(0,0) for y = 1-0» t = 0.5
-4 -3 -2 -1 0 1 2 3 4
4 - 2(-8) 5.5 (-7) 6.40(-6) 2.7S(-5) 4.54(-5) 2.78(-5) 6.40(-6) 5.5 (-7)
3 4(-3) 2.53(-6) 6.10(-5) 5.53(-4) 1.88(-3) 2.40(-3) 1.15(-3) 2.03(-4) 1.40(-5)
I
^ 2 2.45(-6) 1.23(-4) 2.33(-3) 1.65(-2) 4.40(-2) 4.40(-2) 1.65(-2) 2.33(-3) 1.23(-4)
I 5. 29 (-5) 2.OS (-3) 3.08 (-2) 1.71(-1) 3.57(-l) 2.79(-l) 8. 21 (-2) 9.06 (-3) 3.75 (-4)
Y G 3.95(-4) 1.22(-2) 1.41(-1) 6.13(-1) 1.00( 0) 6.13(-1) 1.41(-1) 1.22(-2) 3.95(-4)
-1 1.02C-3) 2.46(-2) 2.23C-1) 7.59(-l) 9.70C-1-) 4.65(-l) 3.38(-2) 5.66(-3) 1.44 (-4)
-2 9.12(-4) 1.72(-2) 1.22(-1) 3.25(-l) 3.25(-l.) 1.22(-1) 1.72(-2) 9.12(-4) 1.81(-5)
-3 2.S2(-4) 4.17(-3) 2.32(-2) 4.33(-2) 3.73(-2) l.ll(-2) 1.22(-3) 5.0S(-5) 7.9 (-7)
-4 3.02(-5) 3.49(-4) 1.52(-3) 2.48(-3) 1.5K-3) 3.49(-4) 3.02(-5) 9.8 (-7) l{-8)
-------
(b) second-order finite-difference computation on a
48x48 point grid (called run (b)),with one and
one-half times greater horizontal resolution
than run (a);
(c) Second-order finite-difference computation on a
96x96 point grid (called run (c)),with three times
greater horizontal resolution than run (a).
The equations used in the finite-difference and spectral models
to simulate (5.4) are virtually identical to those described for the
color problem, e.g.,equations (5.2) and (5.3), the only differences
occurring in the representation of the advective terms. That is,
the velocity field U(x,y) is replaced by +(ay), and the velocity
field V(x,y) is replaced by the constant -c. The time-stepping
schemes used here are identical to those described for the color
problem.
Table 5B displays the values of the ratio F (x,y~) /F (0 ,0)
obtained by tabulating expression (5.5), for equally-spaced points
(x,y) spanning the region of maximum concentration amplitude.
Table 5C shows the absolute errors made by the spectral run, i.e.,
the differences between the analytic values of the ratio F(x,y)/
F(0,0) and those obtained from the spectral computation. In like
manner , Tables 5D, 5E display the absolute errors of the 48x48 point
and the 96x96 point finite difference runs, respectively.
Upon examination of these tables, we may draw several signi-
ficant conclusions. First, we observe that the errors made by the
run (b) are consistently about four times greater than, and uni-
formly in the same direction as, the errors made by run (c). Second,
we see that, in the region of maximum amplitude, both the (a) and
-70-
-------
TABLE 5B
Analytic Value of F(x,y)/F(0,0)
(y = 4.0, f = 0.5)
t
-
-
-
_
- 1
1/2
0
1/2
1
3/2
2
.019
,141
.368
,332
.104
,011
- 1/2
,104
.613
1.252
.885
.216
.018
X
0
,216
1.000
1,600
. 885
..169
.011
+ 1/2
.16-9
,613
,767
,332
.050
.003
+ 1
.050
.141
.138
,047
.005
<.001
-71-
-------
TABLE 5C
Absolute Errors in Spectral Run
- 1
- 1/2
+ 1/2
+ 1
t
-
-
-
-
1/2
0
1/2
1
3/2
2
<.001
<,001
+ .004
+ .007
+ ,003
<,001
-.001
<.001
+ .007
+ .013
+ .008
+ .003
-.002
-
+ .011
+ .01,9
+ .001
+ .005
-.003
-.001
+ .012
+ .016
+ .007
+ .002
-,004
-,004
+ .004
+ .009
+ .004
+ .001
-72-
-------
TABLE 5D
Absolute Errors in 48*48 Point Finite Difference Run
- 1
- 1/2
1/2
1/2
0
1/2
1
3/2
2
-.003
-.006
-.023
-.039
-.007
+ .004
-.013
-.017
-.056
-.085
-.006
+ .007
+ .020
-
-.028
-.063
<.001
+ .005
-.015
-.017
-.036
-.034
-.003
<.001
-.007
.-,006
-.009
-,006
<.001
<.001
-73-
-------
TABLE 5E
Absolute Errors in 96*96 Point Finite Difference Run
x
- 1
- 1/2
-I- 1/2
t 1
1/2
0
1/2
1
3/2
2
+ .001
-.002
-,006
-,010
-.002
t.OOl
-.003
-.004
-.014
-.022
-.001
+ .002
-.005
-
-.007
-.016
<,001
+ ,001
-.004
-.004
-.008
-.009
-.001
<.001
-.002
-.002
-.002
-.002
<.001
<.001
-74-
-------
(c) computations are of about equal accuracy and capture equally
well the spatially local peak structure (although the finite-differ-
ence approach is "local" and the spectral approach is "global'1 in
nature). Third, with some reference to the more extensive tabula-
tions from which these tables have been extracted, it is seen that,
in the wings of the distribution, the analytic result and the finite-
difference calculations decay to zero very rapidly. By contrast,
the spectral computation makes a small oscillatory error of absolute
magnitude less than 3*10~ that of the peak amplitude. This small-
amplitude residual oscillatory error is typical of spectral repre-
sentations for strongly peaked fields. However, as we have noted
previously, the error made by truncating a spectral representation
at N terms decreases more rapidly as N is increased than any
power of 1/N. Hence, for only modest additional increase in spec-
tral resolution, the magnitude of this characteristic oscillation
can be reduced to negligible proportions.
It is important to note that, for this two-dimensional compu-
tation, therefore, the core storage advantage of the spectral method,
in this case, is almost thirty-fold. This makes it practicable, from
the point of view of available core, to pursue a much larger class
of calculations than are presently possible with finite-difference
methods.
Again, as regards execution time, the spectral method is more
than competitive. Run (a) required about 36% more computation
time than did run (c). However, as before, the spectral calcula-
tions were made in fully alias-free manner, and an immediate four-
fold reduction in computation time is possible if we relax the
requirement for fully alias-free spectral interactions, by per-
forming the calculations is pseudospectral manner. Thus, for
-75-
-------
comparable accuracy, the spectral approach is to be preferred, on
economic grounds, for solution of advection-diffusion problems
typified by Neuringer's problem. We believe, therefore, that
spectral methods merit serious attention for prospective air-
quality simulation efforts.
-76-
-------
6.0 Spectral Modeling of the Turbulent Diffusion of a
Passive Scalar
This chapter describes a modeling effort, as yet unsuccessful,
using state-of-the-art spectral methods to study the turbulent dif-
fusion of a passive scalar source in a fully three-dimensional,
turbulent, sheared velocity field, without recourse to closure
arguments. The impetus for this effort stemmed from the expecta-
tion that such a model would provide an idealized set of turbulent
diffusion data, free of closure parameterization assumptions, against
which various air-quality turbulent diffusion models might be tested.
Moreover, it was anticipated that the model might permit the compila-
ation of model turbulent diffusion data, suitable for the evaluation
of eddy Austauch coefficients, for a variety of cases for which observa-
tional data are lacking, or for which existing data are unreliable.
These expectations seemed well justified, at the initiation
of this study effort, by the very real successes achieved by earlier
spectral numerical models of this type in reproducing the essential
dynamics of three-dimensional turbulent flows in a number of cases
(Orszag and Patterson 1972, Orszag and Pao 1974), and by the demon-
strated utility of the numerical results so obtained in testing the
validity of various formal turbulence theories. The model experi-
ments contemplated for the current effort presupposed, therefore,
the availability of computational resources (core storage and CPU
time) comparable to that employed in the previous successful efforts.
It was not until well along in the current contract study, and only
after the various computer codes had been written and checked out,
that the computational resources available to us were thenceforth
severely restricted. The critical limitation was that of core
storage. Th^ experiments had been planned on the premise that the
-77-
-------
codes could be run with (32)J modal resolution, as had previously
been the case at the Lawrence Berkeley CDC facility. Unfortunately,
new and unforeseen user demands upon that facility had the practical
effect of restricting our model experiments to only (16) modes.
Previous experience with homogeneous turbulence codes suggested
that (16) modal resolution wais but marginally adequate for the
evolution of physically realistic turbulent flows, and we had but
little hope for the prospective success of a (16) model for strongly
sheared flows. Nevertheless because of our problems with Lawrence
Berkeley we had no choice but to proceed with our experiments at
the reduced resolution. Initially we hold out some hope for partial
success because the largest scales of motion (which are the most
important for turbulent eddy transport) could still be adequately
resolved. However, after the fact, the difficulty of achieving
physically credible small-scale equilibrium statistics in the pres-
ense of a strongly sheared mean flow at the (16) resolution pre-
cluded our getting satisfactory results. (The attainment of equili-
brium is necessary to make possible a meaningful comparison of our
numerical model with the corresponding predictions of a turbulence
closure model based upon equilibrium turbulent conditions.)
The remainder of this chapter is devoted, first, to a thorough
description of the model we have developed. Then, a discussion is
given of the difficulties encountered, by reference to figures which
depict the evolution of the velocity and vorticity fields during
the course of one of our unsuccessful experiments.
-78-
-------
6.1 The Numerical Model
In this section we define the Bouss'inesq equations describing
the physical system to be simulated, and develop the spectral
forms of these equations which are actually implemented in the
computer model which has been constructed.
We begin from the familiar Boussinesq system of equations
for shallow layer convection in an effectively incompressible
fluid. The scaling arguments underlying the derivation of these
equations are presented in Spiegel and Veronis (1960) and Calder
(1967), and will not be restated here. The Boussinesq equations
may be written as
(6.1)
at
= -JJL
(6.2) « V-V = 0
(6.3)
oo at ioo
where
A is a unit vector in the vertical direction
V is the•three-dimensional velocity field
p is the deviation of pressure from its equilibrium
value P
T is the deviation of temperature from its equilibrium
value T , which is a linear function of height z,
Te = T00(1-Z/H)
T00' P00 are the surface equilibrium temperature and
density, respectively,
-79-
-------
2
N , the square of the Brunt-Vaisala frequency, =
ST
cr el
—^— ( • • + —), assumed constant.
P00 3Z cp
v, the kinematic viscosity, = y/Pnn/ where y is the
coefficient of viscosity
a, the thermal viscosity = cr/P , where a is the ther-
mal conductivity
Equation (6.1) is the momentum equation, incorporating the con-
tinuity condition (6.2). Equation (6.3) is the heat equation.
For convenience, we will define the variables P = p/p ,
02? -U
• = T"^ • , whence
00
(6.4) |=s + v-VV = - VP + 0A +vV2V
;\cy, o 9
(6.5) •£• + we + irw = KV .
The term 6A in (6.4) provides the convective coupling of tem-
2
perature to the vertical velocity field, and the term N W in
(6.3) represents the effect o.f basic-state stratification upon the
evolution of the thermal field.
Equation (6.4) may be written in the rotational form
(6.6) % - V x w = - V(P + ~ V-V) + 0X + VV2V
i dt — — 4 — —
where w = Vx v is the vorticity field. .
Equation (6.5) may be written, using the incompressibility condi-
tion, as
(6.7) + V-(V0) + N2W =
-80-
-------
Equations (6.6), (6.7) and (6.2), together with appropriate ini-
tial and boundary conditions, define the system to be simulated
spectrally.
To integrate this system with time, we use truncated Fourier
spectral methods. That is, we represent the velocity and tempera-
ture fields in truncated discrete Fourier series (Appendix 2) of
the form
(6.8) V(x,t) = 7 U(k,t)e1-'-
|k|
-------
where a,3,y may each range over the three indices 1, 2, 3
(i.e., the three orthogonal directions in wavespace). The symbol
6 is the Kronecker delta, and summation over repeated Greek
indices is implied. The pressure term has been eliminated by
use of the incompressibility condition (see also the discussion
in section 2.1.3}. Note, however, that the form (6.10) requires
that a convolution sum be evaluated for each pair of indices
(3,y). Such an approach is quite costly computationally,
especially for a fully three-dimensional system. Instead, we
have employed the pseudospectral approach to integration of the
momentum and heat equations. Returning to equation (6.6), we
see that the source of the convolution term is the nonlinear
vector product term VXUK To avoid the convolution requirement
we use the following procedure:
(a) Beginning from the Fourier transform of. the velocity
field, U(k,t), create the Fourier transform of the
vorticity field,n(k,t) by performing the wavespace
operation J2(k,t) = i k x U(k,t).
This is the direct analogue of to = VxV. Note that because
the vorticity field is evaluated in wavespace, it is free of
first-differencing phase and amplitude errors.
(b) Recover the ir^al-space vorticity field u(x,t) by inverse
discrete Fourier transformation upon ft(k, t), namely
u)(x,t) = I g^tje1^'-
|k|
-------
(c) Evaluate, in discrete real-space,the vector cross
product term V(x,t) *w(x,t) - , say, £(x,t). This
product is, of course, fully aliased unless the pro-
cedure outlined in Appendix 3 is used.
(d) Generate the discrete Fourier transform of ^(Xft), de-
noted by V_(k't)
l(k,t) - I ^(Xftje"1-*- •
x
2
(e) Represent the diffusion term vV V(x,t) in wavespace as
2
- vk U(k,t) -
Then equation (6.6) may be written in prognostic spectral
form as
(6.11) (~ + vk2)U(k,t) - nk,t) - 6(k,t)A.
at — — — — —
The spectral form of the pressure head term, -V (P + -^ V-V) , need
not be included in (6.11), since the only effect of this term is
to guarantee that the velocity field remains incompressible. In-
stead, incompressibility can be enforced, at each successive time
step, by requiring the replacement of U(k_,t) by
(6.12) u' (Jc,t) = U(k-,t) - -=*-(k-U(k,t))
It
whence ik-U'(k,t) = 0, thus maintaining the incompressibility
(see section 2.1.3) . The nonlinear term V- (V0) in the ther-
mal equation (6.7) may be treated in a similar pseudospectral
manner, by the following procedure:
Ca) From the real-space velocity field V(x,t), and -the
real-space temperature field 6(x,t),, from the real-space
vector V(x,t)6(x,t) =, say, Z_(x,t).
-33-
-------
(b) Create the Fourier transform of Z^(x,t), which will
be denoted by C_(k,t):
Uk,t) - I Z^xjtje""1^ .
x
(c) Evaluate the divergence term V- (V9 ) , directly in wave-
space, as ik'C_(k.»t). (The divergence operation is ik-
in wavespace).
Thus, (6.7) may be written in spectral form as
(6.13) IE + Kk2)8(k,t) - - ikk?(k,t) - N2U3(k,t)
where U_(k,t) designates the Fourier transform of the vertical
velocity field w.
Leapfrog time differencing was used for the advective terms
and Crank-Nicolson time differencing for the diffusion terms.
Thus, if the superscript ()n designates the value of a field at
the n time step (t =nAt) the time evolution of the system can
be explicitly represented as
(6.14) (1 +
en(k,t)X]
(6.15)
- At [ik'
(6.16) Un+1(k,t) - Un+1(k,t) -
-84-
k
-------
Finally, we have included in our model a passive scalar
field, whose turbulent evolution is governed by the velocity
field v(x,t). The scalar conservation equation may be written
as
(6.17) || = vcV2c + s(X,t)
where c(x,t) is the scalar field, v is the scalar diffusion
coefficient, and s(x,t) is a source term for the scalar field.
If we let C(k,t) denote the Fourier transform of c(x,t),
S(k_, t) denote the Fourier transform of s(x, t), and M(k,t)
denote the Fourier transform of the product V(x,t)c(x,t),
then the temporally discretized version of (6.17), by analogy
with (6.15) is just
v k2At . v k2At
(6.18) (1 + _£_—)Cn+1(k,t) = (1 ^-= )Cn~i(k,t)
- At[ik-Mn(k,t) ]
+ At Sn(k,t).
The prognostic equations (6.14, 6.15, 6.18), together with the
replacement (6.16) define the numerical model. Initial and
boundary conditions will be discussed shortly.
-85-
-------
6-2 Physical Parameters
An attempt was made to design an experiment to simulate the
turbulent diffusion of a passive scalar field in a sheared, tur-
bulent flow, in such manner that reasonable comparison could be
made with the results of one or more turbulence closure models for
the same physical problem. To this end, the geometrical dimensions
of our numerical model were fixed as 9.6 km in each of the two
horizontal directions and 0.96 km in the vertical direction. The
horizontal and vertical grid intervals were 600 m and 60 m res-
pectively. All fields were taken to be periodic in the two hori-
zontal directions. At the upper and lower boundaries of the sys-
tem, free-slip boundary conditions were imposed on the horizontal
velocity components (that is U = V =0), while the vertical
Z Z
velocity W and the scalar field C were constrained to be zero.
For the initial experiments, we chose a neutrally stratified
flow, with no convective coupling between the thermal and vertical
velocity fields, because it was felt that the inclusion of the
stratification and Boussinesq forcing terms from the outset would
impose an unwanted degree of complexity in the interpretation of
the results. Accordingly, the terms At9n in (6.14) and N U3
in (6.15) were suppressed. Since the temperature field no longer
had any dynamical significance, equation (6.15) was in fact ignored.
A linearly sheared mean-velocity profile U(z) was chosen of the
form:
-86-
-------
The mean flow field was perturbed initially by a small-
amplitude, random three-dimensional field, the construction of
which is explained in Appendix 5. The mean shear profile U(z)
was maintained constant throughout the evolution of the flow by
fixing in time the corresponding spectral components 0^(0,0,K ).
Of course, fixing the mean wind field implies a continuous source
of kinetic energy for the perturbation wind field (by Reynolds-
stress conversions of the form (u'w'U ). It was hoped that the
z
turbulent perturbation wind field would come to approximate energy
equilibrium, the Reynolds-stress excitation being matched by suf-
ficiently vigorous viscous dissipation at the smallest scales.
(In fact, during each of a number of experiments over a range of
Reynolds numbers, the system never tended toward equilibrium in
any physically meaningful sense. At large Reynolds numbers, there
would result a strong cascade of energy toward high wavenumber,
which could not be dissipated sufficiently rapidly. If the Rey-
nolds number was reduced to the point where energy dissipation
could keep up with the generation of kinetic energy, the pertur-
bation flow structure typically degenerated into a set of large-
scale vortices ,iri.-a totally unphysical manner. Some examples of
these characteristic difficulties will be shown in a later sec-
tion. In retrospect, this behavior may be understood as a conse-
quence of the very limited modal resolution with which the experi-
ments were attempted, together with the strong mean^shear profile
adopted. )
-87-
-------
6* 3 An Unsuccessful Experiment
In this section we follow the evolution of the turbulent
flow fields in one experiment. The objective was, as previously
discussed, to create a turbulent flow structure for which the
generation of kinetic energy (by conversion of the mean shear U )
would come to equilibrium with the viscous dissipation of kinetic
energy. The difficulties encountered, which proved insurmountable
at present, are illustrated.
To aid in visualization, all velocity fields are pictured
in the form of contour plots. Positive-valued contour lines are
shown as solid lines, negative-valued contour lines are shown as
dashed lines. The minimum and maximum contour values and the
contour interval are indicated for each figure. These contours
have been constructed by an automatic contouring package.
In order to depict the three-dimensional spatial variation
of the velocity fields u, v, and w, the fields are displayed
in each of three orthogonal planes. That is, for example, we
show the variation of the u velocity component in the y-z
plane for some fixed x-coordinate, say x , the variation of u
in the x-z plane for fixed yn, and the variation of u in
the x-y plane for fixed z . The choices XQ, y z are
arbitrary; because the fields are turbulent, the individual con-
tours are to be understood only as representative of the charac-
teristic spatial structure in the given plane.
In the experiment to be discussed, the perturbation kinetic
energy spectrum was chosen to have the spectral power law depen-
1 222'1/2
dence E'(k)
-------
wavenumber. The rms magnitudes of the initial perturbation
velocity fields, compared to the rms magnitude of tho mean
velocity field U(z), were
= 2.2 x 10~2; (U')rms ~ 0.2 m/sec
?' i /•>
'
< (\7 ' } > ' — 9
_J -, p. - = 2.4 x 10 ; (V')rms - 0.2 m/sec
'
'
-3
= 5.6 x 10 ; (W')rms - 0.05 m/sec.
The ratio of total perturbation energy to mean energy was initially
about 5 x 10 . The viscosity coefficient used in the vertical
direction was fixed as 1/100 that of the viscosity coefficients
used in the two horizontal directions, consistent with the (square
of the) ratio of vertical to horizontal spatial intervals. The
2
particular choice made in this experiment was vtrrvmAT = 10 m /sec.
VVERTICAL = °"1 m /sec- (In other experiments we took values
of v which were tenfold higher and lower) . The equivalent large-
scale Reynolds number, based upon a characteristic velocity U=10
m/sec. , a horizontal length scale L = -^— km, and the value of
2rt
VHORIZONTAL' was thus ~ 160°- Tne characteristic eddy circulation
time L/U was ~ 10 sec. The time step, as prescribed by the
Courant advection condition for a spectral scheme,
UAt <
AX -
was taken as 10 sec., which is a slightly conservative choice
-89--
-------
Figure 6.1A U Field (y-z projection), t =0
Umin = 5'4 m/S6C
U =10.2 m/sec
max '
Contour Interval = 0.3 m/sec
-90-
-------
Figure 6.IB U Field (y-z projection), t = 10 sec
Umin = 5'4 m/sec
Umax = 10"2 m/sec
Contour Interval 0.3 m/sec
-91-
-------
Figure 6.1c U Field (y-z projection), t = 2 x 10 sec
U . =4.8 m/s'ec
mm
U =10.8 m/sec
max
Contour Interval 0.3 m/sec
-92-
-------
Figure 6.2A U Field (x-z projection), t = 0
°
min
m/sec
Umax = 10'2 m/sec
Contour Interval =0.3 m/sec
-93-
-------
Figure 6.2B U Field (x-z projection), t = 103 sec
Umin = 5'
Umax = 10'2 ra/sec
Contour Interval =0.3 m/sec
-94-
-------
Figure 6.2C U Field (x-z projection), t = 2 x 10 sec
U . =5.1 m/sec
mm '
U =10.8 m/sec
max '
Contour Interval 0.3 m/sec
-95-
-------
Figure 6.3A U Field (x-y projection), t = 0
U . =7.5 m/sec
mm '
U =8.35 m/sec
max
Contour Interval 0.05 m/sec
-96-
-------
Figure 6.3B U Field (x-y projection), t = 10 sec
U . =7.1 m/sec
mm '
U = 9.2 m/sec
max '
Contour Interval = 0.1 m/sec
-97-
-------
Figure 6.3C U Field (x-y projection), t = 2 x 103 sec
U . =5.4 m/sec
mm
U =10.2 m/sec
max '
Contour Interval 0.3 m/sec
-98-
-------
Figure 6.4A V Field (y-z projection), t = 0
>/" A. "" \
V . = -.4 m/sec
mm '
V = .8 m/sec
max
Contour Interval = 0.08 m/sec
-99-
-------
Figure 6.4B v Field (y-z projection), t = 103 sec
Vmin = -'48 m/sec
Vmax = '80 m/sec
Contour Interval = 0.08 m/sec
-------
Figure 6.4C V Field (y-z projection), t = 2 x 10 sec
_.•" i;- '.'.'.I' ;•"" i"" •:[ ••"
V . = -.64 m/sec
mm
V
max
m/sec
Contour Interval 0.08 m/sec
-101-
-------
Figure 6.5A W Field (y-z projection), t = 0
Wmin = ~'
Wmax = -1
m/sec
m/sec
Contour Interval = .01 m/sec
-102
-------
Figure 6.5B W Field (y-z projection), t = 10 sec
Wmin = ~'24
Wmax - -12
Contour Interval
= .02 m/sec
-103
-------
Figure 6.5C W Field (y-z projection), t = 2 x 10 sec
W . = -.30 m/sec
mm '
W = .27 m/sec
max '
Contour Interval .03 m/sec
-------
The typical structure of the initial velocity fields is
represented in figures 6.1A, 6.2A, 6.3A, 6.4A, 6.5A. Note in
particular that the mean vertically sheared wind field U(z)
is prominent in 6. 1A and 6.2A, although perturbed slightly by
the initial turbulent u' components. In figure 6.3A, which
represents an x~y slice of the initial u field, we see that the
horizontal structure of the u field is, indeed, initially quite
random, with a p1 fluctuating component of magnitude 0.4 m/sec
/
(the mean U value is about 8 m/sec). Figures 6.4A, 6.5A.show,
respectively, typical slices in the y-z plane, of the small-ampli-
tude velocity components V and w1 respectively. The structure
of the v> field is quite random spatially, while the w* field
is rather more organized into vertical rolls. (This is a conse-
quence of the specific procedure by which the perturbation fields
were synthesized; it reflects the fact that the vertical extent
of the system is one-tenth that of the horizontal extent.)
In figures 6.IB ... 6.5B we display the same fields after the
system was permitted to evolve for 100 steps (or 10 seconds),
corresponding roughly to the large-eddy circulation time. " The
mean shear structure of U(z) is still strongly evident, although
the u1 perturbation velocity field has more or less doubled in
amplitude, as pan be inferred from the range U -U . in
figure 6.3B. However, by comparing figure 6.4B with figure 6.4A,
and figure 6.5B with figure 6.5A, we see that the small-scale
structure of the perturbation velocity fields v1 and w' has
been damped appreciably, which suggests that, in this experiment,
turbulent energy is being dissipated at the higher wavenumbers
more rapidly than it is being cascaded downscale by nonlinear
-105-
-------
exchange with low-wavenumber modes.
After another 100 times steps of evolution, the situation
is worse yet. Although the perturbation velocity fields have
roughly doubled again in magnitude, this growth is almost en-
tirely, at the very lowest wavenumbers as may best be seen in
figures 6.1C and 6.5C.
Since, after two eddy "circulation times, the turbulent vel-
ocity fields were clearly further from equilibrium .than before,
it was pointless to continue the time evolution of this experi-
ment. We therefore started anew, in a second experiment, this
time with ten-fold, les.ser viscosity. Although the decrease in
viscosity tended to preserve relatively more of the small-scale
structure, the growth of total perturbation energy was unrealisti-
cally rapid, and. the contour plots showed the mean wind field was
soon masked by the magnitude of the perturbation velocity fields.
Alternatively, we also tried a ten-fold increase in viscosity,
which resulted in a rather more constant magnitude of perturbation
kinetic energy, but at the price of even stronger suppression of
all but the very lowest harmonic components. After further exper-
iments with somewhat difficult initial conditions, it was reluc-
tantly concluded that, with the limited resolution and cpu time
available, it would not be possible at this time'to evolve a
turbulent flow in reasonable equilibrium with the desired mean
shear profile.
-106-
-------
1'® Relaxation of Anisotropic Turbulence to Jsdtropy
One of the critical parameters of turbulent transport models
is the rate at which an initially non-isotropic flow decays to
its asymptotic isotropic state. Return-to-isotropy terms appear
exnlicitly, for example, in the turbulence closure model of
Donaldson (1973) and that Lumley and Khajeh-Nouri (1973), in
the form of the pressure-strain correlation expression
3u' .
pi I _ i. i
<
This represents> for incompressible flow, a redistribution of tur
bulent kinetic energy among the three velocity components. A
variety of closure prescriptions have been given for this term.
For example, Donaldson has chosen the Rotta model (Rotta, 1951)
- . 2
3uT 3u p^ _ ikq
/„ i n i _ - - ^ - \
(u u '
" X3x, -ax. ' A T V i k 3
k i 1
. 1/2
where q = (u .u1.) is a scalar measure of the turbulent
velocity field, Xi is a length scale introduced for dimensional
consistency (and requires independent specification), and p.. is
the density.
Instead of specifying the rate of return to isotropy in ad-hoc
fashion, 'one may have recourse to numerical simulation for an ex-
perimental measure of this parameter. It is important to use
direct spectral methods to ensure high accuracy. To perform such
a simulation, we require a fully three-dimensional closure-free
Navier-Stokes model. The basic model described in section 6.1
-107-
-------
serves especially well with the following minor modifications:
(a) the^ thermal and scalar concentration equations are
ignored;
(b) the vertical boundary conditions are modified to per-
mit general periodic symmetry conditions in the vertical
(rather than the particular cosine-like or sine-like
conditions previously specified).
To best examine the relaxation-ta-isotropy process, it is
useful to start from an axisymmetric turbulent flow, i.e., a tur-
bulent three-dimensional flow field with one axis of symmetry.
The initial flow conditions for such a relaxation problem may
be constructed in the following manner. We begin from a velocity
field V(x,t) which is statistically homogeneous and axially
symmetric about an axis whose direction is described by the
vector n. Suppose we define the second-order two-point moments
of V, namely:
<\7 I v +- ^ V I v' +• ' \ ~> — V I Y v1 +• -h M
<-v.v.xfujv.ix_,t. ) > — T . . ix/ x_ ,u,t ;.
The Fourier transform of ¥ is given by
¥. . (K,t,f ) = / exp (IK- (x-x1 ) Y . . (x,x' , t, t' )d(x-x' ) }.
13 — ' — — i j
For the axisymmetric case, this may be decomposed into
^ (K,t,t') = $(1) (K,t/t')e[1)(K)e.f1) (K)
+ $(2) (K,t,t')e(2)(K)e.|2) (K)
— i — j —
-108-
-------
whore the two vectors e and e_ are
e(1)(K) = (K x n)/| K x n|
e(2) (K) = K x(K x n)/|K x (K x n) | .
Here, e is the direction transverse to the axis of symmetry,
and e is orthogonal to both K and e/ '. The functions
<£>(1)(K) and $(2)(K) are given by
$ (1) (K) = <| (K x n) -U(K) |2>/|K x n|2
$(2)(K) - <|n-U(K)|2>/|K x n|2
where U(K) is the Fourier transform of V(x). The function
$ (K) represents the modal kinetic energy in the direction
(2)
transverse to the symmetry axis n, and $ (K) represents
the modal kinetic energy in the direction orthogonal both to
wave vector K and to n. If the statistical properties of
the initial velocity field are invariant to rotation about the
symmetry axis, then $ (K) and $'2'(K) are invariant as
well.
A useful 'experimental' measure of the rate of return to
isotropy is provided by the non-dimensional parameter
E2(t) - E1(t)
where
E, (t) = 2ir/ dK K2$(1) (K,t)
1 0
-109-
-------
JL ln{E2(t)- E^l)}
Figure 7.1 Time History of R (t) for E.. (0) = 0 and
E^O) = 1/2 E2(0)
-110-
-------
,-<» 2 (7 )
E (t) = 2-n I dk K $v '(K,t)
0
represent the kinetic energies in the orthogonal directions
e(1)(K), e(2)(K), respectively, E(t) = E-^t) + E2(t) is
the total kinetic energy, and
e(t) - 4TT2 /°°dk K4[${1)(K,t) + $(2)(K,t)]
0
is the total viscous energy dissipation rate. The ratio
E(t)/e(t) describes the rate of fractional decay of total
energy, while the ratio -r—
E2 (t)-E1(t) |/[E2(t)-E1(t)| describes
the time-dependent fractional rate of return to isotropy. R (t)
a.
may be considered equivalent to the Rotta constant.
To study the time-dependence of this parameter upon the extent
of anisotropy initially present in a turbulent flow, we performed
two independent experiments. In the first case, the axisymmetric
turbulent velocity spectrum was chosen so that initially all the
(2)
modal kinetic energy was in the direction e_ (K) , thus E (0)=0,
E (0)^0, and the initial flow field was totally anisotropic. In
the second case, the initial conditions were chosen as E (0)=yE2(0)
so that the flow was fairly isotropic to begin with,
Rotta1s phenomenological argument suggests that R (t) should
ci
be constant in time with a magnitude ~ 2. However, this argument
is predicated upon small departures from isotropic conditions. It
is interesting, therefore, to see how well this prediction holds
for flows which are initially quite anisotropic. Figure 7.1 dis-
plays the time history of the expression R (t) for each of the
cl
two cases studied (the time t is normalized by the large-eddy
circulation time). Somewhat surprisingly, in view of Rotta's
-111-
-------
argument, we sec that for the second case E1(0) = ^ E2(0),
which represents fairly isotropic initial conditions, the rate of
return to isotropy is considerably less constant than that observed
in the first case E, (0)=0, which began with completely anisotropic
initial structure. This is indicative of large statistical fluc-
tuations in the first case, and suggests that even fairly iso-
tropic turbulent flows remain in strong disequilibrium for long
periods of time. That is, isotropic relaxation of turbulence is
a very slow process. Thus, turbulence closure models which suppose
self-similar asymptotic conditions are strongly suspect in this
regard.
We may note, finally, that these numerical results are roughly
consistent with calculations based upon the Direct Interaction Ap-
proximation turbulence theory of Kraichnan (Herring 1974).
-112-
-------
8.0 Recommendations for Future Research
Although the primary research objective of direct turbulence
simulations was not achieved in the current contract because of
limited computer resources, we remain confident that the use of
spectral methods holds considerable promise for operational
improvement of air quality numerical modeling in a number of
specific applications areas. These include:
(a) more accurate treatment of advection and viscous terms
by means of spectral representations which preserve
better than do finite-difference schemes the true phase
and amplitude structure of such terms;
(b) accurate representation of the advection-diffusion of
point-source and line-source pollutants, as suggested
by the discussion in chapter 5>
(c) efficient spectral schemes for accurate representation
of vertical structure and boundary layers, in particular
the use of Chebyshev spectral methods to provide effec-
tively increased resolution at viscous boundaries;
(d) the use of spectral methods, in conjunction with mapping
techniques, for accurate treatment of topography, as,
for example, the modeling of pollutant advection-diffusion
from a highway situated in a valley;
(e) the use of spectral methods for the accurate simulation
of lateral inflow and outflow boundary conditions, with-
out requiring the imposition of ad-hoc one-sided differ-
encing schemes.
We recommend, therefore, that careful investigation be made
of the prospective benefits to be gained by implementation of
these spectral methods within existing air-quality models of
-113-
-------
importance to EPA, in particular, the photochemical air-pollution
model developed at Systems Applications, Inc., and the air-pollu-
tion model developed at IBM.
As to the problem of simulating the turbulent diffusion of
passive contaminants in fully three-dimensional turbulent sheared
flows by the direct spectral approach of the current effort, we
do not recommend further pursuit of this study at present for two
reasons.
The first of these concerns the timeliness of application of
this research. It does not appear likely that extension of the
present study would provide to EPA, in the near term, the capability
for better estimation of eddy Austauch coefficients for turbulent
diffusion of pollutants in the place of field observations. We
do believe, however, that the prospects for success of a longer-
term research effort remain very good. We do not know of any
fundamental theoretical impediments to the successful numerical
simulation of such turbulent flows by direct spectral methods,
especially in the light of previous recognized achievements in
this area by Orszag and collaborators.
The second reason concerns the prohibitive economics of the
required calculations. Neither the computer core time necessary
for adequate spectral resolution, nor the computer time required
to establish realistic turbulent dynamics, will be available to
us in the near future. At present, such calculational efforts
can possibly only be undertaken by university-affiliated groups
with access to the NCAR computer facility or a comparable installa-
tion. When the computing-resource situation improves, as may be
anticipated, it would be most worthwhile, we believe, to pursue
this study anew.
-114-
-------
9.0 References
Arakawa, A. (1966) J. Comp. Phys. _1, 119-143.
Arakawa, A. (1970) Numerical Solution of Field Problems in
Continuum Physics 24,, American Mathematical Society,
Providence, R.I.
Bourke, W. (1972) Mon. Weath. Rev. 100, 683-689.
Calder, K. L. (1967) QJRMS, 88-92.
Collatz, L. (1960) Numerical Treatment of Differential Equations,
Springer-Verlag, Berlin.
Cooley, J. W. and Tukey, J. W. (1965) Math. Comp. 19, 297-301.
Crowley, W. P. (1968) Mon. Weath. Rev. 96, 1-11.
Cunnold, D., Alyea, F., Phillips, N. and Prinn, R. (1975) J.A.S.
32^ 123.
Deardorff, J. W. (1970) J. Fluid Mech. 41, 453.
Donaldson, C. duP. (1973) "Atmospheric Turbulence and the Dis-
persal of Atmospheric Pollutants," EPA report EPA-R4-73-016a.
Dorr, F. W. (1970) SIAM Rev. 12, 248.
Eliassen, E., Machaver, B. and Rasmussen, E. (1970) "On a Numer-
ical Method for Integration of the Hydrodynamical Equations
with a Spectral Representation of the Horizontal Fields,"
Institute for Teoretisic Meteorologi, Univ. of Copenhagen,
Report No. 2.
Fornberg, B. (1972) "On High Order Approximations of Hyperbolic
Partial Differential Equations by a Fourier Method." Dept.
of Computer Sciences, Uppsala Univ. Sweden, Report No. 39.
Fox, L., and Parker, J. B. (1968) Chebyshev Polynomials in Num-
erical Analysis, Oxford Univ. Press, Oxford
Fromm, J. E. (1963) "A Method for Computing Nonsteady Incom-
pressible Viscous Fluid Flows," Los Alamos Scientific
Laboratory, Report LA-2910.
Fromm, J. E. (1970) An Introduction to Computer Simulation in
Applied Science, ed., F. F. Abraham, W. A. Tiller. IBM
Data Processing Division, Palo, Alto, Calif.
GARP (1974) "The GARP Programme on Numerical Experimentation,"
Report of the Int. Symp. on Spectral Methods in
Report No. 7, Copenhagen.
-115-
-------
Gentleman, W. M. and Sande, G. (1966) Proc. Fall Joint Computer
Conference, 1966, 563-578.
Gordon, A. and Stern, W. (1974) "Spectral Modelling at GFDL,"
on GARP (1974) (see above).
Grammeltvedt, A. (1969) Mon. Weath.Rev. 97, 384.
Harlow, F. H. and Welch, J. E. (1965) Phys. Fluids 8, 2182.
Herring, j. R. (1974) Phys. Fluids 17., 859-872.
Herring, J. R., Orszag, S. A. Kraichnan, R. H. and Fox, D. G. ,
(1974) J. Fluid Mech. ^6_, 417-444.
Hockney, R. W. (1970) Methods in Computational Physics 9, 135.
Israeli, M. (1971) "Time Dependent Motions of Confined Rotating
Fluids," PhD. Thesis, Dept. of Applied Mathematics, M.I.T.,
Cambridge, MA.
Kreiss, H. O. and Oliger, J, (1972) Te&lus 24, 199.
Kreiss, H. O. and Oliger, S. (1973), Methods for the Approximate
Solution of Time Dependent Problems, GARP Publication Series
No. 10, WMO, Geneva.
Lanczos, C. (1956) Applied Analysis, Prentice Hall Pub. Co.,
Englewood, N.J.
Lilly, D. K. (1965) Mon. Weath. Rev. 93, 11.
Lumley, J. L. and Khajeh-Nouri, B. (1973) "Computational Model-
ing of Turbulent Transport," Proc. Second IUGG-IUTAM
Syrap. on Atmoshperic Diffusion on Environmental Pollution,
Academic Press, N.Y.
Merilees, P. E. (1973) Atmosphere rL, 13-20, (1974) Mon. Weath.
Rev. 102, 82-84.
Metcalfe, R. W. (1973) "Spectral Methods for Boundary Value
Problems in Fluid Dynamics," PhD. Thesis, M.I.T. Cambridge
MA.
Molenkamp, C. R. (1968) J.A.M. T_, 160-167.
Neuringer, J. L. (1968) SI AM J. Appl. Math 16_, 4.
Orszag, S. A. (1969) Phys. Fluids (Suppl. 2) 12_, 250.
-116-
-------
Orszag, S. A. (1971a) Stud, in Applied Math 50, 293.
Orszag, S. A. (1971b) J. Fluid Mech. 49, 75
Orszag, S. A. (1971c) Stud, in Applied Math 50, 395.
Orszag, S. A. (1971d) Phys. Rev. Letters 26, 1100.
Orszag, S. A. (1971e) J. Fluid Mech. 50, 689.
Orszag, S. A. (1972) Stud, in Applied Math 51, 253.
Orszag, S. A. and Israeli, M. (1974) Ann. Rev. Fluid Mech. _6_, 281,
Orszag, S. A. and Pao, Y. H. (1971) "Numerical Computation of
Turbulent Shear Flows," Advances on Geophysics ISA, Academic
Press, Inc., NYC 225-236.
Orszag, S. A. and Patterson, G. S. (1972) Statistical Models
and Turbulence, ed. Rosenblatt, Springer-Verlag,Berlin, 127.
Patterson, G. S. and Orszag, S. A. (1971) Phys. Fluids 14,
2538-2541.
Phillips, N. A, (1959) The Atmosphere and the Sea in Motion.
Rockefeller Institute Press, NYC.
Piacsefc, S. A. and Williams, G. P- (1970) J. Comp. Phys. 6, 392.
Roache, P. J. (1972) Computational Fluid Dynamics, Hermosa
Publishers, Albuquerque, NM.
Rotta, J. C. (1951) Z. Phys. 129, 547.
Spiegel, E. A. and Veronis, G. (1960) Astrophys. J. 131, 442-447.
Swartztrauber, P. (1972) A Direct Method for the Discrete Solu-
tion of Separable Elliptic Equations (unpublished).
Williams, G. P., (1969) J. Fluid Mech. 37, 727.
-117-
-------
Appendix 1
The Data Content of Discrete-Grid Representations
It is conventional practice / in the numerical solution of
many fluid motion problems, to perform such calculations with
respect to a discrete spatial grid and suitably spatially dis-
cretized versions of the governing differential equations whether
by finite-difference of finite-spectral methods. In this section,
we discuss the underlying assumptions and limitations which inhere
to the discretization method.
Imagine that we wish to estimate the frequency spectrum of a
continuous spatial waveform, f(x), which is periodic with wave-
length L. A harmonic analysis of this waveform might be performed
as follows :
(1) partition the interval L into N equal segments by
the discrete points x =pAx (p=0, ... N-l) ; Ax=L/N,
where for convenience, N is a power of 2);
(2) tabulate the string of values f (,x ) at each of
these points;
(3) employ the Discrete Fourier Transform algorithm
(Appendix 2) to compute the Fourier modes 9(^q)
corresponding to the discrete wavenumbers
9-°' i- ••• N/2 •
Note that there is an upper wavenumber cutoff (the "folding" or
Nyquist, frequency N/2) so that we possess no data whatever about
-118-
-------
harmonic components higher than q = - N/2. If we now attempt to
reconstruct the waveform f(x) by inverse discrete Fourier trans-
form upon the truncated set of modes g(A ), q=0, ... 2' the only
data recoverable about the spatial behavior of f(x) is precisely
its values at just the discrete points x (p=0, ... N-l). That
is, the data content of the truncated discrete spectral represen-
tation of the waveform is precisely equivalent to that contained
in the N distinct functional values f(x ); P=0, ... N-l. There-
fore, the spectral representation of this data contains precisely
N distinct degrees of freedom. To see why this is so, consider
the Fourier series representation of a real function periodic with
fundamental wavelength L, sampled at the uniformly spaced points
x (p=0, ... N-l), viz:
f V = aO + q^ Uq C03^ + bq Sin^) '
If this system of equations is to be uniquely soluble for the set
of Fourier coefficients a's and b's, then we can only specify
a total of N Fourier coefficients. The most general choice of
coefficients which are determinable for the data string f (x ) ,
... f(xNt_-i) in the absense of any additional information about
the function f, is precisely the set
V V V a2' b2 ••••
Note that the discrete Fourier representation of the function
incorporates exactly as much data about the spatial structure
of the function f (x) as does the discrete string f (x ) ... .
How then do we interpret the correspondence between a real
physical variable which is a continuous spatial function and
-119-
-------
its point-wise representation in a discrete-grid numerical model?
One approach is to assume that the physical variable is sufficiently
smoothly varying; that is, its assigned value at a discrete point
x is representative of its behavior over some neighborhood K-XO
of the point in question, although the notion of neighborhood is
not made precise. In fact, for highly discontinuous processes,
as, for example, in the propagation of shocks, or small-scale tur-
bulent cascade processes such an assumption is extremely dangerous.
A more satisfactory approach is to consider the discrete-point rep-
resentation as only a " best-fit" to the data contained within the
numerically resolvable scales of the variable. The finite cutoff
wavelength forces us to surrender any detailed data about small-
scale structure which cannot be determined entirely from the re-
solvable larger scales. This interpretation does not require a
notion of smoothness.
-120-
-------
Appendix 2
/
Discrete Fourier Transform Representations
The bulk of the computational effort in the model we have
'developed is expended in performing discrete Fourier transforms
from discrete three-dimensional real space to discrete three-
dimensional wave space, and vice versa, using the fast Fourier
Transform (FFT) as discussed in Appendix 4. In this appendix
the discussion is limited only to the discrete three-dimensional
Fourier transform. The discrete three-dimensional Fourier trans-
form pair is defined as follows:
Let f be a three-dimensional field defined over the
Jl' J2' J3
N, x N_ x N real-space grid:
J2=l, ... N2
J.,=l , ... N_ .
Let F be a three-dimensional field defined over the
T' 2 ' i
NI x N x N_ discrete wave-space grid
(kn, k_, k_) = (K,Ak K^Ak K Ak) K =1, ... N
L £• J L - I £• i 3 1- -L
K =1, ... N .
Then f and IT - „ are a discrete Fourier trans-
Jl" J2' J3 Kl' K2' 3
form pair if (with one choice of normalization)
. N, N0 N,, 9
1} FR. R = i r y1 y2 y3 e"27T
Nl
(cont'd)
-121-
-------
(A2.1)
fcont'd)
(J2-
1)(K2-1)
N2
(J -1) (K -1)
+
J
J
(A2.2)
J-i , J,, , J
(
+
N N
yl y2
I L
3 K^=l K2=l
Kl-1) (J -I)
+
N2
?3 e+2,,1
z
K'=l
(K^-l) (J1-l)
Nl
(Kl-1) (Jl-1)
3 3
N !<[, K' K'
The basis functions are orthogonal in the sense
Na +0 .(J -1)(K -K1)
y ± 2 TT i a a a
j' = l
a
N
a
N 6M (K ,K'
a N a' or
a
(a=l,2,3)
(3=1,2,3)
wh e re
'l, for modulo[(K-K1),N]=0
0, otherwise
delta function, modulo N.
Writing the functions in vectorial form
, i.e., the Kronecker
J = (J,, J2, J3)
K = (K;L, K2, K ) ,
-122-
-------
wo have
2-iri j. (K-K1 )
Nl N2
If f
1, 2,U3
is a real array, then F,
is a conjugate-
symmetric array (the conjugate symmetry being denoted by a star)
FK1 K K, = |FK, K, K, 1 , where
1,2,3 I ^1' 2' 3
K!=l if K.=l
if
Noting that the wave-space coordinates K^, (i = l,2,3) are indexed
relative to one rather than zero and folded relative to the Nyquist
frequency N/2, the following table provides the correspondence
between the computational wave-space coordinate K. of the algorithm
and the associated non-dimensional physical wavenumber:
K. of DFT algorithm
1
2
corresponding physical
wavenumber
0
1
N/2
N/2 + 1
N/2 + 2
N/2-1
N/2
-(N/2 -1)
N
-1
-123-
-------
Appendix 3
Elimination of Aliasing in Discrete Real-Space Multiplications
This appendix describes, in a simplified two-dimensional con-
text, the algorithm developed by Orszag (1971a) for the efficient
alias-free evaluation of real-space products. In order to help
the reader become fully familiar with such techniques, we have
included all relevant mathematical details.
Let us define the two-dimensional real-space grid (Xj , Yj )
-L «£
for 0 _< J , J <_N-1; that is, X = J^^x, y = J Ay, Ay = Ax.
For compactness we will use the vector notation J=(J1,J?).
Similarly, we define the two-dimensional wave-space grid K=(k ,k ),
where - — _< k , k < j. These are equivalent to the notation of
Appendix 2, but for a shift in coordinate origins. Consider now
two real functions of the variables X , Y These will be
Jl J2°
denoted as
gT = g(X , YT )
— 1 9
hj = h(Xj , Yj )
defined over the real-space grid. Their product
W(XT YT ) = w, = q h
1' 2 — 3 ^Z
is fully aliased in the sense that all spectral components | KJ >'•_- appec
falsely, within the spectral domain |K| £ j. The problem is to
construct a scheme for evaluation of W such that only those
spectral components for which |K| < l^- are properly included.
-124--
-------
Define, as in Appendix 2, the two-dimensional discrete
ourier transform pair(1'2).
(A3.1) f = £ p e——- K-J = say, DFT,
«^ „ i\ JM — ~-~ J\
|K| <_ N/2
where £ , |K| ^ N/2 is taken to mean - ^ < K. , K < ~
K 2-122
The inverse DFT is
, _ 2-rri K-J
(A3. 2) F = -L- Jf e ~N = o'1 fT
K 2 '±. J J
— N J — —
where 7 is taken to mean 0 < Jn , J_ < N-l .
j -12-
Then the DFT of the fully aliased local product W,. is
j
27TJ K-J
g h e
N j
W = D W = -^ I g h e N
2
2iri K1 -J , 2iri K"-J _ 2TTJ K-J
(A3. 3) = -~ J I G_ e + N I H „ ie N e N
N2 J K- K K" ^
By rearrangement of order of summations we have
,-,„. w lVVr. „ v -2TTi(K-(K'+'K")) -J
(A3. 4) W = -^ I I G , H „ I e -^- - - - -
- N K' K" - - J
where
GK =
HK = ° hJ '
(1) The notation K-J is understood to mean the scalar dot product
K1J1+ K2J2 + K3J3-
(2) The notation D is shorthand for the inverse Fourier transformation
operation applied to PK-
-125
-------
The orthogonality relationships can be written as
N-l 2-ni (K -(K' + K")) J
~~
(A3.5)
J1=0
Since
N -1- -1- •*• -1- = N, if K = (K' + K") modulo N
= 0, otherwise .
< ^-, IK" I < 5L then IK' + K: I < N and the
j. • — 2 ' 1 ' — I i z ' —
summation (A3.5) - N, if K± = (K| + Kj) or (KJ + KJ + N), or
(K' + K" - N). Similar orthogonal relationships hold for summa-
tion over index J?. Therefore, the summation over index J^.
Therefore, the summation
\ ^L It J_ * T_ , . ii % \
e N v- v-
2 r
N , if
H
Ki
Ki
Ki
-r j\ /
+ K" =
+ KJ =
+ K" =
( • J
Kl
K -N
K +N
K2
and K'
K2
+ K2 -
+ K'2 =
+ K2 =
,
K2
K2-N
K +N
^-- >
and
= 0, otherwise
(A3.4) may therefore be written out fully as
+
K2
K2
-126-
-------
whence
(A3. 6) WT, T, = W,
-KUK +N
where
WK K = G H
a' b K
The first term on the right hand side of (A3. 6), W „ „ , is the
K-i- t K.p
— 1
desired alias-free convolution of GT, and H^, i.e., the D
K — K
transform of the alias-free product of g_ and h . Each of
- - j j^
the other terms in A3. 6 implies the presence of aliased contri-
bution of the form K ±N , etc. The task now is to isolate
the first term from the other eight terms on the right hand side
of (A3. 6). To achieve this in economical fashion, Orszag has
introduced the "shif ted-grids" procedure. There are a number of
variants of this method, but we will choose the "quartergrid"
scheme as particularly convenient. Let us define the four shifted
grids
s S
(p=lf "4)
where
^ ='(J]
-P
+ %
' J2 '
S2
•• -/p)
,
P P
S = (1,1);
-127-
-------
(T'hc shifted grids are introduced for the purpose of developing
expressions which permit the isolation of the alias-free term.)
For each shifted grid (p=l,--4), define
2TTJ K-J
(A3.7) g ) = GKe N Sp
2TJ K-J
(A3.8) h,T = I He N -p
^q ' V *
O XV
-p -
and define
(A3. 9) W^ = D'1 g -h( )
K1'K2 (-S ) -S
-P -P
= -^ I I I I GK, K, H „ I I exp <^i) [(KJ+K--K )
N^ K| K^ K^ K2' Kl' K2 Kl' 2 J1 J2 N 111
Again using the orthogonality relationships
I I exp
S!
J, J. ^ N
1 "2
S2 Sl S?
^— T ^ - -L_ £•
.- ^ -2 -2, ,_2- 4 , _, „ c N
where
K* , K* = 0, ± N only
Thus, evaluating the quadruple sums ^ ... £ for each of the
Ki
shifted grids in turn, we find:
-120-
V II
K2
-------
(A3.10a)
§1
KlfK2
-1 W
KlfK2-N
+ i W.
+ W
(A3.10b)
W
K ,
1W
K1-N,K2
iW,
(A3..10C)
WT
, «
- w
(A3.10d)
-4
+ i W.
-129-
-------
(A3.10d) (Cont'd)
- i W
+ i W
— W
K1+N,K2+N
Thus :
1, 2
'^2
S _
+ w~4 = 4v\r
K1'K2 Kl, K2
whence we have the remarkable result
(A3.11)
W
K1>K2
= T I
W,
p= 1 , • • 4
Finally, by performing one further DFT
(A3.12)
W.
= DW
K1'K2
we recover the alias-free product of g and h .
J J
As outlined above, the grid-shifting procedure requires
twelve DFT operations, namely
K
-P
H
K
-P
(p=l, •• 4)
W
~p
-130-
-------
Orszag shows that it is possible to save half of these operations
by a simple modification to the above scheme, as follows:
(a) use only one diagonal pair of shifted grids, say
J and J ;
-1 -2
(b) perform evaluations (A3*7, 3.8);
(c) replace (A3.11) with
2 s
W1 = - I W~p
K1'K2 2 p=l K1'K2
where the prime signifies that W1 is not yet the alias-
1' 2
free convolution of G^, and H ;
K K
(d) truncate the field W to the wave-space region
0 £ (K2 4- K2
and let
W = fw1 )
KlfK2 V KX,K2'TRUNCATED'
(e) obtain the alias-free real-space product field by
performing operation (A3.12).
The corresponding algorithm for three-dimensional systems
is considerably more complex in implementation, but is concep-
tually identical.
-131-
-------
Appendix 4
The Fast Fourier Transform
The basic algorithms underlying the Fast Fourier Transform
are well documented in the literature (Gentleman and Sande 1966,
Coo ley and Tukey 1965) .
The definition of the discrete Fourier transform (DFT) de-
fined over the one-dimensional string of N data points D. (j^O,
. . . , N-l) can be written as
i ^~ ^ ' v ? /N
(A4.1) T^ = — I D. WDK , W=e^ 17 , k=0 , ..., N-l.
K /N" j = 0 J
(where we have incorporated /N normalization in each direction,
in contrast to the one-way 1/N normalization used in Appendix 2)
Explicitly for N=8,
/8 TQ = D W° + D-j^W0 + D2W° + . . .
/8 T, = D.W0 + D..W1 + DW2 + ...
, . .. 0
(A4.2) -1 ° 1 2
A «
Q
Q
A « »
/8 T = DW + D-j^W + D2W
/8 T = DW° + DW3 4- DW6
Naively, to perform the evaluation of T , one would compute
K
directly from (A4.2). This would require 64 complex multiplications
and 56 complex additions. However, the trigonometric function
T,jk 2irijk/N 2-irjk , . . 2-rrjk ^ . ,_
WJ = e J - cos — rj* — + i sin — ^ — possesses two important
properties namely:
(1) Periodicity ... WN - e2TTiN/N = 1^ SQ thafc fQr N=g^
W8 = W° = 1, W9 = W, W10 - W2, etc.
(2) Symmetry ... For N a power of two, we may utilize the
-132-
-------
reflections of the cosine and sine functions in the
2345
eight octants: W = i, W - iW, W = -1, W = -W, etc.
Hence, expressions (A4.2) reduce to 32 multiplications and
56 additions:
(A4.3)
/8 TQ = DQ + DI + D2 + D3 + D4 + D5 + DS + D?
/8 T, = D_ + D,W + D^i + D_iW - D. - D._W - Dri - D_iW
10123 4567
/8 T2 = DQ + DI± - D2 - D3i + D4 + D5i - D6 - D?i
/8 T- = D + D iW - D,i + D3W - D - D iW + D i - D_W
/8 T4 = DQ - DI + D2 - D3 + D4 - D5 + D6 - D?
/8 T = D - D-j^W + D i - D3iW - D^ + D^ - DQi - D iW
^ T6 = °0 - °li - °2 + °3i + °4 - V * °6 + °7i
/8 T = D - D iW - D i - D_W - D. + D-iW + D i + D W .
/ U J- ^-•J^tJ O /
By adroit grouping of terms in (A4.3), it is possible to
reduce the amount of computation still further. Proceeding sys-
tematically, define:
A0 = DQ//8
Ax = D4//8
A2 = D2//8
A3 = D6//8
A4 - D1//8
A5 = D5//8
-133-
-------
(A4.4) (Cont'd)
A6 = D3//8
- D?//8
Next:
(A4.5!
Bo = Ao + Ai
Bi = Ao - Ai
B3 = (A2-A3)i
B4 = A4 + A5
B5 = A4 - ^
B6 = A6 + A7
B? = (A -A?)i
C0 = B0
C2 = B0 - B2
C3 = Bl - B3
C4 = B4 + B6
C6 = 'W1
C? = (B5-B?)iW
= C + C
0 0 4
= r + c
1 Cl °5
T3 = C3 + C7
T4 = CQ * C4
-134
-------
(A4.5) (Cont'd)
T —
6
T =
- C
C
— C
C
Excluding the normalization of (A4.4), we have accomplished
the DFT in only 5 complex multiplications and 24 complex additions.
In general, when N is a power of two, the Fast Fourier Transform
(FFT) performs a DFT in Only j log2N complex multiplications
(for N >_ 32) and only N log0N complex additions. Comparatively,
.£ ~
2 2
the naive approach to (A4.2) requires N multiplications and N
additions.
The flow diagram below displays the FFT technique of (A4.4)
and (A4.5). A circle of the form It/ y means to add (subtract)
all the values flowing in from the left and then to multiply by
w". An empty circle means do nothing to the incoming value.
-------
Notice that the first stage of the calculations, (A4.4),
D'Tforms a peculiar permutation of the data.
Index before Index after
0 0
1 4
2 2
3 6
4 1
5 5
6 3
7 7
Writing these indices in binary notation, we see why this
process has received the name "bit reversal":
Index before Index after
000 ' 000
001 100
010 010
Oil 110
100 001
101 101
110 Oil
111 111
This sorting is required for all the major variations of
the Fast Fourier Transform algorithm. It can be absorbed into
the computational stages of the FFT, but only at the cost of
doubling the amount of computer storage. The time required for
this permutation is generally about 25% of the computational
stage time. Thus, it is normally better to trade time for
storage and perform bit reversal explicitly. (Notice that
no additional storage is needed to perform bit reversal, since
all data are exchanged in pairs.)
It is easy to generate the sequence of integers 0. 1, 2,
etc., each from the previous one. In binary notation, the process
-136-
-------
may be stated as: add one to the low order bit, and carry upward
if it overflows (becomes larger than one), repeating with higher
order bits until no further overflows occur. Similarly, beginning
with 0, which is certainly the bit reversal of the integer 0, we
generate the reversal of the successive integers by an analogous
process: add 1 to the highest bit of the reversed index; if it
overflows, carry one to the next-lower bit, repeating with lower
order bits until no further overflow occurs.
E.g., for N=8, we start with 000 and add one to the high-
order bit:
100 = 4.
Generate the next reversed integer by adding a high order one:
200+010 = 2.
Again add one:
110 = 6.
And again
210^020-^001 = 1.
And so forth.
In a high-level language like Fortran, an integer may be
tested for the value (0 or 1) of the high-order bit by comparing
it to N/2; the high bit is 1 if the integer is greater than N/2.
If the high bit is 0, then we must test the second highest bit
by comparing to N/4, and so on.
-137-
-------
Appendix 5
Construction of Initial Perturbation Flow Field
The initial flow fields specified in this model, as described
in chapter 6, are of two types:
(a) the mean velocity field, with a suitable shear structure
U(z) ,
(b) the perturbation velocity field, which is comprised of
a large ensemble of small-amplitude incompressible modes,
distributed in statistically homogeneous and isotropic
manner, and obeying a designated spectral power law E(k).
This Appendix describes the process by which such a per-
turbation field is constructed in wave space.
Let us define the two random processes = <8 . (k ) B .(k1' )>=0
i 3 — i — ] —
(kyk1) where <••> designates ensemble averaging, and k is the
discretized wavenumber vector. The random processes a and 3
can be obtained from any of a variety of pseudo-random number
generation algorithms. Next, define the vector potential field
in wave space as
A(k) = { [E(k)]1/2/k2} + | 2 SLn a (k) |1/2 e27TlS (-} .
Then, we construct the initial perturbation velocity field
in wave space U(k) as the curl, of A(k)
U(k) = ik x ,A(k)
-138-
-------
which automatically guarantees incompressibility- Finally, the
real-space perturbation velocity field V(x) is obtained from
U(k) by discrete Fourier transformation.
-139-
-------
TECHNICAL REPORT DATA
(Please read Instructions on the reverse before completing)
REPORT NO!~
EPA-600/4-76-007
2.
4. TITLE AND SUBTITLE
SPECTRAL MODELING OF ATMOSPHERIC FLOWS AND
TURBULENT DIFFUSION
8. PERFORMING ORGANIZA IION CODE
7. AUTHOR(S) " " '
Arthur Bass, Steven A. Orszag
8. PERFORMING ORGANIZATION REPORT NO.
3. RECIPIENT'S ACCESSION-NO.
6. RFt'OHT DATE
January 1976
9. PERFORMING ORGANIZATION NAME AND ADDRESS
Flow Research, Inc. (N.E. Division)
1 Broadway
Cambridge, MA 02142
10. PROGRAM ELEMENT NO.
1AA009
11. CONTRACT/GRANT NO.
EPA 68-02-1297
12. SPONSORING AGENCY NAME AND ADDRESS
Environmental Sciences Research Laboratory
Office of Research and Development
U.S. Environmental Protection Agency
Research Triangle Park, North Carolina 27711
13. TYPE OF REPORT AND PERIOD COVERED
Final
14. SPONSORING AGENCY CODE
EPA-ORD
15. SUPPLEMENTARY NOTES
16. ABSTRACT
This report presents a survey of discrete spectral and pseudospectral
numerical methods to simulate atmospheric flow and turbulent diffusion.
Some applications of these methods to air quality simulation modeling are
presented. A three-dimensional spectral incompressible numerical model is
described in detail. Computational resource limitations precluded successful
evaluation of eddy Austauch coefficients. Some numerical results are presented
for the rate of relaxation of anisotropic flows.
Recommendations and suggestions for further research are made concerning
the prospective utility of these spectral methods for air quality simulation
modeling.
17.
KEY WORDS AND DOCUMENT ANALYSIS
DESCRIPTORS
b.lDENTIFIERS/OPEN ENDED TERMS
c. cos AT I Field/Group
* Numerical analysis
Turbulence
* Turbulent diffusion
Air flow
Air pollution
Environmental simulation
Mathematical models
12A
20D
04A
13B
14B
8. DISTRIBUTION STATEMENT
RELEASE TO PUBLIC
19. SECURITY CLASS (This Report)
UNCLASSIFIED
21. NO. OF PAGES
151
20. SECURITY CLASS (This page)
UNCLASSIFIED
22. PRICE
EPA Form 2220-1 (9-73)
------- |