EPA-650/4-74-044
NOVEMBER 1974
Environmental Monitoring Series
'••«
'"•*'•*«
'"• *V,
>•'•*•
'•'•".
>••*.•,
*•"•».%
*••'.*»
"•"•'•
*.
•V •
,V-',
•*•'.'.
:j$:j:bjfj:j3f$±§ji:R^
c^^A^^t^jjrjillic-ciJEjpjcii-l^ji-^^ft];^
-------
EPA-650/4-74-044
LABORATORY
AND NUMERICAL SIMULATION
OF PLUME DISPERSION
IN STABLY STRATIFIED FLOW
OVER COMPLEX TERRAIN
by
Jung-Tai Lin, Hsien-Ta Liu, Yih-Ho Pao,
Douglas K. Lilly, Moshi Israeli, and Steven A. Orszag
Flow Research, Inc.
1819 S . Central Avenue
Suite 72
Kent, Washington 98031
Contract No. 68-02-0800
ROAP No. 21ADO
Program Element No. 1A1009
EPA Project Officer: William H . Snyder
Meteorology Laboratory
National Environmental Research Center
Research Triangle Park, North Carolina 27711
Prepared for
OFFICE OF RESEARCH AND DEVELOPMENT
U.S. ENVIRONMENTAL PROTECTION AGENCY
WASHINGTON, D.C. 20460
November 1974
-------
EPA REVIEW NOTICE
This report has been reviewed by the National Environmental Research
Center - Research Triangle Park, Office of Research and Development,
EPA, and approved for publication. Approval does not signify that the
contents necessarily reflect the views and policies of the Environmental
Protection Agency, nor does mention of trade names or commercial
products constitute endorsement or recommendation for use.
This document is available to the public for sale through the National
Technical Information Service , Springfield, Virginia 22161.
11
-------
LABORATORY AND NUMERICAL SIMULATION OF
PLUME DISPERSION IN STABLY STRATIFIED FLOWS
OVER COMPLEX TERRAIN
by
Flow Research, Inc., Kent, Washington 98031
Forward
Flow Research was under contract with the Environmental Protection
Agency (Contract No. 68-02-0800) to investigate plume dispersion in
stably stratified flows over complex terrain. Two tasks are involved:
Task I - Laboratory simulation of plume dispersion in stably
stratified flow over a complex terrain, which is covered in Flow
Research Report No. 29.
Task II - A feasibility study of numerical simulation of plume
dispersion in stably stratified flows over complex terrain,
which is covered in Flow Research Note No. 40 and Flow Research
Report No. 30.
These three reports, each of which can be separated as an
individual report, are bound together under the title Laboratory and
Numerical Simulation of Plume Dispersion in Stably Stratified Flows
over Complex Terrain. The order of the binding is task oriented.
-------
TASK I
LABORATORY SIMULATION OF PLUME DISPERSION IN
STABLY STRATIFIED FLOWS OVER
A COMPLEX TERRAIN
BY
JUNG-TAI LIN
HSIEN-TA Liu
YiH-Ho PAD
-------
Laboratory Simulation of Plume Dispersion in Stably
Stratified Flows over a Complex Terrain
by
Jung-Tai Lin, Hsien-Ta Liu, and Yih-Ho Pao
Flow Research, Inc., Kent, Washington 98031
Abstract
Laboratory simulations were conducted in a stably stratified towing
tank to investigate the effects of stability and terrain on plume
dispersion under extreme atmospheric conditions. A three dimensional
model was towed through the tank to simulate flow over complex
terrain. The terrain model contained the essential features of a complex
terrain, namely, mountains, mountain ridges and valleys. A stack was
located in the valley upstream of the mountain ridge between two mountain
peaks. The simulation scaling factor was 1:2500. The stack Froude number
was 3; the stack Reynolds number was 530t and the internal Froude number
F/ based on the mountain ridge height ranged from 1 to 3.
The characteristics of plume dispersion were investigated
using shadowgraph and dye visualization. With reference to plume
dispersion over a flat plate, the terrain effects revealed from
visualization results are:
1) The plume rise and spread increase as a result of the reduction
of the local flow speed— the upstream blocking phenomenon;
2) The lateral spread near the upstream side of the ridge is
enhanced significantly by the lateral flow around the terrain
model;
3) Accumulation of pollutants occurs in a semi-stagnant region just
upstream and below the ridge;
4) Downstream of the ridge, pollutants carried around and over the
terrain model remain aloft for small F/, but are carried downward
along the lee surface as F, increases.
The turbulent velocity and density in the plume were measured with hot-
film and conductivity probes. The normalized turbulent intensities and flatness
factors of velocity and density fluctuations are approximately self-similar.
-------
ii
The mean and root-mean-square pollutant concentrations, using temperature
as a tracer, were measured on the upstream slope of the terrain model.
Highest concentrations occur at stations on the ridge.
-------
iii
Table of Contents
Page
Abstract i
Table of Contents iii
Nomenclature iv
1. Introduction 1
2. Principles of Laboratory Simulation 3
3. Experimental Facilities and Methods 6
3.1 The Stratified Towing Tank 6
3.2 The Terrain Model and Platform 6
3.3 The Effluent Injection Device 7
3.4 Flow Visualization Techniques 8
3.5 Hot-Film Anemometer, Conductivity Gauge, and Temperature Gauge 8
4. Flow Visualizations 10
4.1 Plume Dispersion in Stratified and Nonstratified Flows over a
Flat Plate 10
4.2 Plume Dispersion in a Stratified Flow over the Terrain Model 12
5. Velocity, Density and Pollutant Concentration Measurements 15
6. Conclusions 17
7. References 19
Table 1 21
Table 2 22
Figures 23
-------
IV
NOMENCLATURE
Dimensions of each term are given in parenthesis: I- = length, t= time,
T = temperature, m = mass.
a_^, i = 0, ..., 5 = Coefficients defined in Eq. (3.1)
D = Stack Diameter (£)
F - Buoyancy Flux of Plume (/£3), g -^ D S
Si
w
o
F = Stack Froude number (Dimensicmless),
F = Internal Froude number based on the stack height (Dimensionless)
V*
=r ——
NH
F, = Internal Froude number based on the ridge height (Dimensionless)
"
U,
_2
g = Gravitational Acceleration (Lt )
h(x,y) = Terrain height (£)
h = Mountain peak height (£)
H = Stack height (£)
K = Ratio of effluent velocity to free stream . velocity (Dimensionless)
W
N = Brunt-Vaisala frequency (cycle t ) , TJ-T \M?~ ~r~
W D a
G
R^ = Stack Reynolds number (Dimensionless) ,
R, = Reynolds number based on the ridge height (Dimensionless)
2-2 2
s = Stability parameter (cycle t ), (2irN)
t = Time (£)
U = Instantaneous velocity component in x direction (It )
-------
U = Mean velocity component in x direction (JUL )
U = Mean velocity component in x direction at z = H (jit )
3.
U = Free stream velocity (tt )
OO
u = Fluctuating velocity component in x direction, (JUt )
W = Stack effluent velocity (t£~ )
s
x,y,z = Cartesian coordinate axes in longitudinal, lateral, and vertical
directions, respectively (£)
x ,y = x and y coordinates of the saddle point on the mountain ridge (£)
Z = Plume rise determined from visualization results (based on the stack
exit level) (£)
Z ,Z = Plume rise defined as the plume trajectories at which V u^ = vu^
and v cr = vcr > respectively (£)
UlclX
a = Longitudinal entrainment coefficient (Dimensionless)
g = Cross-flow entrainment coefficient (Dimensionless)
AY»AZ - Lateral and vertical spreads of plume determined from visualization
results (I)
Az , AZ - Vertical spreads defined as the width between the plume boundaries
at which
/ u2 12 and v®^ = a
' resPectively •
T = Adiabatic lapse rate (T£~ )
T = Lapse Rate (T£~ )
k = Mountain ridge height (£)
_3
p = Density (mt )
— —3
p « Mean density of the ambient fluid (m£ )
_3
p = Density of ambient fluid at stack exit level (mi )
3.
_3
p = Stack effluent density (mi )
S
_3
Ap = Ps ~ Pa (ml )
£ t
v = Kinematic viscosity (£ t )
_3
E = Reduced density (mt ), (p -1) x 1000
L = Reduced density of the effluent (m£~3), (p - 1) x 1000
s3 S
-------
vi
O = Fluctuating reduced density
0 = Temperature (T)
0 = Temperature of ambient fluid at stack exit level (T)
3.
0 = Effluent temperature (T)
s
9 = Fluctuating temperature (T)
-------
Laboratory Simulation of Plume Dispersion in Stably
Stratified Flows over a Complex Terrain
by
Jung-Tai Lin, Hsien-Ta Liu, and Yih-Ho Pao
Flow Research, Inc., Kent, Wash. 98031
1. Introduction
Serious air-pollution episodes usually occur in the areas where the
terrain features are dominantly irregular and complex, and during the
periods when the atmosphere is stably stratified. In a stably stratified
flow over a complex terrain, individual flow phenomena, namely, blocking,
velocity concentration, flow separation, lee waves, drainage, and heated-
island effects, influence the transport of pollutants in their own
particular ways. Therefore, the combined effects of the flow phenomena
on air-pollution problems are expectedly complicated. The present study
attempts to understand through laboratory experiments terrain effects,
not including drainage and heated-island effects, on the pollutant
transport in stably stratified flows.
For simulating the plume dispersion in a neutral atmosphere over a
complex terrain, Hino (1968) conducted experiments in a wind tunnel.
Veenhuizen et al. (1973) conducted similar experiments in a towing tank.
Orgill et al. (1970) simulated the aispersion of pollutants in stably
stratified flows over a complex terrain in a wind tunnel by cooling with
dry ice the flow near the terrain surface. In the present investigation,
the dispersion of a bent-over plume in stably stratified flows (with 1Lnear
density profiles in a towing tank) over an idealized three-dimensional
terrain model were studied. Both flow visualizations and probe measure-
ments were conducted in an attempt to gain further insight into the problem.
In this report, the principles of laboratory simulation are investi-
gated in Section 2 to determine the basic simulation parameters for the
plume dispersion in a stably stratified flow over a complex terrain. The
experimental facilities are described in Section 3. Flow visualization
results of plume dispersion over a flat plate and over a 3-D terrain are
-------
presented in Section 4; the experiments with the flat plate were conducted
to provide a reference for investigating the terrain effects. The results
of velocity and density measurements in the plume and pollutant concentration
measurements on the terrain surface are presented in Section 5. Conclusions
are given in Section 6.- In addition, flow visualizations are presented
in Flow Research Movie No. 6 as a part of the report.
-------
o
2. Principles of Laboratory Simulation
Geometric, kinematic, and dynamic similarities must be satisfied in the
simulation of plume dispersion in a stratified flow over a complex terrain.
The equation for plume rise in a stably stratified atmosphere, and
the equation for the wind field are, respectively
tf
s
and
Z = f(x,y,D,H,h(x,y),W ,a,a,U ,U,Y-T,g,v> (2.1)
o £> d OO
U = g(x,y,z,D,H,h(x,y),Y-r,U ,g,v) (2.2)
where Z is the plume trajectory (the centroid of the plume cross-section) above
the stack, x,y are the respective longitudinal and lateral coordinates referenced
to the stack base, D is the stack diameter, H is the stack height , h(x,y) is the
terrain profile, W is the effluent speed, 0 is the temperature of the effluent,
s s
Q is the ambient fluid temperature, U is the free stream wind speed,Y is the laps
SL
rate of atmosphere, F is the adiabatic lapse rate, g is the gravitational accelera-
tion, and V is the kinematic viscosity of the air. Figure 1 shows a definition
sketch of the most relevent variables.
Based on dimensional arguments, the following nondimensional equations
can be derived from Eq. (2.1) and Eq. (2.2),
Z/D = f(x/D,y/D,H/D,h(x/D,y/D)/D,U/Uco,K,FD,Fh,R1)) (2.3)
u/uoo = g(x/D,y/D,z/D,H/D,h(x/D,y/D)/D,F^,P^) , (2.4)
where
W
K = TT^- is the velocity ratio,
00
W
F = —/ ~. n. — ^s the stack Froude number,
* 6a
U
OO
F^ = — 1 > p \ "~ ls tne internal Froude number,
" "" / _r^ " *"" /
W D
R = —— is the stack Reynolds number,
u h
CO
R = —— is the terrain Reynolds number;
and ^ is a characteristic height of the terrain.
-------
Geometric similarity is satisfied when Z/D, H/D, and h(x/D,y/D)/D in
the laboratory have the same values as those in the field. When K and U/U^ in
the laboratory are equal to those in the field, kinematic similarity is
satisfied. Plume dispersion is dynamically simulated when F , F, , R , and
R in the laboratory are eflual to those in the field.
In the laboratory simulation, the fulfillment of the dynamic similarity
of F and F, is usually not a problem. The dynamic similarity of the stack
Reynolds number R cannot be fulfilled as the Reynolds number in the laboratory
is usually about 3-4 orders of magnitude smaller. However, experience (Hoult
et al. 1972) has shown that the Reynolds number effect is not important as long
as the plume is turbulent in a nonstratif ied fluid . The Reynolds number effect on plume
dispersion in the stably stratified fluid will be investigated in Section 4.1.
In the dimensional analysis expressed in Eqs. (2.3) and (2.4), the
turbulent structure in the atmospheric boundary layer is not considered.
In the case of a neutral atmosphere, the boundary layer may have a thickness
of about 300m which is higher than the average stack height of about 150m,
and thus, in the laboratory simulation of plume dispersion, turbulent
characteristics in the atmospheric boundary layer must be simulated properly.
In the case of a stably stratified atmosphere, the boundary layer is thin
and the atmosphere is basically not turbulent (Slade, 1969). Hence, in the
simulation of plume dispersion in a stably stratified flow, the boundary
layer and the terrain Reynolds number need not be critically simulated as
long as the stack height is greater than the boundary layer thickness and
the plume dispersion occurs outside the boundary layer.
The laboratory simulation of plume dispersion was conducted in a
stratified towing tank in which a wide range of K, F , and F, can be
achieved to meet the various field conditions. A salt solution was used
for both the effluent and the ambient stratified fluid. The stack Froude
number F and the internal Froude number F, are then redefined as
FD =
PS
and Uoo
(2.6)
-------
where Ap = p - p , p is the density of effluent at the stack exit, p is
S SL S 3.
the density of the ambient fluid at the level of stack exit and N is the Brunt-
Vaisala frequency (cycle/sec) defined as N = T~ /- j , where ~
4r \/ P dz dz
is the mean density gradient of the ambient fluid.
A modeling example is given in the following by considering a prototype
stack of 127 m height and 7.94 m in diameter. The effluent temperature Q
is 160 C; the effluent velocity W Ls 20.4 m/sec; and the ambient temperature
S
at the level of the stack exit Q i.s 5 C. The stack Froude number is
a
therefore F = 3.09. Further consider a stable atmosphere having a lapse
rate y = .025 c/m, a wind speed U - 6 m/sec, and a characteristic height
of the terrain h = 407 m. The internal Froude number is therefore F< =
2.64; and the velocity ratio is K = 3.4. If a scaling factor 1:2500 is selected
for the laboratory simulation, the model stack height and diameter will
be 5.0 cm and 0.32 cm, respectively; and the characteristic height of the
model terrain is 16.3 cm. To satisfy the simulation of the internal Froude
number F, , 11^= 6 cm/sec and N = .14 cycle/sec are used. To meet the
requirements of the velocity ratio K and the stack Froude number F , W = 20.4
3 3 D s
cm/sec, p = 1.175 gm/cm , and p = 1.015 gm/cm are used.
S 3.
-------
3. Experimental Facilities and Methods
3.1 The Stratified Towing Tank
The plume dispersion experiments were performed in a stratified towing
tank, 18.3m long, 1.2m wide, and ,91m deep. See Flow Research Report No. 4
(Pao et al. 1971) for a detailed description of the stratified towing
tank. In the following, some special features of the tank are described
briefly.
The towing tank was designed and constructed for investigating
turbulent wake and internal wave phenomena in stably stratified fluids.
1) It has transparent side, bottom, and end walls to allow flow visualiza-
tion from all directions. 2) It has a filling system to provide a stratified
fluid at a predetermined stratification within four hours. Salt (NaCl)
solution is the working fluid for making a stratified fluid. The temperature
of the fluid is controlled within 0.1°C. 3) It has a smooth oil-lubricated
carriage to allow for accurate measurements of turbulent velocity,
temperature, and salinity fluctuations, and the carriage can be used- to tow
the terrain model in the tank. 4) It has a mini-computer system for direct
on-line data acquisition and analysis. The mini-computer system can handle
a maximum data rate of 44,000 samples per second at 13 bits accuracy. The
system has a memory core of 20,000 words and a disc storage of 2.5 million
words at 16 bit accuracy. In conjunction with this mini-computer system,
a series of computer programs has been developed to obtain statistical moments
up to fourth-order.
3.2 The Terrain Model and Platform
o
A terrain model was made of polyurethane foam (9.6 Kg/m density) mounted
on a 3.66 m x 1.14 m x .64 m acrylic plastic sheet reinforced with aluminum
angles and bars on the back side. The mountain terrain section was inter-
changeable with a flat surface. The surface of the plastic sheet was covered
with a sheet of .64 cm polyurethane foam. The foam surface was studded
with artificial roughness elements, .32 cm in diameter and having an average
2
height of .32 cm, with a density of about 1 element/cm . The entire model
was painted flat white to improve photographic quality.
-------
The vertical profile h(x,y) of the terrain model is described by the
sum of three Gaussian functions
22 22
h(x,y) = ao{exp[a1(x-xQ) + a2(y-yQ) ].+ exp[a1(x-xQ) + a2(y+yQ) ]}
+ a3{exp[a4(x-xo)2 + a5y ]} , (3.1)
where the coordinates are referenced to the stack base. The following
coefficients, aQ = 17.80, a.^ = -.0008513, a = -.01197, a = 16.0,
a, = -.01171, a = -.002314, x = 61, and y = 18.82 (length dimension
in cm) were selected. The terrain model had the essential features
of a complex terrain, namely, mountains, mountain ridges, and valleys.
Cross-sectional views of the terrain model are plotted in Fig. 2; and
contour map is plotted in Fig. 3. The characteristic length scales of
the terrain and the stack are summarized in Table 1.
In the stratified towing tank, the terrain model and the stack were
actually mounted upside down under the oil-lubricated carriage, and the
effluent, heavier than the ambient fluid at the stack exit, was injected
downwards from the stack exit into a stably or a neutrally stratified fluid.
The terrain model was towed at a speed specified for a particular experiment.
This arrangement was used to simulate the hot effluent rising upwards from
a stack into a stable ora neutral atmosphere.
3.3 The Effluent Injection Device
Concentrated salt (NaCl) solution was used as the effluent fluid.
The effluent device consisted of a pressurized storage chamber for the
salt solution, a portable tank of compressed air, a needle valve for
regulating the effluent discharge, and a flow meter for measuring the
flow rate. A dye container filled with a blue food dye was inserted
between the flow meter outlet and the stack. The stack was made of
stainless steel tubing with an inside diameter of .318 cm. The effluent
flow was tripped with a small ring (.16 cm I.D.) installed inside the
tubing about 1 cm from the stack exit, and the plume was fully turbulent
at the stack exit. Figure 4a shows the effluent injection device.
To measure pollutant concentrations on the terrain surface, temperature
was used as a tracer. The plume was heated by a portable boiler, which was
made of a stainless steel filter case (Cat. No. BRX10-3/4B, Commercial
-------
8
Filter Corporation) wrapped with a 400 watts heating tape (Cat. No. 33745-120,
Van Waters and Rogers, Inc.)> and was connected to the effluent feeding line
between the dye container and the stack. By adjusting a voltage regulator
(a Variac Transformer) connected to the heating tape, the plume temperature
could be controlled. The temperature of effluent at the stack exit was
measured with a 10 gauge thermo-couple probe (SCXSS-010G-36" T/C Omega
Engineering, Inc.) inserted inside the stack. The plume temperature was
kept at 77°C + 2°C for ground concentration measurements.
3.4 Flow Visualization Techniques
Shadowgraph and dye methods were used to visualize the plume dispersion
in stratified and nonstratified fluids. Since the shadowgraph intensity is
directly proportional to the second derivatives of the refraction index,
the small structure of turbulent motions in a plume of a salt solution and
in the turbulent boundary layer of the stratified flow over the terrain surface
can be visualized from shadowgraph pictures. For the present experiments, a
parallel light beam of diameter 30.5 cm was obtained by projecting a point
light source at the focal point of a parabolic mirror of 30.5 cm in diameter
and 2.5m in focal length. A Xenon-Mercury D.C. light source (1000 watts)
equipped with condensing lens and pin-holes was used as the point light
source. When a parallel light beam passes through the turbulent stratified
plume, a shadowgraph is produced.
Dye visualization methods were used to visualize the large scale
plume motions and to trace the plume • Movie and still cameras
were used to record visualization data during the experiments.
3.5 Hot-Film Anemometer, Conductivity Gauge, and Temperature Gauge
The mean and fluctuating velocities in the plume were measured with a
10-channel constant temperature anemometer (Model 1051-10, Thermo-Systems, Inc.)
using conical and cylindrical hot-film probes (TSI Model 1230S and 1290AK)
as sensors. The anemometer has a frequency response up to 1 KHz. Figure
5 a shows a typical velocity calibration curve plotted as the output
voltage vs velocity. The sensitivity is about 2 volts per cm/sec.
-------
Flow Research has developed a fast response temperature gauge which
measures resistance changes of a metallic film corresponding to variations
of the fluid temperature.The gauge employs the sensor as one arm of a balanced
bridge driven by a 10 KHz oscillator which supplies a low current (several
MA) to the bridge. The bridge is initially balanced by means of an adjustable
resistor. As the resistance of the sensor varies, the signal of imbalance
modulates the 10 KHz carrier frequency, and is then amplified and synchronously
detected.
Cylindrical hot-film probes (TSl Model 1290AK) were used as the
sensors which had a temperature coefficient of .0026 f2/ C and a natural
frequency (thermal) of 850 Hz. A calibration circuit which changes the
current through one side of the bridge by .26% (equivalent to 1 C temperature
change) was built in the gauge. The output of the bridge during calibration
was 15Hz square waves whose peak-to-peak amplitude corresponded to a 1C
temperature change. The linearity of the gauge was within 0.1%, The
temperature gauge was calibrated against a mercury thermometer with .1 C
increment in the water bath. A typical temperature calibration curve is
shown in Figure 5 b ; the sensitivity is about 1.5 V/ C. The equivalent
noise level in temperature measurements was .004 C.
Ten conductivity gauges (Model 1110) developed by Flow Research, Inc.
were used to measure density profiles in the plume. This gauge has a
frequency response from DC to 1 KHz. The sensor was a single electrode
conductivity probe which had a .025 mm stainless steel tip platinized with
platinum black solution. The sensitivity of the conductivity gauge is about
1 volt per 0.025 gm/cm . Figure 6 a shows a typical calibration curve of
the conductivity probe where E is a reduced fluid density defined as
£ = (p- 1) x 1000, where p is the specific gravity of the fluid.
Figure 6b is a typical density profile which has a fairly constant gradient
over about 60 cm. The experiments were conducted at a depth of about 15 cm
below the water surface.
-------
10
4. Flow Visualizations
4.1 Plume Dispersion in Stratified and Nonstratified Flows over a
Flat Plate
The plume spreads (R = 530) in stagnant nonstratified (F = 2.84)
and stratified fluids (F = 2.92) were observed with shadowgraph pictures
shown in Figs. 7a and 7b. The plume in a nonstratified fluid spreads
continuously as it rises. The plume in a stratified fluid (N = 0.14 cycle/sec)
spreads in the initial stage as the plume in a nonstratified fluid does, but
in the stratified fluid, the plume rise is limited. The shadowgraph shows
that the plume spreads outward in horizontal directions after it reaches
its terminal rise. The plume is turbulent except in the outer region of the
spreading disc.
The entrainment coefficient a was measured using R = 2az, where R is the
radius of visible plume boundary, and has a value ex = 0.08. The entrainment
coefficient observed by the others, are, a = 0.093 by Morton, et al. (1956),
a = 0.075 by Briggs (1968), and a =- 0.08 by Ricou et al. (1961).
The dispersion of bent-over plumes (K = 1.7, R^ = 530) in nonstratified
(F = 2.84) and stratified (F = 2.92) flows is shown in Figs. 8a and 8b.
The corresponding internal Froude numbers F = U^/NH are °° and 1.69,
respectively. The bent-over plume in a nonstratified flow is turbulent and
spreads continuously as it rises. The shadowgraph picture shows that the
upper boundary of the plume is rougher than the lower boundary because in-
side the plume the mean density stratification is unstable in the upper half
of the plume and is stable in the lower half .
Figure 8b shows that both the plume and the boundary layer, which has
a thickness about 2.5 cm, in the stratified case are turbulent. In the near
downstream of the stack, the plume in the stratified flow spreads and rises.
But further downstream, the plume reaches a terminal height and the vertical
spread is limited. The trajectories of plumes, which are defined as the
mid-point of the upper and lower plume edges, are determined by averaging
five ensembles. They are plotted as ZuJ^F vs xU^/F in Fig. 9, where F
2
is the buoyancy flux of the plume defined by F = W (p -p )gD /4p .
o o a. a-
The plume in a nonstratified flow rises continuously and its trajectory
is correlated by the equation
3 / 3 \ 2/3
ZU J J '
°r ^ = 1.6 KFD"'J ( r r/J- (4.2)
-------
11
The coefficient in the right hand side of Eq. (4.1) was observed to
be 1.8 by Briggs (1969), 1.44 by Hewett, at. al. (1970), and 1.6 by
Turner (1973). In the stratified flows, the vertical rise of plumes
is limited by the stable stratification, and the plume oscillates with
x after its maximum rise is reached.
Trajectories of plumes in stratified flows are re-plotted as
I/"} -1 /2 '
Z/(F/Uoos) vs x/(UOTs ) in Fig. 10, where s is the stability
2
parameter defined by s = (2irN) . Other laboratory data (Hewett et al.
1970) and field data from TVA plants (see Briggs 1969) are also
plotted in Fig. 10. For high F , the maximum plume rise occurs at
-1/2
about x/(Uoos ) = 4 or Nx/U^ = 0.65, and the terminal plume trajectory
has a height about 2.2. When F decreases, the plume rise increases
n
as possibly the result of the reduction of cross-flow entrainment with
increasing stability. This could explain partially for the scatter
of the field data plotted in Fig. 10.
The laboratory and field data of terminal plume rise in stratified
1/S — 3/8
flows are plotted as Z vs 5 F S in Fig. 11. The terminal plume
rise is defined as the final height of the plume top. Our experimental
data of the plumes in stagnant stratified fluids (FT1 = 0) and in stratified
n
flows (F = 3.11 and 6.22) are in good agreement with the data obtained
n
by the others (Briggs 1969). The agreement shown in Figs. 10 and 11
would then support the arguments (Section 2) that the stack Reynolds
number need not be simulated critically in a stably stratified flow
as long as the plume is turbulent at the stack exit.
The vertical and lateral spreads of plumes in a nonstratified flow
3 3
over a flat plate are plotted as AZU^ /F and AYU^ /F in Figs. 12 and 13,
respectively, where the plume spreads are determined by averaging 5
ensembles. The vertical and lateral spreads are correlated by
-------
12
(4.4)
and AYU 3 /xU 3 \ 2/3
00 I 00 I
—T~ = 1-5 ^-jT-j (4.5)
xU3
00
in the range —g—> 6. The lateral spread is greater than the vertical
r
spread, and the aspect ratio of the plume cross section is about 1.2.
A top view of a dyed plume in a nonstratified flow is shown in Fig. 14 a
The plume appears to split into two parts, which indicates the formation
of the twin vortex. The twin vortex has been observed in the field under
some conditions (Scorer, 1968).
Figures 12 and 13 show the vertical and lateral spreads
of plumes in stratified flows. In the region where the buoyancy dominates
the plume rise, Eq. (4.5) provides the limiting spreads of the plume in
stratified flows. The vertical spread is limited by the stable stratification
and appears to oscillate with x after it reaches its maximum value. A top view
of a dyed plume in a stratified flow (FTI = 6.22) is shown in Fig. 14 b .
ri
The twin-vortex structure is not observed because the vertical motion of
the twin vortex is suppressed in a stratified flow.
4.2 Plume Dispersion in a Stratified Flow over the Terrain Model
A side view and a top view of a dyed plume at F^ = 3.IS, !1 = 5.3, and
R_^ = 530 in a stratified flow (F, = 0.97) over the terrain model are
presented in Fig. 15. The flow direction was from right to left at
U^ = 3 cm/sec. The stack was located at x/D = 192 upstream of the ridge
between the two mountain peaks. The plume appears to be fully turbulent
and tends to relaminarize as it levels off. Before the plume levels
off, it tends to overshoot. The residue momentum of the plume which
dissipates as the entrainment process progresses, initiates the
overshooting. Further downstream, the plume remains aloft in both
upstream and downstream sides of the terrain. The top view of the dyed
plume shows the spread of the plume in the lateral direction. The pollutant
accumulation is seen around the I-shaped bright area (Fig. 15 b ), which
corresponds to the two mountain peaks and the ridge between them (See Fig. 3).
-------
13
The effects of the flow speed on the plume dispersion over the terrain
were examined by increasing the flow speed from U^ = 3 cm/sec to 6 cm/sec
but leaving the other flow conditions unchanged (Figs. 16 a and 16 b ).
When the wind speed increases, the plume bends more toward the downwind
side of the stack, and the overshoot of the plume becomes less pronounced.
Further, in the lee side of the mountain ridge, the plume is carried toward
the ground surface by the downslope flow resulting from the increase in flow
speed.
Two additional cases with weaker stratifications (N = .14 cps) than
those described above are shown in Figs. 17 and 18. It can be readily
observed that the terminal rise increases with decreasing stability while
the gross features of the plumes remain unchanged.
The terrain effects are investigated by comparing the vertical and
lateral boundaries of a plume over the terrain with those of a plume over
a flat plate (Fig. 19). With the presence of an obstacle, the stably
stratified flow upstream is blocked. Velocity profiles presented in
Section 5 indicate that the upstream blocking extends beyond the distance
where the stack is located. The upstream blocking which reduces the mean
advection speed U at the stack exit effectively increases the local speed
cl
ratio W /U . As a result, the plume rise increases with increasing
S 3.
blocking.
In a stably stratified flow, a semi-stagnant or plugged region (Turner
1973) usually develops just upstream and below the obstacle. The
accumulation of pollutants observed in Fig. 16 confirms the existence of
such a region.
When the obstacle is three-dimensional as is the terrain model used
in the experiments, the stably stratified fluid tends to flow around
rather than over the obstacle because the vertical motion is inhibited.
The induced lateral current carries a portion of the pollutants around
the mountain peaks and therefore effectively enhances the lateral spread
of the plume as observed in Fig. 19.
Plots of the plume trajectories and the plume spreads AZ/D and
AY/D vs Nx/U^ are presented in Figs. 20, 21, and 22, respectively.
The rise and the spreads for plumes over the terrain model
are always greater than those over the flat plate provided that the induced
lateral current is relatively small. Under such a circumstance, the up-
stream blocking contributes mainly to the increase of the rise and the
-------
14
spreads. As the plume approaches the mountain ridge, the induced lateral
current becomes pronounced. It increases the lateral spread considerably
and consequently reduces the vertical spread as can be seen in Fig. 22.
Both the maximum rise and the maximum vertical spread tend to occur about
Nx/U^ = .65 whether the flow is over the terrain model or over the flat
plate. Apparently, these two features only depend on the buoyancy of
the plume and the stability of the ambient fluid.
The dispersion of bent-over plumes from a stack located downstream
of the ridge was also investigated. Figure 23 a shows a dyed plume
released from a stack F = 2.09 and F^ = 1.32. The plume rises
straight upwards and spreads horizontally in both upstream and downstream
directions, as the maximum rise is reached. However, at a higher internal
Froude number F, = 2.64 (Fig. 23b), the plume is bent over toward the
downstream direction and is carried downstream by the downslope flow.
The formation of lee waves traced by the dispersing plume and the turbulent
layer near the ground surface can be clearly observed in Fig. 23b.
-------
15
5. Velocity, Density and Pollutant Concentration Measurements
The turbulent plume in a stratified flow was investigated by measuring
the velocity and density profiles in the plume from a stack which was
located in the valley and at 192 D upstream from the ridge between the
mountain peaks. The stack exit was at an elevation x/D = 16, and the
mountain ridge had a height k/D = 510. The free stream velocity U^ was
6 cm/sec, the Brunt-Vaisala frequency was N = .14 cycle/sec, and the
internal Froude number was .therefore F = U^/NL = 2.64. The stack Froude
number F was 3.09, the velocity ratio K was 3.4, and the stack Reynolds
number was R^ = 530.
The conductivity probes were aligned vertically at the plume axis
and the hot-film probes were 1 cm away from the plume axis (see Fig. 4b) .
The velocity and density measurements were conducted at three downstream
stations from the stack, x/D = 32, 64, and 96. The plume trajectories
indicated by the position of the maximum velocity intensity \/u max
l~2~
by the maxium density fluctuation ^ max are plotted in Fig. 24. They
are in fair agreement with the plume trajectory determined from dye
visualizations .
The mean velocity profiles U/U^ in the upstream of the terrain ridge
are plotted vs z/D for x/D = 32. 64, 96. The significant variations in
the velocity in the vertical direction is probably not caused by the
boundary layer but rather by the terrain because the turbulent boundary
layer has a thickness only of about 10 D. Further, the velocity varies
significantly with the distance from the mountain ridge and the effects
of terrain on a stably stratified flow (upstream blocking) are therefore
clearly demonstrated.
The maximum values of the longitudinal velocity intensity vu' „„,
Tiltl A. -
at 3 downstream positions are plotted vs x/D in Fig. 26. The decay of
velocity intensity is approximated by
'max/Ws = 2'8 (X/D^"1' t5-1)
f=f
The maximum values of the density intensity^a r^x/^ are also plotted
us x/D in Fig. 26 where £ = (ps -1) x 1000 is the reduced density of
s
the effluent at the stack exit and a is the turbulent fluctuation of
the reduced density. The maximum density intensity decays in the range
30 < x/D < 60 but incresase instead at x/D = 96.
-------
16
The normalized velocity intensities \/u2/v u2 are plotted vs
max
2(z-Z )/AZ in Fig. 27 where Z as the vertical coordinate of the maximum
velocity intensity, Az is the width between the plume boundaries at which the
velocity intensity is half of its maximum. The normalized profiles are not
symmetric with the plume trajectory but approximately self-similar. Also plotted in
Fig. 27 are the normalized density intensities va2/N/a2 Vs 2(z-Z )/Az
where Z^ is the vertical coordinate of the maximum intensity of density
fluctuations, and AZ is the width between the boundaries at which the intensity
of density fluctuations is half of its maximum.
/ 99
The flatness factor of the velocity fluctuations u /(u ) is plotted
vs 2(z-Z )AZ in Fig. 28a. Inside the plume the flatness factor has an
average value of about 2.8, and near the plume boundary, the flatness factor
has a_sharp variation. The flatness factor of the density fluctuations
a /(a ) is plotted vs 2(z-Z )/Az in Fig. 28b and it has an average value
of about 3 inside the plume but has a sharp variaton near the plume boundary.
It is noted here that it has been observed from the velocity and density
measurements in a turbulent wake in a stratified fluid (Pao et al. 1973, and
Lin et al., 1973) that the flatness factor remains fairly constant inside
the wake but has a sharp variation near the wake boundary.
The pollutant concentrations on the ground level of the terrain were
measured using temperature as a tracer. The locations of the five ground
stations on the upwind slope of the terrain are shown in Fig. 3 and are
summarized in Table 2 where x,y, and z are the coordinates referenced to
the stack base. The stack Froude number was F = 3.16, the stack
Reynolds number was R = 530, and the plume temperature at the stack exit
was 0 = 77 + 2 C. The ambient fluid was stratified with a Brunt-Vaisala
s
frequency N = 0.14 cycle/sec, and the ambient temperature was Q = 19 C.
3.
The towing speeds were at U =3 cm/sec and 6 cm/sec, and the corresponding
oo
internal Froude numbers were F = 1.32 and 2.64,respectively. The mean and
rms pollutant concentrations determined from the mean and rms temperature
measurements are tabulated in Table 2b. For F, =1.32, the mean concentration
n
at the station 1 is the highest and the rms concentration at all stations
is not significant. However, as F^ increases to 2.64, the mean
concentration at station 1 is three times reduced but the
rms concentration increases by a factor of about 9. For the terrain investigated
here, the highest pollutant concentration (mean and rms) always occur at the
mountain ridge.
-------
17
6. Conclusions
i. The fair agreement among the present results of plumes over a flat
plate in the towing tank and other laboratory and field data (Figs. 10 and 11)
suggests that the stack Reynolds number IL need not be critically simulated
for plume dispersion in a stable atmosphere. It is required, however, that
the plume should be turbulent at its exit for a correct simulation in both
the near and the far fields.
ii. With the presence of the three-dimensional terrain model, the stably
stratified flow (F^ = 1 to 3) upstream of the mountain ridge is blocked. The
velocity data (Fig. 25) verify the existence of the upstream blocking. As a
result, the plume rise and spread increase because the local speed ratio
W /U increases. In one case (Fig. 16), a significant amount of pollutants
£> 3. *•
is trapped in the semi-stagnant region just upstream of and below the mountain
ridge.
iii. The stably stratified fluid tends to flow around rather than over
the terrain model because the vertical motion is inhibited. Near the mountain
ridge, the induced lateral current is pronounced. The lateral spread increases
considerably and consequently the vertical spread reduces.
iv. At a small Froude number (F/ = 1.32) the pollutants carried around
and over the terrain model remain aloft. As the Froude number increases
(Fi = 2.64), pollutants downstream of the terrain model are carried closer
to the ground surface by the downslope flow.
v. When the stack is located in the lee side of the terrain model, the
trajectory or the plume is strongly governed by the internal Froude number.
For small F^ (approx. 1.32), the plume, spreads into both the upstream and
downstream sides of the stack because the flow field adjacent to the downstream
side of the mountain ridge is nearly stagnant. However, at a higher Froude
number (F = 2.64), the plume is advected downstream by the downslope flow.
h
vi. The plume trajectories determined as the loci of the maximum rms
velocity and rms density in the plume are in fair agreement with that determined
from the visualization results.
vii. The rms velocity and rms density profiles normalized by the respective
maximum values are not symmetric with che plume trajectory. These profiles
-------
18
(Fig. 27) plotted against the height normalized by the respective half
radii, appear to be approximately self-similar. The normalized flatness
factor of the velocity and density fluctuations in the plume (Fig. 28;
are also approximately self-similar. The flatness factors remain almost
constant but have a sharp variation across the plume boundary .
viii. The pollutant concentrations in the upstream surface of the terrain
were measured by using temperature as the tracer. The mean and rms con-
centrations at stations in the mountain ridge are higher than
those at other stations.
-------
19
Preferences
Briggs, G. .A, (1968). "Penetration of Inversions by Plumes", Paper
presented at 48th Annual Meeting of the American Meteorology Society,
San Francisco.
Briggs, G. A. (1969). Plume Rise, U.S. Atomic Energy Commission
Critical Review Series, TID-25072.
Bringfelt, B. (1968). "Plume Rise; at Industrial Chimneys," Atmospheric
Environment, 2_, 575-598.
Fay, J. A., M. Escudier, and D. P. Hoult (1969). "A Correlation of Field
Observations of Plume Rise", Fluid Mech. Lab. Publication No. 69-4, Dept.
of Mech. Engineering, MIT.
Hewett, T. A., J. A. Fay, and D. P. Hoult (1970). "Laboratory Experiments
of Smokestack Plumes in a Stable Atmosphere", Fluid Mech. Lab., Dept.
of Mech. Engineering, MIT.
Hino, M. (1968). Computer Experiment on Smoke Diffusion over a Complicated
Topography. Atmospheric Environment, 2_, 541-558.
Hoult, D. P., and J. W. Weil (1972). "Turbulent Plume in a Laminar Cross
Flow", Atmospheric Environment, ^, 513-531.
Lin, J. T., and Y. H. Pao (1973). "Velocity and Density Measurements in the
Turbulent Wake of a Self-Propelled Slender Body in Stratified and Nonstratified
Fluids", Flow Research Report No. 34, Kent, Washington
Morton, B. P.., G. I. Taylor, and J. S. Turner (1956). "Turbulent Gravitational
Convection from Maintained and Instantaneous Sources". Proc. Roy. Soc.
London, A234, 1-23.
Orgill, M. M., J. E. Cermak, and L. 0. Grant (1970). "Laboratory Simulation
and Field Estimates of Atmospheric Transport-Dispersion over Mountainous
Terrain," Technical Report CER70-71MMO-HRC-LOG40, Colorado State University.
Pao, Y. H., J. T. Lin, R. L. Carlson, and L. P. C. Smithmeyer (1971). "The
Design and Construction of a Stratified Towing Tank with an Oil-Lubricated
Carriage". Flow Research Report No. 4, Kent, Washington.
Pao, Y. H., and J. T. Lin (1973). "Velocity and Density Measurements in the
Turbulent Wake of a Towed Slender Body in Stratified and Nonstratified
Fluids", Flow Research Report No. 12, Kent, Washington.
Ricou, J. P. and D. B. Spalding (1967). "Measurement of Entrainment by
Axialsymmetrical Turbulent Jets", J. Fluid Mech., 11, 21-32.
-------
20
Slade, D. H. (1969). "Low Turbulence Flow in the Planetary Boundary
Layer and Its Relation to Certain Air Pollution Problems", J. of
Applied Meteorology, 8^, 514-522.
Scorer, R. S. (1968). Air Pollution, Pergamon Press, London.
Turner, J. S. (1973). Buoyancy Effects in Fluids. Cambridge University
Press.
Veenhuizen, S. D., J. T. Lin, Y. H. Pao, D. W. Peecher and G. L. Hiatt.
(1973). Laboratory Simulation of Plumes From Kennecott Copper Smelters
in Garfield, Utah: Neutral Atmosphere. Flow Research Report No. 9,
Kent, Washington.
-------
TABLE 1. STACK AND TERRAIN LENGTH SCALES
Scale factor 1:2,500
Descriptions
Stack Diameter D
Stack Height H
Mountain Peaks Height h
Mountain Peak Separation
Mountain Ridge Height h.
Prototype
7 . 94 (m)
127.0 (m)
609.6 (m)
941.0 (m)
406.4 (m)
Model
.32 (cm)
5.08 (cm)
24.38 (cm)
37.64 (cm)
16.26 (cm)
-------
22
Table 2 SUMMARY OF GROUND SURFACE TEMPERATURE MEASUREMENTS
a. Coordinates of Ground Stations
Station
1
2
3
4
5
Station
1
2
3
1
2
3
4
5
x(cm)
61.0
56.3
53.8
61.0
56.3
b.
"co C
k * N£ "s*1
1.32 77
1.32 77
1.32 77
2.64 77
2.64 77
2.64 77
2.64 77
2.64 77
y(cm)
0
0
0
6.4
6.4
Experimental Data
6/\
^ff
c> oa( c) e< c> 0^
19.55 19.76 3.7xlO~3
19.55 19.58 5.2xlO~4
19.55 19.58 5.2xlO~4
20.65 20.73 1.4xlO~3
20.65 20.68 5.3xlO~4
20.65 20.70 8.9xlO~4
20.38
20.38 20.40 3.5xlO~4
z (cm)
16.3
12.7
8.9
16.3
12.7
V7
7.0xlO~5
-4
1.2x10
1.7xlO~5
6 . OxlO~4
1.6xlO~4
-4
1.1x10
4 . 9xlO~4
1.4xlO~4
-------
23
•o a* <-.
c 01 o
TH a. >~
li O^ IX
C
O
iJ
O
.*:
•r
C
a>
-------
24
>?
o
_
2
<
H
w
H
o
fl
+
o
CM
I
O
ro
I
01
T3
O
O
•CM
+
O
r-t
+
• O
o
. ,__(
1
x— v
B
a
s — '
o
X
1
X
w
u
S3
<
H
C/3
M
Q
>J
<
a
t— i
o
E3
H
[^--j
O
»
O
^J
•H
efl
J— 1
S-l
QJ
H
0)
s:
4-1
4-1
O
en
2
QJ
•H
>
i— 1
a
c
o
• H
4J
U
(!)
en
1
en
en
o
J_4
O
D
Ml
•H
(nio)z
-------
• Stack
X Stations tor Velocity and
Conductivity Measurements
Stations tor Surface
Temperature Measurements
x
i
x
w
u
C/D
M
o
<:
M
I
z
o
70
60
50
40
30
20
10
-10
-20
-30
-40
-50
-60
-70
-40 -30 -20 -10 0 10 20 30 40
LATERAL DISTANCE y(c-m)
Figure 3. A Contour Map of the Terrain Model
-------
26
To Stack
Flow Meter
Needle Valve
Pressure
a. Effluent Injection Device
Probes
b. Experimental Set-up for Velocity and
Conductivity Measurements
Figure 4. The Effluent Injection Device and Arrangement for
Probe Measurements
-------
2 r-
o
>
w
o
H
^J
o
-2
-4
1
L
6 8
VELOCITY (cm/s)
10
12
a. Velocity Calibration of a Hotfilm Probe with a Constant
Temperature Anemometer
£• o
o
>
w
-1
-2
22
23
24
TEMPERATURE ( C)
h. Temperature Calibration of a Thinfilm Probe with a Temperature
Gauge
Figure 5. Typical Calibrations of Velocity and Temperature Probes
-------
28
10
e
u
H
P-.
20
30
40
50
60
70
80
DENSITY (p-1) x 1000 (gm/cm )
60
40
20
T
T
T
a. Calibration Curve of a
Conductivity Probe
b. Density Profile
10 20 30 40
DENSITY (p-1) x 1000 (gm/cm3)
2.4
1.6
w
<
H
-.8
-•-1.6
50
Figure 6. Typical Calibration Curve of a Conductivity Probe and
the Density Profile in the Towing Tank
-------
29
a. Nonstratified Fluid (N = 0 Hz)
F - 2.84
b. Stably Stratified Fluid (N = .14 Hz)
18
16
14
12
10
8
6
4
2
0
0
O
EC
I
N
18
16
14
12
10
8
6
4
2
0
/-x
o
^/
EC
1
N
FD = 2.92
Figure 7. Shadowgraphs of Vertical Plumes in Stagnant Fluids
-------
30
(U10)
(mo) H - z
o
0)
CO
CO
Csl
II
fe
o
01
•H
4-J
CO
C
O
o
0)
CO
CM
CN
Q
cfl
QJ
•H
l-M
CO
C
O
CO
0)
3
l-(
PM
II
ffi
Pn
^5
O
rH
0)
•H
M-l
•H
4_J
cfl
J_l
4J
CO
>,
iH
.0
CO
4-)
CO
.
rQ
O
1
C
0)
pq
CO
M_j ^5
O O
CO PL,
p, T^
nj CD
M -H
&0 M-l
S -H
O •^-J
TD CO
cO J-i
r; 4J
CO CO
t
CO
-------
31
o
o
O
o
o
Ol
•r~1
u-4
•H
o
o
O
S-i
4-1
C/0
T3
P
CO
•a
•H
U!
c
o
CO
0;
B
Oi
> 0;
O 4-1
I ^
C P-.
CQ 4->
CO
O &-.
co co
0)
•H &-J
O >
4-1 O
o
ai co
co o
(-. ^H
H fe
c^
-------
32
en
•H
3
T)
0)
S-i
4-J
CO
C
V)
01
CD
C CO
4-1 Q
O
CU T3
•r-1 iH
Cfl 0)
J-i -H
E~^ CT |
,
>4-l O
O 4J
Cti
C O
O XI
en co
to ao
ex c
6 o
o 6
u co
o
3
00
-------
33
10000
1000
tsj
w
PH
C
H
100
10
.1
.1
D
Davies (Long Beach)
Vehrencamp Ambrosio, Romie
(Mojave Desert)
Crawford and Leonard (Ice Ring)
Morton, Taylor and Turner (Tank)
10
100
1000
10000
Figure 11. Laboratory and Field Data of Terminal Plume Rises in Calm and
Stable Fluids (See Briggs 1969).
-------
34
sa
fa
oo
ro
CN
CN
vO CO
OO iH
O
fa CN ro
o a <
o
o
T3
(U
o
-a
c
tfl
OJ
•H
C
•H
cn
0)
3
CO O
O
M O
01 iH
> fa
D
oc
•H
fa
ro
-------
o
o
oo
-------
36
1
0
1
10
1 1
20
1 1
30
1
40
1
50
x (cm)
a.
u
H
and F = 2.84
10
x (cm)
20
30
40
I
b. FH - 6.22 and FD - 2.92
50
I
Figure 14. The Lateral Spreads of Bent-over Plumes
in Nonstratifled and Stratified Flows
over a Flat Plate. K = 3.4 and R = 530.
-------
37
r- 40
-30
- 20
-10
- 0
o
a. Side View
U = 3cm/s
,--30
--20
.-10
- 0
- 10
20
30
b. Top View
Figure 15. Dye Visualizations of Plume Dispersion in a Stratified
Flow over the Terrain Model. F = 3.16, F, = .97,
K = 6.8 and R^ = 530. D h
-------
38
i- 40
-30
- 20
- 10
- 0
0
o
^x
ffi
i
N
a. Side View
U = 6cm/s
I--30
--20
- -10
- 0
- 10
- 20
L 30
b. Top View
Figure 16. Dye Visualizations of Plume Dispersion in a Stratified
Flow over the Terrain Model. F = 3.16, F, = 1.94, K = 3.4,
and R = 530 D "•
-------
39
'
/•— \
E
U
3J
1
N
40 -
30 -
20 -
10 -
0 -
a. Side View
U = 3cm/s
5
o
30
20-
10-
-10-
-20-
-30-
b. Top View
Figure 17. Dye Visualizations of Plume Dispersion in a Stratified Flow
over the Terrain Model. F = 3.09, F/ = 1.32, K = 6.8, and
RD = 530.
-------
40
e
o
S—X
ffi
I
N
B
o
30
20 -
10
0
30 -,
20 -
10 -
0 -
-10 -
-20 -
-30 -1
a. Side View
*«<,* <*•**' *.
U = 6cm/s
00
b. Top View
Figure 18. Dye Visualizations of Plume Dispersion in a Stratified Flow
over the Terrain Model. F = 3.09, F, = 2.64, K = 3.4, and
RD = 530. U n
-------
41
z_
n
Plume over the Terrain Model
Plume over a Flat Platt
x/0
+O-4 ~
+ 32 -
-32 -
-96 -
! I I I I I I I I 1 I
I I I I I I I
128 160 192
256 288 32
L' = 6cm/s
4- Mountain Pea
-)- Mean cam Pea
Figure 19. Vertical and Lateral Spreads of Plumes
In Stratified Flows over a. Plate Plate
and the Terrain !Icdel. F = 3.16,
F;, = 2.64, F = b.22, n =D3-n, and
RD = 530.
-------
42
—
-
-
-
oo
0
fj^ P!
— I— t
1 1 1
Lower Boundary
Terrain model
— w
j— [
CD
III!
O
' ' •£ ' «Q U l I i i I i i
jS «C D
\ •< D
$ •< n
? ^ % £
* «3 cP ^
* « n e
^t tJ v 1 o —
F* PS 0
. « ^'— d>
^ • < D
D
• < a
^ • ,
co
4J
Q
X
Oi
E
te
ce
0;
c
C
c
CO
C
1—I
tu
o
eg
4)
M
M
-------
43
D
w
OJ
o
tn
O
m
cx
o
4J 13
Oi
C •<"•'
O U-l
•H
en 4-j
4-1 Cfl
CJ !-i
01 4-J
ca
13
§
o
CQ
-------
44
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
V
0 3
\ • I *<
• ' < R
^ • ^ u
"" «& • < a ~
., D
** • < D
• < D ~
* a
.fc • < O
D
• < a
*
oo -a- oo «* >r D
^ • » • » F^
-Si CM -3" CM -s _
>-l rH iH »*
cd cu a)
~ (3 O O -i w w
U M >-i cd cd
O CU CU rH ^-t Jy
J H H fn Fn T*
^ D -
j§ ^ • < a
1 i i i i i i i i 1 i i i i i i
00 Or
O O ^H
>T
x^\
vO
r-l
m
r-l H
\~s
CO
I
rH
P4
M-l
O
13
cd
cu
£
rH
cd
CU •
4-1 CO
X 8 cd &
H
01 PK
CU
a -H
O T4-I
•rl
CO 4J
4-1 Cd
r-l OH
01 4J
14-1 C/3
>4-l
•H "cd
cd -t-i
CM
CM
1
r<4
I— 1
O
CN
-------
45
a.
f= 1.32
UOT = 3cm/s
r50
-40
-30
— 20
-10
— o
/— •>
jl
i
N
b. F, = 2.64
UOT = 6cm/s
r50
-40
-30
-20
-10
— 0
w
I
N
Figure 23. Dye Visualizations of Bent-over Plumes (FD = 3.09) from a Stack
Located in the Lee of the Terrain Model.
-------
o o
r-i •
ft* ro
O
l-i C
14-1 P-i
0) CO
C *J
•g S
£ B
01 01
4-> M
0) 3
P CO
rt
CO 0)
CU
M
00
| P
-------
47
VD
Q
X
00
C -T
•H <—'
tt •
(L ~-^
^ >
O ft
E -
C -
o
-------
48
max
W
.1
.09
.03
.07
.06
.05
.04
.03
.02
I I I
10
i l l
. . I
20 30 40 50 60 80 100
JC
D
y/D - 3.14
200 300 400
max
.015
.010
.008
.006
.005
.004
.003
.002
10
\ -
Slope
O
I l
l l l i l 1
20 30 40 50 60 80 100
y/D - 0
200 300 400
D
Ficure 26. Profiles of Maximum RMS Velocity and Density Fluctuations in the Plume
O _ _ _ , _ —. i TI *-»rvr\T«— l£/. V ^ < /
with Downstream Distance from the Stack. FD - 3.09,
and R - 530.
2.64, K - 3,4 ,
-------
z-Z
u
(AZu/2)
-1
-2
_0
AZ (cm) Z (cm) / ^ (cm/.'-
i, ' 11 mo v
max
.8
1.0
max
1.73
.85
0.60
a.
z-Z
(AZa/2)
-1
-2
-3
y/D = 0
AZ (cm) Z (cm)
a o
'A
a (cm/sot
max
1.45
.65
1.05
.6
.8
1.0
b.
Figure 27. Profiles of Normalized RMS Velocity and Density
Fluctuations in the Plume. F
K - 3.4, and R - 530.
3.09, F^- 2.64,
-------
50
z-Z
(AZu/2)
-1
-3
y/D = 3.14
x/D AZ (cm) Z (cm) / u (cm/sec)
u u max
14.3
15.7
16.3
1.73
.85
0.60
10
a.
z-Z.
(AZ0/2)
1
0
-1
-2
-3
y/D = 0
AZ0(cm) Za(cm)
—T"
a (cm/sec)
max
1.45
.65
1.05
12
16
20
a4/(a2)2
b.
Figure 28. Flatness Factors of Velocity and Density Fluctuations in the
Plume. F
D
3.09,
2.64, K = 3.4, and R = 530.
-------
TASK II
A FEASIBILITY STUDY OF NUMERICAL SIMULATION
OF PLUME DISPERSION IN STABLY STRATIFIED
FLOWS OVER A COMPLEX TERRAIN
-------
A Feasibility Study of Numerical Simulation of Plume
Dispersion in Stably Stratified Flows over Complex Terrain
by
Flow Research, Inc., Kent, Washington 98031
I. Introduction
Flow Research, Inc. has considered several alternative approaches to
the numerical modelling of realistic flow fields with simple source
diffusion in the vicinity of complex terrain. Our approach to this
problem has followed several diverse paths, all of which we believe add
to give a comprehensive modelling approach that is most appropriate to
guide future work.
In Part A (Flow Research Note No. 40, "Calculation of Stably
Stratified Flow Around Complex Terrain" by D. K. Lilly), a method is
proposed for the calculation of strongly stratified flow in complex
terrain. This method is expected to be useful for modelling single
source diffusion of pollutants when the source lies below the mountain
tops. Extensions of this model are proposed to handle the effects of a
turbulent atmospheric boundary layer, as well as mountain lee waves.
The method proposed in Note No. 40 involves a Froude number expansion
of the steady state Boussinesq equations of motion. The Froude number
e is defined by
e = U2/h2N2
where h is the mountain height, U is a characteristic horizontal velocity,
and N = - TT^ is the square of the Brunt-Vaisala frequency. If the
Po a Z
potential energy required to lift fluid to the mountain top is much larger
than the kinetic energy of the fluid then £ is very small. It seems that
in cases of strong stratification, it suffices to truncate the Froude number
expansion at first order. This expansion technique applies directly to the
flow field. Using K-theory and the resulting flow field, it is then possible
to calculate the single source diffusion of a pollutant in this strongly
stratified environment. This method will be explored further under the
current EPA Contract No. 68-02-1293.
-------
We have also studied the feasibility of pt-forming spectral computations
in complex terrains. A brief description of a possible technique for doing
this follows. Suppose the terrain is given by z = f (x,y), where z is the
vertical height of the ground. We suppose that the top of the computational
domain is z = z™' independent of x,y (for simplicity). We then introduce
new coordinates
y' = y
zT~f (x,y)
so that 0 j^ z' <_ 1. In terms of these coordinates, derivatives are
transformed into
3f(x1.y<)
3F(x,y,z) = JL _ - — l-x!__ UL
_ -
3x 3x' z -
3f(x',yl)
3F(x.y,z) _ 3F 3y' 3F
f\ I \J-Zy c/l TN^T
3y 3y' z_-t(x ,y') dz
1
3 F(x.ytz) = ZT 3F
8z zT-f(x',y') 3z'
Since the functions
ZT -
and
can be easily stored a priori in their spectral forms, it follows that
evaluation of derivatives, nonlinear terms, etc., can be easily implemented
using the techniques already extensively developed by Orszag and others.
Spectral simulations of the type discussed here, when applied to the
Boussinesq time-dependent equations of motion, have several advantages
over finite-difference simulations of the same flows, especially the absence
of phase errors in the propagation of linear theory waves. These advantages
-------
appear whether or not f(x,y) is smoothly continuous, but they are more
obviously realized when f(x,y) is smooth. Such spectral simulations have
the advantage over the Froude number expansion discussed above of
applying in both stable and unstable environments. In addition, the
effects of atmospheric turbulence can be explicitly taken into account
by means of turbulence transport models or sub-grid scale nonlinear
eddy viscosity coefficients. These full spectral simulations together with
associated single source diffusion calculations can provide a complete
set of atmospheric data for the evaluation of environmental impact of
proposed installations.
Flow Research, Inc. has also made a careful study of techniques for
the imposition of radiation boundary conditions in numerical calculations.
In Part B (Flow Research Report No. 30, "Numerical Simulation of Radiation
Boundary Conditions" by M. Israeli and S. A. Orszag), several methods are
discussed for the imposition of radiation conditions. This problem is
particularly important in the stratified atmospheric environment in several
respects. First, numerical models of atmospheric motions are generally
performed with some kind of lid affixed to the top of the atmosphere.
If this lid reflects internal waves back into the computational domain, the
computations can be degraded and serious inaccuracies can result. The
problem has previously been treated satisfactorily only for the over-
simplified case of monochromatic waves. In the general case, various
prescriptions for imposing "spongy" boundary conditions are used — it is
hoped that the progress reported in Report No. 30 will be useful in the
solution of this vexing problem in atmospheric dynamics. Second, the
lateral boundaries can also reflect internal waves in a stratified fluid giving
spurious results. Again, the appropriate radiation boundary condition should
be applied. In summary, it is exceedingly important for numerical simu-
lations of stratified flows to make use of the proper boundary conditions,
otherwise improper effects of the stratification, like downstream as well
as upstream blocking, may result. With the boundary conditions of the
type discussed in Report No. 30, these difficulties are expected to disappear.
-------
TASK II
PART A
CALCULATION OF STABLY STRATIFIED FLOW
AROUND COMPLEX TERRAIN
BY
DOUGLAS K, LILLY
-------
A-l
Calculation of Stably Stratified Flow
Around Complex Terrain
by
Douglas K. Lilly
Flow Research, Inc., Kent, Washington 98031
The purpose of this note is to derive and demonstrate a method for
calculating the flow around and through complex terrain when thermal strati-
fication is sufficiently strong that a part of the flow cannot pass over
the terrain. It is suggested that this method might be useful for deter-
mining the transport of pollutants introduced into the atmosphere above
the boundary layer but well below the tops of terrain obstacles. The
application may be somewhat more general if a turbulent boundary layer
model is incorporated into the lower part and a simplified wave model can
be fitted to the top. From a more fundamental point of view the results
show how the difficulties associated with blocking of flow across two-di-
mensional obstacles essentially disappear when a more realistic three-di-
mensional terrain is introduced.
It is known that when the kinetic energy of a stably stratified fluid
stream impinging upon an obstacle is smaller than the potential energy re-
quired to lift it over the obstacle, the flow will be blocked to some
degree upstream of the obstacle. The basis of the present method is
the assumption that each layer of the blocked flow resembles two-dimensional
(horizontal) potential flow around the contours of the terrain at that
level, and that the deviations from two-dimensional potential flow can
be deduced as higher order terms in a convergent series expanded about a
small parameter.
In the following we develop the formal series solution and illustrate
its application by means of an analytic solution for a simple terrain shape.
More realistic solutions would require numerical methods and also should
take account of the separation of the quasi-potential flow from the back
side of the mountain as well as boundary friction and heat transport effects.
-------
A-2
The steady state incompressible inviscid Boussinesq equations of
motion, continuity, and conservation of buoyancy may be written in dimen-
sionless form as
£ • ° «
2
h .. , 3w , 3w . 9w
3x 3y 3z
U2 ,
-w = 0 " (5)
The variables have been scaled as follows. Horizontal coordinates x
and y by L, a characteristic horizontal mountain dimension; vertical
coordinate z by h, the mountain height; horizontal velocities u and v
by U, a characteristic horizontal velocity below the mountain top level;
vertical velocity w by hU/L; disturbance pressure p - p (z) by Pc>u2 ;
_ 2
disturbance density p - p (z) by U po/gh, where p is a reference
density. The pressure and density scaling were chosen so that the pressure
disturbance would be of the same nominal magnitude as the horizontal
kinetic energy fluctuations, and the density scaling so that the hydro-
static terms in (3) would be of the same nominal magnitude.
The above scaling introduces two dimensionless parameters. The
first, in Eqn. (3), is simply the ratio of the height to the width of
the mountain, and if that is small it suggests use of the hydrostatic
equation, which we assume to be appropriate. The second parameter, in
Eqn. (5), is an internal Froude number of inverse gross Richardson number,
where N =-%— —H . This Froude number is essentially the ratio of the
Po 9Z
-------
A-3
mean flow kinetic energy to the potential energy required to lift fluid
over the mountain from the ground. In cases of significant blocking the
parameter is small. Equation (5) suggests an expansion solution in
powers of the Froude number, with the lowest order solution a purely hori-
zontal flow. Thus we assume the expansion as follows:
w = e w
(6)
P - P(1) + e
where e - U2/h2N2.
Upon substitution of the above expansion into Eqns. (1) - (5) we obtain
the first order equations as
a
3x
+ p(1)) - -a'1' C (8)
(9)
-------
A-4
If we multiply (7) by u , (8) by v and add, it is seen that the Bernoulli
equation is satisfied along a streamline of the first order velocity. Also by
differentiating (7) with respect to y and (8) with respect to x and subtracting,
we can obtain a similar equation for conservation of first order vorticity.
Since the two-dimensional continuity equation, (10), allows the formation of
a stream function fy , these relations may be specified as
+ P<1> - function OJ,(1)) (ID
and C(1) = vV1} = function Oj,(1)) . (12)
If an upstream velocity distribution were specified it "Touid be possible to
solve (12) for the flow, (11) for the pressure, and (9) for the density. In
a simpler case, however, we simply assume that the incoming free stream flow
is irrotational, so that the solution required is that of potential flow around
the mountain contour at each z-level independently. Note that the free stream
velocity may be an arbitrary function of height, and not dimensionless unity
as implied by the scaling of Eqns. (1) - (5). Note also that the addition of the
Coriolis parameter to the vorticity in Eqns. (7), (8), and (12) need not cause
any significant change unless the horizontal scales are so large that the variation
of the Coriolis parameter with latitude must be considered.
The most important second order equation is that arising from substituting
definitions (6) into (5), the buoyancy conservation equation. This leads to
the relation
w<2) . „(!) 3 + ,(D (13)
From this, the divergent part of the second order horizontal velocity may be
obtained from integration of the continuity equation and the rotational part and
the second order pressure from the second order contributions to (1) and (2).
To show results of this solution from an analytically solvable case we choose
a constant free stream velocity equal to unity and a mountain with circular
contours, i.e.,
2 o
x + y = F(x) at the mountain surface (14)
-------
A-5
Then the irrotational flow solution around such a mountain will be given by the
derivatives of the velocity potential , which may be obtained conventionally
by assuming a source sink pair at the origin x=y=z=0. The potential function
and the velocities are given by
(15)
- cos9 (r + F/r)
- sine
where r is the distance from the origin and 0 the angle in cylindrical coordinates.
From Eqns. (11) (with the function equal to 1/2, that is half the free stream
velocity squared), we then obtain p , and from (9) and (13) p and w as
follows;
pvi' - -?- (cos 26 - F/2r2) (18)
- - (cos 20 - F/r2) (19)
w(2) . 2FJ. (cQ8 3Q _ 31 cose + 2L_ cosQ ) (20)
r3 r2 r
The higher order terms for the divergent and rotational velocities become
considerably more difficult to compute, analytically at least.
2
At the mountain surface F = r • From this we can see the general trend of the
solution. Eqn. (17) shows that the first order flow vanishes at the normal
stagnation point and is twice the mean stream velocity at the edge. The pressure
is high at the stagnation point and low at the edge, as expected. The disturbance
density, and therefore the vertical displacement, vanishes at the stagnation
point (but not upstream of it) and is negative at the edge (if F' < 0, corresponding
to a mountain width decreasing with height). The vertical velocity also vanishes
-------
A-6
at the stagnation point and at the edge, but is first upward then downward upstream
of the stagnation point and downward in the acceleration region.
It is not obvious whether a zone of convergence can be readily established for
this solution. An idea of the probable zone of convergence can be obtained by
looking at the conditions under which the first order stratification,9p /3z,
is of the same magnitude as the mean stratification. At such a location the first
order streamlines may cross and the solution becomes unrealistic. This can be
shown to occur when |9p^/9z = e'1, i.e. when e = r2/ (F" (cos26-F/r2)-F' 2/r2) | .
2
It always occurs first at the mountain surface, F = r , but the details depend
on the precise form of F chosen.
2 2
For an ellipsoidal (spherical in the scaled coordinates) mountain, r 4- z =1.
If e = 0.1 the probable zone of convergence is below z = 0.51. For e = 0.2 it
is below z = 0.27, and no convergence can apparently be obtained if e > 0.25.
2
On the other hand for a conical mountain, r + z = 1, the limiting condition is
2
simply e = (1-z) , so that there is a probable zone of convergence for all values
of e < unity. The meaning of the non-convergence at higher levels is not
completely clear but is must be associated with the fact that a complete solution
of any mountain wave problem requires an upper boundary condition. The
apparent convergence of the above solution close to the ground implies that
the solution there is independent of the upper boundary condition. It may
be speculated that the upper part of the present solution could be used as a
boundary condition for a linearized lee wave solution in the upper regions but
there may be an intermediate region where neither approximation is sufficiently
valid.
-------
TASK II
PART B
NUMERICAL SIMULATION OF RADIATION BOUNDARY
CONDITIONS BY DAMPING
BY
M, ISRAELI
STEVEN A, ORSZAG
-------
Table of Contents
Page
Table of Contents i
1. Introduction B-l
2. Steady Oscillations B-6
3. Initial - Value Problems B-ll
4. Reference B-12
5. Figure B-13
-------
B-l
Flow Research Report No. 30
Numerical Simulation of Radiation Boundary
Conditions by Damping
by
M. Israeli and Steven A. Orszag
Flow Research, Inc., Kent, Washington 98031
1. Introduction
Radiation boundary conditons appear in a wide variety of physical
problems involving wave propagation. With mixed initial-boundary value
problems of hyperbolic partial differential equations, boundary conditions
should be specified for dynamical quantities propagating on characteristics
entering the domain of integration, while no boundary conditions are
necessary for quantities propagating on characteristics leaving the domain.
If the domain of integration is finite, numerical solution of such a problem
may be effected by standard techniques, such as the finite-difference
method. On the other hand, if the desired region of integration is
infinite, numerical solution encounters difficulty because of the necessary
limitation of the computational domain to a finite region. With initial-
value or radiation problems, the appropriate boundary conditions at infinity
are normally "radiation boundary conditions"; that is, the amplitude of
waves entering from infinity are required to be zero, while no conditions
areplaced on outgoing waves propagating to infinity. With scattering
problems, the amplitude of an incoming wave from infinity is specified
and the amplitudes of the "scattered" outgoing waves emanating from the
physical domain are sought; scattering problems can be easily reformulated
as radiation problems. The difficulty with numerical solution of these
problems is transferring the desired boundary condition at infinity to
a condition at the finite boundary of the computational domain. Many
important problems await proper resolution of the above difficulty, in-
cluding flow past a body or obstacle in a stratified fluid.
The mathematical difficulty is best brought out by the simple problem
of the one-dimensional wave equation
utt - Uxx (Olx<">
-------
B-2
with radiation boundary conditions at °°,
u(0,t) = f(t) (t > 0)
and the special initial conditions
u(x,0) = 0
ut(x,0) = 0 .
(0 < x <
The exact solution to this problem is
(0 t < x
u(x,t) = j (2)
(f (t-x) . t >_ x
The difficulty with numerical solution of the problem is that if
arbitrary boundary conditions are placed at any finite boundary used
to simulate the infinite region then arbitrary and wrong results are
gotten. For example, if the boundary condition u = 0 is applied at
x = L, the solution to (1) becomes
u(x,t) = tf(t-x-2nL) - f(t+x-2n+1L)], (3)
n=0
where f(t) is extended to negative t as f(t) = 0 for t < 0. The solution
(3) is identical to (2) only for 0 <_ t < L. Beyond t = L, waves are 'reflected'
off the boundary x = L and there is a finite discrepancy between (3) and (2) .
When numerical solution proceeds using such boundary conditions at x = L,
the simulation of (2) can be faithful everywhere only for t < L, and near
x = 0 only for t < 2L. For t > 2L, the reflected waves reach x = 0 and
completely change the solution there from the desired outgoing solution.
Of course, the numerical solution is useful for t < L so that if L is
large enough, the solution remains useful for a long time. However, in
practice , L can not be taken very large (especially in two and three
dimensional problems) because of limited computer resources so that the
solution obtained within the finite domain is completely unsatisfactory.
For the problem (1), a satisfactory solution to the above difficulty
would be application of the radiation boundary condition
u + u = 0 at x = L . (4)
-------
B-3
With (4), the solution to (1) within 0 £ x £ L is precisely (2);
waves do not reflect from x = L with the boundary condition (4). The
problem (1) may be solved satisfactorily using numerical methods with
(4) applied at x = L. However, this does not solve the general problem
of radiation boundary conditions. For example, with the two-dimensional
wave equation in polar coordinates,
utt = Urr + r Ur + \ U89 '
the generalization of (4) at the finite outer boundary r = L is
u + u = 0 (6)
However, with (6) errors of order 1/r are made in the radiation condition,
so that unless L is very large, the errors due to (6) will limit the
accuracy of the numerical solution to (1). Further difficulty with boundary
conditions like (4) or (6) is anticipated in more difficult problems where
either the linear wave dispersion relation remains intractable or non-
linearities, strong inhomogeneities, etc. are present.
Another method for handling the problem of an infinite domain that
may at first seem appropriate is a coordinate transformation of the in-
finite domain to a finite domain or, nearly equivalently, the use of a
nonuniform grid in finite difference approximation to (1). For example,
the coordinate transformation
z = 1 - exp (-x) (7)
transforms 0 £ x < °° into 0 £ z < 1. Equation (1) becomes
tt zz z
It would seem that numerical solution of (8) on the uniform grid z = n/N,
n = 0, ..., N would give accurate solution of (1) with the radiation
conditions. In effect, this would amount to use of the nonuniform grid
1
xn
n
directly on (1).
-------
B-4
Unfortunately, however, this method cannot work in general. It is
well known (cf. Orszag and Israeli 1974) that difference approximations
to a wave problem require high resolution per wavelength of the solution
in order to be accurate; typically, order 10 grid points per wavelength
are required for, say 10% accuracy. Thus, if the typical wavelength of
the exact solution is X, we require grid spacing Ax < 1/10 X or, from (9),
which is violated when
Thus, the waves cannot be adequately resolved near z = 1. What then
happens to the resolution depends on the numerical scheme; in general,
inadequately resolved waves reflect from the grid. In essence, if the
wave between grid point x and x+h is not resolved because the wave oscillates
too rapidly then the coupling between the grid points is like that between
two weakly coupled mass points one of which is oscillating rapidly. The
oscillation in one does not effectively drive an oscillation in the other;
for the present problem there is a direct analog to the boundary condition
u = 0 at x = L or (1), so that reflected waves result.
The trouble with mapped coordinates can be seen already from (8) where
z = 1 appears as a singularity of the spatial part of the equation. Mapping
just sweeps the radiation condition problem under the rug, where it can
still befoul numerical solution.
The purpose of this paper is to investigate the application of still
another technique for imposing radiation boundary conditions, namely damping.
If damping is imposed on all waves propagating near the outer boundaries
of the finite computational domain, then the outgoing waves are reduced in
amplitude when they hit the outer boundary and reflect, while the reflected
waves are reduced in amplitude as they propagate in from the outer boundary
to the region near x = 0 of most interest. By properly 'tuning' the damping
to the problem, it is hoped to minimize the errors from the reflected waves
and thus, permit faithful numerical solution of the problem. Of course,
-------
B-5
introduction of weak damping is the standard mathematical artifice for
analytic solution of steady state wave problems. With any (arbitrarily
small) damping, incoming waves from infinity are damped to zero amplitude
by the time they reach any finite region of space, while outgoing waves
are but slighly damped in any given finite region if the damping is small
enough. In the limit of zero damping, the solution satisfying the desired
radiation boundary condition is obtained.
-------
B-6
2. Steady Oscillations
The most transparent way to impose damping on (1) is to first rewrite
it as the system
v = w
t x
W = V
t X
where v = u , w = u . With second order (Newtonian viscosity) damping,
L X
the system is
v = w + v(x)v
t x xx (1Q)
W = V
t X
We consider the solution of (10) on the interval 0 <^ x <^ L with boundary
conditions
w(L,t) = 0
ikf-
v(0,t) = sin kt = Im(e ) (11)
A ikt
If we look for solutions of the form v = Im(v(x)e ), it follows that
v + . * , . v = 0 (12)
xx l+ikv(x)
^ —ikx
When v = 0, the outgoing wave solution is v = e
When v is constant, there are solutions to (12) of the form exp(iax)
where
a - + k// 1+ikv = + a . (13)
The boundary conditions (11) are satisfied by the solution
^ _ cos a (L-x)
cos a L
The error relative to the exact outgoing wave exp (-ikx) is
e(x,t) = Im [(v(x) - e
-------
B-7
so that the maximum error at any x is
e (x) = '"' * ~ikx
max
Several limiting cases are of interest:
!_ . 2
(i) Vk « 1, Vk L « 1. In this case a « k - 2 so that
v(x,t) %, sin kt cosk(L-x)/coskL, (15)
which is just the same standing wave obtained without damping.
In this limit, the effect of damping is everywhere negligible.
1 /"k
(ii) vk » 1. In this limit, a % j / — (1-i), so that the wave
A V v
is distorted. In the limit k L2 » kv, it follows that
v(x) * expf|/| (i-l)xj (16)
2
(iii) vk « 1, vk L »1. In this limit, the maximum error is
j 2
e (x) ft/ 1 - e 2 , .
max b , (17)
so that accurate simulation of the outgoing wave can be achieved
on a fixed interval 0 < x < x provided x « L. In fact,
— o o
the conditions of case (iii) can be reinterpreted as follows.
2
First, Vk « 1 can be rewritten vk « k so that the requirement
is that the damping rate be much smaller than the wave frequency.
2
Second, the condition Vk L » 1 requires kL » 1, so that the
number of wavelengths in the computational domain is large. The
latter requirement together with the condition x « L shows that
the constant damping is a very poor way to achieve radiation
conditions. In fact, a more careful analysis of the maximum
error than provided by (17) would show that there are two essential
kinds of error in v(x) which may be termed phase error and damping
error. Phase error is that error induced or the frequency of the
outgoing wave by the presence of the damping, while the damping error
is due both to damping of the outgoing wave and the unwanted
reflected wave. In the present limit, phase error is negligible.
-------
In order to make error of order e in the region 0 £ x
necessary that
while
x , it is
-Vk2L
(18)
(19)
In other words, we require
, vkx ^ 2e
With e = .01, k = x =1, it follows that we must choose L ^ 250,
V ^.02. Thus, 99.6% of the computational domain is "wasted", hardly
a satisfactory situation.
The above analysis can be extended to other damping factors. With
various v(x), (12) can be solved analytically in terms of special
functions like Airy functions, hypergeometric functions, etc. However,
the most general results are obtained by WKB analysis. The WKB solution
to (12) is
l-ikv(x)
l-ikv(O)
Error less than e for 0 <_ x <_ x requires
cos
COS
kf ds
/
L -Wl-ikv(s) J
kf ds
/
Wl-ikv(s)
J v (x) dx < 2 E
o
L
J v(x) dx > In i
(20)
(21)
which are the generalizations of (18), and (19).
For example, if we wish to choose the optimal damping factor with
v(x) = Vxn for given ratio of useful to total computational domain, x0/L,
-------
B-9
then we must satisfy
1 Vk2 xn+1 * 2£
n+1 o
1 , 2T n+1 n , 1
Vk L ^ In —
n+1 e
so that
(22)
In other words, as n increases, the maximum error on 0 < x < x
— — o
decreases (at least while the above approximations are valid). For
example, with x = -^L and n = A then e'v .01. With x = k = 1, it
follows that V^ .06, L ^ 3.
The condition for applicability of the WKB approximation is
dX .^ .. -v 2iT r:—n—-,—r
— «1, X = — /l-ikv(x),
or with v(x) = Vx , nvx « 1. If we require this condition everywhere,
then we must have
nv L11'1 « 1. (23)
This condition is not satisfied by the above example. However, the weaker
condition nVx11 « 1 for validity of WKB on 0 < x < x is satisfied. It
o J — — o
may be shown by lengthy and complicated turning point analysis of the WKB
solution that this weakened condition is sufficient for the present
application.
We have studied the solution to (12) numerically. Rewriting (12) as
vxx + a2(x)v = 0,
we used the fourth-order difference approximation
(A(x+h)+l) v(x+h) + (10A(x)-2) v(x) + (A(x-h)+l) v(x-h) = 0 (25)
2 o
where A(x) = [a(x)] h /12. The boundary conditions on (25) are
•j
9v/3x = 0 at x = L, v(0) = 1, the former being applied with 0(h ) accuracy
by Taylor expansion.
-------
B-10
The calculations were performed with various k and v(x) on the
interval 0 _£ x j£ rr (i.e., L = TV) . The results for the maximum error
£ (x) are shown in Fig. 1. For example, a maximum error of 4% on
max .^ ^
the interval 0 _£ x ^ x = ~TT is achieved by v(x) = .05x when k = 3.
Among all coefficients in v(x) = Vx ,v = .05 gives the smallest maximum
4
error. For k = 6, the error can be reduced to 2% withv= .004x over
the same interval. Again the coefficient V = .004 is optimal. According
to (22), x = -rL should give 1% errors when n = 4.
o j £
For larger values of x /L, the errors are larger. With v(x) = .04x ,
the error is 21% for x /L = — . (The latter example was chosen at
o 2.
random with no effort given to optimization).
The present calculations are meant to be just an illustration of
how best choices for v(x) can be made. More detailed analyses will be
presented later.
-------
B-ll
3. Initial-Value Problems
The above analysis on the effect of damping in steady-state problems
can be misleading for time-dependent computations because of slowly
decaying transient effects.
Consider the case of v(x) = V so that
u = u + Vu (26)
tt xx xxt
with initial conditions
u(x,0) = 0
ut(x,0) = 0
and boundary conditions
u(L,t) = 0
u(0,t) = sin kt .
This problem may be solved easily by Laplace transform in t giving the
Laplace transform
sinh p(L-x)//l+py
„, x _
U(x,p) -
~ „ . //-, . .
p2+k2 sinh pL//l+pV
The poles of U(x,p) at p = + ik give the persistent part of the
solution as t -»• °° , namely (14). However, U(x,p) has other poles at
(n ^ 0) (28)
/1+pv
P = - (Vn2TT2 + /vn\4 - 4L2n2TT2)/2L2
*• * -
TO 2
the latter holding when vn IT « 4L . The decay rates of these transient
92 2
modes are VnTr /2L . In particular, if L = TT, the most persistent transient
persists for the time 2/v, which may be very long. In fact, this effect
gives much larger errors in time-dependent problems than may be expected
on the basis of the steady state analysis in Section 2. We have made
-------
B-12
several different kinds of numerical experiment to verify this behavior.
Indeed, we found much larger errors than could be explained on the
basis of the steady state analysis.
In summary, our analysis indicates the drawbacks and pitfalls encountered
in imposition of radiation boundary conditions. To date, we believe that
the best approach is to use appropriate boundary conditions like (4) and
(6) on radiation boundaries, together with optimized damping, as
explained in Section 2.
Reference;
Orszag, S. A., and Israeli, M. 1974. Ann. Rev. Fluid Mech. j6, xxx.
Figure Captions;
Fig. 1 Spatial dependence of the error with damping factor v(x).
-------
B-13
2w/3
Fig. 1 Spatial dependence of the errrr with damping factor v (x)
-------
DO
r—
CD
CU 03
Q. —
Q. -
D
Is
O CD
* <
(D (V
m 3
- (0
3 2
O. TU
ID
2 in
m -^
Z t>
*>;>
• 32
D
U» i -n
"J m
O m
>0
o
------- |