Skip to main content

Full text of "Numerical solutions of atmospheric flow over semielliptical simulated hills"

See other formats


I 



NASA Contractor Report 3430 




Numerical Solutions of Atmo^pii^ric^i^^ S 
Over Semielliptical Simulated Hills 



Chih Fang Shieh and Walter Frost 



CONTRACT NAS8-32692 
JUNE 1981 



fVIASA 



TECH UBRARY KAFB. NM 

inn 




NASA Contractor Report 3430 



Numerical Solutions of Atmospheric Flow 
Over Semielliptical Simulated Hills 



Chih Fang Shieh and Walter Frost 

The University of Tennessee Space Institute 

Tidlahoma, Ten n essee 



Prepared for 

Marshall Space Flight Center 

under Contract NAS8-32692 



(M/NSA 

National Aeronautics 
and Space Administration 



Scientific and Technical 
information Branch 

1981 



I II nil lllll III 1 111111 I I nil III! I III I II iiHi^iiniiiiMiiB^^^^^^m wnniB^Hn 



ACKNOWLEDGMENTS 

This research was supported under NASA Contract No. NAS8-32692. 
The authors are grateful for the support of A. Richard Tobiason of 
the Office of Aeronautical and Space Technology, NASA Headquarters, 
Washington, DC. Special thanks go to Dennis W. Camp of the Space 
Sciences Laboratory, Atmospheric Sciences Division, NASA/George C. 
Marshall Space Flight Center, Alabama, who monitored the research 
program. 



n 



ABSTRACT 

Analysis of atmospheric motion over obstacles on plane surfaces 
is carried out in the present study to compute simulated wind fields 
over terrain features. Emphasis is on a semielliptical , two-dimensional 
geometry. Numerical simulation of flow over rectangular geometries is 
also discussed. The numerical procedure utilizes a two-equation turbu- 
lence model and the selection of the necessary constant coefficients in 
J 

the model is considered an important part of the present study. 

In the present approach the partial differential equations for 
the vorticity, stream function, turbulence kinetic energy, and turbu- 
lence length scale are solved by a finite-difference technique. The 
numerical solutions have been compared with available experimental data; 
agreement is good. 

It is found that the mechanism of flow separation induced by a 
semiellipse is the same as in the case of flow over a gradually sloping 
surface for which the flow separation is caused by the interaction 
between the viscous force, the pressure force, and the turbulence level. 
For flow over bluff bodies, e.g., solid fences, a downstream 
recirculation bubble is created due to the inability of the flow to 
negotiate with the abrupt change of the surface shape. Increasing the 
aspect ratio and/or increasing the turbulence level results in flow 
reattachment close behind the obstacle. 



m 



TABLE OF CONTENTS 

SECTION PAGE 

1.0 Introduction 1 

2.0 Review of Previous Works 6 

3.0 Turbulence Models 12 

3.1 Mean Flow Equations 12 

3.2 Closure Problem 13 

3.3 Concept of Eddy Viscosity 14 

3.4 The Two-Equation Models 15 

3.5 K-L Two-Equation Model 16 

3.6 Constants of the K-L Model 21 

4.0 Functional Equations in General Orthogonal Curvilinear 

Coordinates 30 

4.1 Continuity Equation 31 

4.2 Conservation of Momentum Equations 32 

4.3 Transport Equations for Turbulence Kinetic Energy and 
Length Scale 39 

4.4 Vorticity and Stream Function Equations 42 

4.5 Analysis of Flow Parameters 44 

5.0 Numerical Algorithm 48 

5.1 Derivation of Finite Difference Equations 48 

5.2 Elliptical Cylinder Coordinate System 55 

5.3 Numerical Coordinate System 58 

5.4 Boundary Conditions 59 

5.4.1 Outer Boundary 59 

iv 



SECTION PAGE 

(Continued) 

5.4,2 Wall Boundary 62 

5.5 Accuracy of Numerical Solutions. 64 

6.0 Results and Discussion 73 

6.1 Characteristics of Flow About a Two-Dimensional 

Semi elliptical Obstacle 73 

6.2 Effect of Aspect Ratio 81 

6.3 Effect of Surface Roughness 84 

6.4 Verification of Results 93 

6.5 Phenomenon of Flow Separation 107 

7.0 Conclusions 113 

Bibliography 115 

Appendix 119 



LIST OF FIGURES 

FIGURE PAGE 

3.1 Effect of the ratio of constants on the distributions of 

K/K and L/L behind a rectangular block 25 

4.1 An orthogonal curvilinear coordinate system 30 

4.2 An axi symmetric coordinate system 33 

4.3 Velocity vectors of a fluid element shown in a constant 

U3 plane 34 

5.1 Grid arrangement with respect to the central node P 50 

5.2 The numerical and physical coordinate system for 

calculating the problem of flow over 

semi ellipses 56 

5.3 Distribution of grid sizes for D/H = 2.0, C-j = 0.025, 

C2 = 1.05 60 

5.4 Streamline, -ip, patterns of flow over a semi ell ipse for 

D/H = 2.0 66 

5.5 Distribution of turbulence kinetic energy, K, for D/H = 2.0. . 67 

5.6 Distributions of turbulence length scale, L, for D/H = 2.0 . . 68 

5.7 Vorticity, w, patterns for flow over a semi ell ipse for 

D/H = 2.0 69 

5.8 Velocity profiles on top of a semiellipse for D/H =2.0 

and z^ = 0.005 71 

6.1 Streamline, ijj, patterns over a semiellipse for D/H = 1.02 

and z = 0.005 75 

6.2 V/Vn isotachs over a semiellipse for D/H =1.02 and 

z = 0.005 76 



vi 



FIGURE PAGE 

6.3 Distribution of turbulence kinetic energy, K, over a 

semiellipse for D/H = 1.02 and z = 0.005 77 

6.4 Distribution of turbulence length scale, L, over a 

semiellipse for D/H = 1.02 and z = 0.005 79 

6.5 Vorticity, w, patterns over a semiellipse for D/H =1.02 

and z = 0.005 80 



6.6 V/Vu isotachs over semi ellipses for z = 0.005 82 

6.7 Velocity profiles on the top of semi ellipses of different 

aspect ratios for z = 0.005 83 

6.8 Velocity profiles on top of a semiellipse with different 

surface roughnesses for D/H = 2.0 85 

6.9 V/Vl, isotachs over a semiellipse for D/H = 2.0 and 

z„ = 0.005 86 

6.10 Streamline, ij;, patterns over a semiellipse for D/H = 2.0 

and z = 0.005 87 



6.11 Distribution of turbulence kinetic energy, K, over a 

semiellipse for D/H = 2.0 and z^ = 0.005 88 

6.12 Distribution of turbulence kinetic energy, K. over a 

semiellipse for D/H = 2.0, z^ = 0.01, and Zj^ = 0.01 90 

6.13 V/Vjj isotachs over a semiellipse for D/H = 2.0 91 

6.14 Wind speed-up ratio on top of semiellipses with different 

surface roughnesses and different aspect ratios 92 

6.15 Wind speed-up ratio over a 90° escarpment 94 

6.16 Turbulence intensity over a 90° escarpment 96 

6.17 Effect of the turbulence of the approaching wind on the size 

of the wake zone behind rectangular blocks 99 

vii 



FIGURE PAGE 

6.18 V/V isotachs downstream of a solid fence, z =0.01 lOO 

6.19 Various shapes of surface obstacles for which either 

analytical or experimental data are compared in 

Figure 6.20 102 

6.20 Vertical velocity profiles on the top of the respective 

obstacles shown in Figure 6.19 103 

6.21 Distribution of dimensionless turbulence kinetic energy, 

K/K , on top of semi ellipses with different z and 

D/H ratios 106 

6.22 Streamline, if, patterns over a semi ell ipse for D/H = 1.02, 

z = 0.001, and Re = 100 110 

6.23 Vorticity, ai, patterns over a semi ell ipse for D/H = 1.02, 

z = 0,001, and Re = 100 112 



A.l Schematic diagram of the relationship between general 
orthogonal curvilinear coordinate system and 
Cartesian coordinate system 120 



vm 



LIST OF SYMBOLS 

a One-half the distance between two focus points of 

semiellipses. Figure 5.2(b) 
a,b,c Coefficient functions in the general elliptic equation. 

Equation 5.1 
C »Cp,jCr»Crj Constant coefficients in Equations 3.9 and 3.14 

U U h K 

C-, ,Cp Constants in Equation 5.17 

d Source term in the general elliptic equation. Equation 5.1 

D Width of semiellipses in the longitudinal direction, 

Figure 5.2(a) 
D|,,D. Diffusion rate of the turbulence kinetic energy and 

turbulence lengtn scale, respectively 
e Unit vector 



F Scalar function in Equation 3.3 

g Body force 

h Metric coefficient of the general orthogonal curvilinear 

coordinate system 

H Height of semiellipses. Figure 5.2(a) 

I Vector of inertia force per unit volume 

K Turbulence kinetic energy 

L Turbulence length scale 

p' Fluctuation component of fluid pressure 

P Fluid pressure 

Pj,,P| Production rate of turbulence kinetic energy and turbu- 
lence length scale, respectively 



IX 



r Radius of curvature. Figure 4.2 

Re Reynolds number, (pHVu)/y 

s Distance between two points on curvilinear coordinates. 

Figure 4.1 
S. Source term of the transport equation of the variable (j) 

t Time 

T Vector of surface stresses 

T^ Turbulence factor. Equation 3.21 

u Coordinate of the general orthogonal coordinates system, 

Figure 4.1 

V Fluctuation component of flow velocity 

V Mean flow velocity 

Vl[ Undisturbed upstream flow velocity at body height H 

Vq Mean flow velocity at a reference altitude z^ 

V* Surface friction velocity 

x,y,z Longitudinal, lateral and vertical coordinates, 

respectively 

z Scale of surface roughness 

^ 

Zn Reference altitude 

a, 3 Angles between the axis of symmetry and the radii of 
curvature, r-, and r^, respectively. Figure 4.2 

V Angle between a constant u^ plane and a reference x-z 

plane. Figure 4.2 
ej,,^. Dissipation rate of turbulence kinetic energy and 

turbulence length scale, respectively 
e Angle between the unit vectors e-, and e^. Appendix, 

Figure A. 1(a) 



K Von Karman constant 

ij Laminar viscosity of the fluid 

y Effective viscosity of the fluid 

V Kinematic viscosity of the fluid 

v^ Kinematic eddy viscosity of the fluid 

p Fluid density 

a The r.m.s. value of the fluctuation component of the 

flow velocity 

0^,0. Schmidt number in Equations 3.7 and 3.11, respectively 

T Shear stress 

cj) Dependent variable in Equation 5.1 

\p Stream function 

w Vorticity 

Superscripts 

Denotes dimensionless quantity 
n Index of numerical iteration 

Subscripts 

Ij2j3 Pertaining to directions 1, 2 and 3 in the curvilinear 

coordinate system. Figure 4.1 
I, J Finite difference indices. Figure 5.2(a) 
n,s,e,w Pertaining to the nodes shown in Figure 5.1 
ne,se,nw,sw Pertaining to the nodes shown in Figure 5.1 
N,S,E,W Pertaining to neighboring nodes which lie, respectively, 

north, south, east and west of node p. Figure 5.1 
NE,NW,SE,SW Pertaining to the nodes shown in Figure 5.1 

xi 



NW Pertaining to a node which lies adjacent to a solid wall 

Denotes undisturbed upstream flow condition 

P Pertaining to a central numerical node at which the 

dependent variable is being calculated currently 
t Turbulent quantity 

W Pertaining to a numerical node on the solid wall 

00 Denotes undisturbed upstream flow condition at an 

altitude z = 10 H 



Xll 



1.0 INTRODUCTION 

Analysis of atmospheric motion over obstacles on plane surfaces 
is carried out in the present study to compute simulated wind fields 
over terrain features. Emphasis is on a semielliptical , two-dimensional 
geometry. Numerical simulation of flow over rectangular geometries is 
also discussed. The numerical procedure utilizes a two-equation turbu- 
lence model and the selection of the necessary constant coefficients in 
the model is considered an important part of the present study. 

Although the terrain features analyzed are highly idealized 
geometries which are not exactly duplicated by natural terrain, the 
results present are sufficiently representative of the real world to 
provide insight into the various hazards which are associated with 
aircraft flight in the vicinity of both natural and man-made objects. 
Controlled flight over terrain, such as landing off the runway environ- 
ment or flying over rising terrain, continues to be a safety issue. 
In recent accidents where terrain was a major factor, the air traffic 
control system on the aircraft was equipped with the capabilities 
that should have prevented the accident. This suggests that there 
were certain unknown factors which may include complex and unexpected 
wind patterns which for this safety issue need re-examination. 
Moreover, operation of low-speed helicopters and V/STOL aircraft in 
large metropolitan areas is also affected by zones of recirculating 
wind fields, regions of large velocity fluctuations, and fields of 
vorticity induced by the atmospheric motion near buildings. 

1 



These complex wind fields can be extremely hazardous during landing and 
takeoff. 

There are many other technological developments for which 
analytical solutions of flow over terrain features have practical appli- 
cations. The dispersion of air pollution is significantly influenced by 
atmospheric motion. In complex terrain, the movement and deposition of 
the pollutant can be significantly altered by a hill or building in the 
vicinity of the stack release. The siting of wind turbine generators 
requires knowledge of regions of high wind speed and avoidance of 
regions where recirculation or highly turbulent flow may occur. Addi- 
tionally, the wind load on structures and wind erosion can be increased 
or decreased by the wind characteristics in the vicinity of surface 
irregularities. These fields of technology also require an understanding 
of the atmospheric flow properties about surface obstacles. 

It is well known that surface wind motion is markedly affected 
by the abrupt change of terrain features [1]. The phenomena of flow 
acceleration, of flow separation, and of augmented turbulence in the 
wind field are observed effects. The effect of shear in the approaching 
flow and the strong pressure gradients near the barriers create zones 
of recirculation both upstream and downstream in the vicinity of the 
obstacle. High turbulence is generated in the shear layer bordering the 
recirculation regions due to the interaction between the turbulence 
stresses and the deformation of the mean flow field. Momentum in the 
shear layer diffuses into the wake and into the quasi-potential flow 



Numbers in brackets refer to similarly numbered references in 



the Bibliography. 



outside the wake, setting the wake fluid into motion and smoothing out 
the sharp velocity discontinuities. Downstream of the obstacles, the 
diffusion of momentum gradually thickens the shear layer until the inner 
flow is blended with the outer flow, forming a new boundary layer. 

The nature of these flow fields has been studied experimentally 
as well as analytically. Reported results of some of these studies are 
reviewed in Section 2.0. Due to the complexity of the flow characteris- 
tics near the obstacles, it is necessary to solve the complete equations 
of motion if a realistic recirculation flow field is to be computed. 
Practical methods, however, for solving such equations require the use 
of numerical modeling techniques. 

The numerical procedure employed in this paper is an elliptical 
finite-difference algorithm initially developed by Gosman, et al. [2]. 
The Navier-Stokes equations have been expressed in terms of stream func- 
tion and vorticity. In order to close the governing equations, the 
Reynolds shear stresses appearing in the time-averaged equations of 
motion are expressed by the product of a scalar eddy viscosity and a 
strain rate. The eddy viscosity is expressed as a function of the turbu- 
lence kinetic energy and the turbulence length scale. These quantities 
are determined, respectively, from two corresponding transport equations 
(two-equation turbulence model [3]). The computation of turbulent flows 
based on current turbulence models in all cases requires the determina- 
tion of one or more constant coefficients. One goal of this study is to 
derive the relationships of these constants from the knowledge of turbu- 
lence characteristics available at the present time. Details of the 
turbulence model employed in this study are discussed in Section 3.0. 

3 



Solving partial differential equations by numerical techniques 
requires accurate numerical presentation of boundary conditions. When 
the boundary is coincident with a coordinate line, the finite-difference 
expression at and adjacent to the boundary can be applied at grid 
points constructed on the coordinates without the need of interpolation 
between grid points. In this study, the governing differential equa- 
tions are written in an orthogonal curvilinear coordinate system. The 
curved coordinate surface can then be fit to most terrain shapes of 
interest. Although this increases the size of the governing equations 
slightly, little additional complication is introduced into the calcula- 
tion procedure. The governing equations as well as the coordinate system 
used in this study are described in Section 4.0. The numerical proce- 
dures and the imposed boundary conditions are given in Section 5.0. 

The computation has been carried out for atmospheric flow over 
two-dimensional semiell ipses, due to the resemblance of their geometries 
to the shapes of hills. An analysis of the effect of the geometry of 
the obstacles on the wind field is carried out in the present study by 
varying the aspect ratio of the semiell ipse. The atmospheric flow is 
considered neutrally stable. The concept of a neutrally stable atmo- 
sphere is meaningful in strong wind conditions [1,4], which is of primary 
importance in this study. Good agreement between the numerical predic- 
tions and the experimental data obtained from both atmospheric boundary 
layer wind tunnels [5,6,7] and full-scale field measurements [8,9] is 
shown in Section 6.0. The comparisons are carried out for flows over 
solid fences, rectangular blocks and escarpments, in addition to the 
case of flow over semiell iptical bodies. 

4 



Parameters influencing the flow field in the vicinity of surface 
obstacles are also discussed in Section 6.0. It is found that the 
primary flow parameters affecting the flow field are surface roughness 
scale, the shape of the obstacles as well as the undisturbed upstream 
turbulence level. Increasing the scale length of the surface roughness 
results in reduced wind speed in the lower layer of the atmosphere due 
to the increase in surface frictional drag. Thus, the streamline is dis- 
placed to a greater height when the flow passes over the obstacle with 
larger surface roughness. For the flow conditions investigated in the 
present study, flow separation is not obtained for flow over semiellipses 
when the molecular Reynolds number. Re, is larger than 10 . Conditions 
with Re = 100, however, show that separated flow is induced by semi- 
elliptical obstacles. The mechanism of flow separation induced by a 
semiellipse is the same as in the case of flow over a gradually slop- 
ing surface for which the flow separation is caused by the interaction 
between the viscous force, the pressure force, and the turbulence level. 
For flow over bluff bodies, e.g., solid fences and rectangular blocks, 
a large downstream recirculation bubble is created due to the inability 
of the flow to negotiate with the abrupt change of the surface shape. 
It is found that the aspect ratio of the block and the turbulence level 
of the approaching flow cause significant effect on the size of the 
downstream recirculating flow region. Increasing the aspect ratio 
and/or increasing the turbulence level result in flow reattachment 
close behind the obstacle. Details of flow separation for semi- 
ellipses as well as for bluff bodies are discussed in Section 6.5. 



2.0 REVIEW OF PREVIOUS WORKS 

The purpose of this section is to review the existing theoret- 
ical analyses and experimental results concerning the aerodynamics of 
surface obstructions in turbulent flow fields. Frost [10] has compiled 
an extensive survey of flow properties around man-made obstacles to the 
wind. The obstacle exerts a drag force on the wind field and distorts 
the flow approaching and passing over it. Recirculation regions are 
created due to the strong pressure gradients near the obstacle. The 
separation and reattachment points in two-dimensional laminar flow fields 
correspond to the location where the normal velocity gradient is zero; 
thus the wall shear is zero. In a turbulent flow field, this is also 
assumed to be true in terms of the ensemble time-averaged quantities [11] 

Various approaches have been used to theoretically describe 
certain features of turbulent flow over barriers. Tani [12] used a 
diffusion equation to investigate the velocity distributions behind a 
permeable fence. Velocities in the vicinity of the fence were not 
accurately predicted. Plate [13] and Chang [14] indicate that the 
velocity distributions in the shear layer bordering the recirculation 
bubble behind an obstacle agree with prediction from mixing-layer 
theory. No flow separation, however, is taken into account. The 
agreement may be due to the arbitrary adjustment of coefficients used 
in the theory. 

Counihan , et al . [15] apply a small perturbation technique to 
study flow over a low block. Their conclusion that similarity prevails 



and that the velocity deficit behind the block correlates with a single 
curve is also supported by their experimental work. Some evidence [5], 
however, is given that such similarity in general may not exist. It is 
unlikely that details of the flow properties in the recirculation 
regions near the barrier can be accurately predicted by using small 
perturbation techniques. 

Several investigators, e.g., Taylor [16], Jackson and Hunt [17], 
and Deaves [18], are studying turbulent flow over gently rolling, low 
hills. The application of their solution to large or steep hills must 
be treated with caution. Frost, et al . [19] extended a turbulent boundary- 
layer theory to approximate the atmospheric motion over semi elliptical 
hills. The Reynolds stresses are assumed to be important only close to 
the ground, and the pressure distributions driving the boundary layer 
are estimated from in viscid flow theory. This assumption is supported 
by the work of Jackson and Hunt [17] for flow over low hills. 

There has been a fair amount of work recently on the problem 
of flow over circular or elliptical cylinders, most of which assumes low 
Reynolds number flow, e.g., Lin and Lee [20] and Haussling [21], or low 
turbulence high Reynolds number flow, e.g., Guven, et al . [22]. None of 
these studies are easily extended to the analysis of turbulent flow 
over barriers. 

At present, most understanding of turbulent flow about surface 
obstructions is obtained from empirical observations. Frost and Shi eh 
[1] have made an extensive survey of the available literature. The 
following summarizes their survey. 



The presence of a bluff body in the flow field retards the 
approaching wind and creates a separation bubble in front of the 
obstacle. The pressure on the front face of the barrier is strong due 
to the conservation of energy. Outside the upstream separation bubble 
these strong pressure gradients displace the flow and accelerate it 
as it passes over the obstacle. A lee separation bubble is created 
if the flow cannot negotiate the change in slope of the barrier. Mass 
flow into and out of the recirculation zone is conserved. At the 
reattachment point the flow splits, part flowing into the bubble and 
part flowing downstream. An equal amount of mass is in turn entrained 
into the free stream from the bubble through the turbulent transport 
processes in the shear layer bordering the recirculation zone [23]. 
Downstream of the reattachment point, a new boundary layer develops. 
The retarded wind profile readjusts to the local surface conditions, 
thus recovering to a new, fully developed profile at a distance far 
downstream. The behavior of this profile has been studied with 
boundary-layer theory [24]. 

The flow properties near an obstacle depend on the geometry of 
the obstacle, the local surface conditions, and the nature of the 
upstream wind. Behind the obstacle, a closed recirculating flow region 
(recirculation bubble) bounded by a stagnation streamline is generally 
accepted as a description of the separation zone in two-dimensional 
flow. The concept of a closed recirculation bubble can also be 
applied to the case of thick boundary layers, e.g., atmospheric 
boundary layer, if the turbulent transport process is more dominant 
than the mean convective transport process [25]. Strong turbulence 

8 



in the undisturbed flow causes a smaller recirculation zone on the lee 
side of most obstacles [26]. 

One of the characteristics of the wake generated by an obstacle 
is the velocity deficit in the wake. The decay rate of the mean 
velocity deficit is independent of the obstacle size, but the magnitude 
of the disturbance is greater for larger obstacles; hence the wake 
persists farther downstream [11]. 

The shape of a barrier has a definite effect on the lee side 
velocity distributions. A solid vertical sharp-edged fence generally 
creates a stronger disturbance than a porous fence, although the dis- 
turbance may be carried farther downstream behind a porous fence. 
Acceleration of the flow occurs in a region located a distance two to 
three fence heights downstream and slightly above the lee side recir- 
culation bubble. Also, large reductions in velocity occur close to the 
ground immediately behind the fence for solid vertical plate fences. 
Increasing the permeability of the fence causes less reduction of the 
flow behind the fence as well as weaker wind speed in the region of 
flow acceleration over the bubble. The area of reduced wind speed 
behind a porous fence persists greater distances downstream due to more 
momentum being convected into the leeward area through the fence 
openings. The windward flow patterns, however, are insignificantly 
affected by the barrier porosity. 

The maximum wind acceleration expressed in terms of the ratio 
of the local wind speed relative to the undisturbed value measured at 
the same height above the local surface occurs at the crest of a 
surface feature in the absence of flow separation [7]. When flow 

9 



separation occurs, the location of the maximum speed lies a small dis- 
tance downstream of the barrier and slightly above the separation region, 
depending on the barrier shape, the local surface conditions and the 
turbulence of the approaching flow. 

VJhen the wind approaches a three-dimensional feature, flow 
acceleration occurs not only on top of the feature but also in regions 
near the sides. In a stably stratified atmosphere, the wind tends to 
pass by an obstacle rather than passing over it. Due to this effect, 
the velocity on top of the obstacle is significantly reduced for flow 
over terrain features with small aspect ratios. However, for flow over 
a long ridge the velocity at the center portion of the ridge is not 
significantly affected by atmospheric stratifications due to quasi - 
two-dimensional flow conditions generally prevalent at that portion 
[7]. Additionally, the flow separation point is observed to move down- 
wind under stable conditions [1]. 

High turbulence generally occurs in a shear layer originating 
from the sharp edge of a bluff body. This high turbulence diffuses 
vertically as flow is convected downstream, thus thickening the shear 
layer. The velocity fluctuation in the streamwise direction is not as 
sensitive to the barrier shapes as is the mean velocity [27]. High 
turbulence intensities, because of their definition, occur in regions 
of low wind speed, normally lying close to the rear face of an obstacle. 
Absolute turbulence, however, is generally higher in the regions of high 
wind speed. In the center of a long fence, where quasi-two-dimensional 
conditions are achieved, the streamwise fluctuation decreases with 
the increasing angle between the wind direction and the fence. 

10 



In general, the magnitudes of the vertical fluctuations are as large 
as the lateral fluctuations [27]. 



n 



3.0 TURBULENCE MODELS 

3.1 Mean Flow Equations 

In the study of atmospheric flow over surface obstacles, the 
flow associated with separation and reattachment is so complicated that 
the general equations of fluid motion, the Navier-Stokes equations, 
must be solved if realistic results are to be obtained. Atmospheric 
flows are generally treated as incompressible turbulent flow, i.e., the 
fluctuation of fluid density is negligible in comparison with the mean 
density in the statistical steady state. This is especially true when 
the atmosphere is neutrally stable. Based on the assumption of 
incompressible flow, the mean velocity and pressure are governed by 
the equations: 

Continuity equation: 



9x.. 



= 



(3.1) 



Momentum equation: 



8V. 9V. 
1 + V ^ - 



SP 



dt J 3x. 



p 9x 



3X. 



V 



3V. 
Bx" 



- V .V . 
1 J 



+ g. 



(3.2) 



where q. is the external force acting on the fluid element in the i 
direction and v is the kinematic viscosity of the fluid. 



12 



For a stable or an unstable atmosphere where thermal effects are sig- 
nificant, the influence of density variations due to temperature changes 
is often approximated with the Boussinesq approximation. The density 
effect thus appears only in the force term of Equation 3.2. The 
temperature distribution is then governed by a scalar transport equation 
for turbulent flow [28]: 

Scalar equation: 



3F 
9t 



^^^ = l^l:<-p^■f)^^F <'•'> 



where F is a scalar function, e.g., temperature; f is the fluctuation 
component of F; and Sr is a volumetric source term. Since the effect 
of heat transfer is negligible if the wind is strong [4], which is 
assumed in this study, the thermal effects on the flow are neglected, 
and therefore, the governing equations in this study are Equation 3al 
3o2. The external force term of Equation 3.2 is neglected due to its 
small effect compared with the shear stresses. 

3.2 Closure Problem 

The main problem in solving Equations 3.1 and 3.2 and thus cal- 
culating turbulent flows is the determination of the Reynolds stresses, 
-pv.v.. One can derive an exact transport equation for the Reynolds 
stresses from the Navier-Stokes equations. This derivation, however, 
simply introduces a third-order correlation term, v.v.v. , [28], which 

still does not permit the system of equations to be closed. Continuing 

13 



the process by deriving exact equations for higher order velocity 
correlations results in still higher order correlations, and thus the 
set of equations can be closed. Thus, at some point a turbulence model 
must be introduced which approximates the turbulent correlations at 
some order in the hierarchy as functions of the lower order correlations 
and/or mean flow quantities. Turbulence is thus simulated by a turbu- 
lence model ranging from a simple algebraic model of mixing length to a 
more complicated differential transport equation model for the higher 
order correlations [3]. The additional equations, differential and/or 
algebraic equations, used for turbulence modeling together with Equations 
3,1 and 3.2 form a closed set of equations. The turbulence model used 
in this study is described in the following. 

3.3 Concept of Eddy Viscosity 

The eddy viscosity concept of Boussinesq [23] is adapted to 
express the Reynolds stresses as: 



■pv.v. = pv^ 



9V. 8V. 
— - + — ^ 



3x. 8x^ 



-|k6.. (3.4) 



where 6.. is the Kronecker delta (6.. = 1 if i = j, otherwise 6.. = 0), 



and K is the kinetic energy of turbulence, l/2(v.v.). The eddy vis- 
cosity, Vj., used in Equation 3.4 is based on the assumption of an 
analogy between momentum transport by turbulence and momentum transport 
by molecular motion. Although this concept has frequently been criti- 
cized as physically unsound, it is found to work well in several flows 
[29]. 

14 



The concept of eddy viscosity has been used in isotropic turbu- 
lent flows by assuming a scalar quantity for v. . For non isotropic 
turbulence fields, different eddy viscosities are sometimes used for 
turbulent transport of momentum in different directions [29]. In certain 
flows, e.g., wall jets and asymmetric wall shear layers, negative 
values of v^^ occur where the Reynolds stresses have the opposite sign 
to that of the velocity gradients (see Equation 3.4). This, of course, 
is not physically meaningful. In most cases, however, the regions with 
negative eddy viscosity are small, and therefore the error induced in 
those areas is of little practical importance [3]. The application of 
Equation 3.4, however, simply shifts the requirement of Reynolds 
stress modeling to the problem of eddy viscosity modeling. 

To overcome this problem, v^ is expressed as being proportional 
to a velocity scale, /K> and a length scale, L, characterizing the 
turbulent motions. Prandtl and Kolmogorov [3] propose: 

V. = C /K L (3.5) 

t y ^ ' 

where C is a constant. The kinetic energy of turbulence, K, and its 
length scale, L, are computed from transport equations. This model, 
called a two-equation model, is employed in this study and is discussed 
in the next section. 



3 . 4 T he Two-Equation Mod els 

Launder and Spalding [3] Indicate that a variable A constructed 
by any combination of K and L, i.e., A = k'^l'", where m and n are real 



15 



numbers, can be used in the transport equations for the two-equation 
models. Rodi and Karlsruhe [29], however, suggests that a transport 
equation for L should not be used due to the difficulties associated 
with deriving a generally valid formula. References [3,30,31] counter 
this and indicate that a common form for the L transport equation does 
exist. A transport equation for L having a general form like Equation 
3.3 has been used in References [2,11,32,33]. Evidence that realistic 
results can be achieved by using a K-L model for complicated separating 
and reattaching flow phenomena is reported in [11,32]. One subject of 
this study is to demonstrate the validity of the K-L model for calcul- 
ating turbulent recirculating flows. 

Forms of the two-equation model other than that used here are 
reported [3]. A common characteristic of all two-equation models is 
that a set of constant coefficients which appears in the transport equa- 
tions must be determined empirically or by inductive reasoning. One 
semi -empirical procedure for determining the constants assumed the 
turbulence is always in local equilibrium [3,32]. This approach, how- 
ever, can break down in regions where the turbulence production rate 
dominates the dissipation rate. This study proposes an alternate 
method to determine these constants based on presently available 
knowledge of turbulent motions. This method is discussed in detail 
in Section 3.6. 

3.5 K-L Two-Equation Model 

Derivations for the transport equation for K and for L are given 
in Hinze [28] and Mellor and Herring [30]. A brief description of their 

16 



development is given here. In the flow field, the convective and diffu- 
sive transport of turbulence, i.e., the history of the turbulence, is 
accounted for by solving two differential transport equations. These 
are the transport equations for turbulence velocity scale and for 
turbulence length scale. The physically most meaningful velocity scale 
is /K, where the turbulence kinetic energy is defined as half of the 
sum of the mean square values of the velocity fluctuation components. 



(l/2)v.v^. The transport equation for K can be derived from the Navier- 
Stokes equations [28]: 



3K , y 9K _ 
8t i ax.. 



9x. 



1 



■,^- 



9V. 



- ^"j 33^^^B3r: 




II 



III 




IV 



I == D|, = Diffusion 



II = p E Production 
K 



III = Viscous transport 
IV = Gj, = Dissipation 



(3.6) 



The viscous transport term in Fquation 3.6 can be written as 



"s^'i 



3V. 



3v. 

3xj ax^ 



= V 



r 2 



8X.2 ^" ^x.^xj'i'j 



17 



since the terms 





and 


9^v. 

V. ^ 
1 3X,.3X. 


3v^ 3V. 

9x. 3X- 



are zero for incompressible fluid. Since the viscous term is generally 
small except when spatial gradients of the Reynolds stress components 
are extremely large, it is generally neglected. The diffusion term in 
Equation 3.6 resulting from the interacting of velocity and pressure 
fluctuations is assumed to be subject to gradient diffusion. The 
diffusion term is then written 



D. = ^ 



K dx. 



t 8K 



^K ^^* 



(3.7) 



where the Schmidt number a^ is a constant which is discussed in 
Section 3.6. Since the energy cascade process is independent of 
the molecular viscosity except at the final stage where the turbulence 
kinetic energy is converted into heat by viscous dissipation, the 
dissipated turbulence energy is modeled by 

e,< = -CpK^^V"' (3.8) 

where Cr, is a constant which is discussed in Section 3.6 . Equation 
3.8 involves Kolmogorov's assumption [31] that the dissipation rate 
of K is determined by the energy-containing motion at high Reynolds 
number. The length scale, L, characterizes the size of the large, 
energy-containing eddies of the turbulence field. 

18 



The K equation based on the above assumptions is thus modeled as 



oK , w oK 
9t ^j dx. 



9X 



V 



t 3K 



J i 



\ ^"jj 



r9v, 



+ V, 



9V_. 



ax. 



8x. 



9x_. 



,3/2 



- C 



D L 



(3.9) 



D 



K 



K 



K 



The length scale characterizing the size of energy-containing 
eddies, in turn, is assumed to be subject to transport processes and is 
also determined from a differential transport equation. As mentioned 
in the previous section, the L equation has the general form of equation 
3.3. The rate of change of L is balanced by the convective transport, 
the diffusive transport, D. , and the production and dissipation rate of 
L in the energy cascade process, i.e., P. and e. , respectively, hence 



St J 9x, ^L ^L ^L 



(3.10) 



It is assumed that the diffusive transport of L is also gradient-driven 
The D. term is therefore modeled as 



D. = 



9x. 



t 9L 



a, 3x 



(3.11) 



where a^ is a constant which will be discussed in Section 3.6. 

Considering that increasing the turbulence kinetic energy dissipation 

increases the fatality rate for small eddies and thus effectively 

increases the eddy size, the production term for L, P. , should be 

related to the dissipation of turbulence kinetic and some characteristic 

time scale. Hence, it is argued that 

19 



Pl " ^K K 



which reduces to 



P^_ ^ C^/K 



(3.12) 



where Cr is a constant which will be discussed in Section 3,6. 
The vortex stretching connected with the energy cascade reduces the 
eddy size; thus, the rate of dissipation of L, e, , is subject 
to the interaction of the mean motion and the turbulence fluctua- 
tions. 



^L = - ^R 



9V. 
't 9x; 



L K 



-1 



(3.13) 



where Cn is a constant which will be discussed in Section 3.6. 
With Equation 3.5 and Equations 3.11 through 3.13, Equation 3.10 
becomes 



9L 



3t ^j 8x 



w: lf§: ^^E^- Vr- 

J L J 



/R" 



3V 



3V 



L^^j 



9x 



ij 



9V. 
9x" 



D. 



(3.14) 



The length scale equation. Equation 3.14, traces the historical develop- 
ment of the length scale at a given location and therefore is strongly 
dependent on the initial value of L. 

Both Equations 3.9 and 3.14 involve several constants which must 
be determined before these equations can be solved. The determination 

of these constants is discussed in the next section. 

20 



3.6 Constants of the K-L Model 

The K-L model described 1n the previous section contains six 
coefficients. They are aj., a,, C , C^, C^ and C^^. Conceptually, these 
coefficients are functions of dimensionless parameters, e.g., Reynolds 
number, which govern the motion of the flow [2]. In high Reynolds 
number flows, however, the coefficients are generally treated as con- 
stants. For fully developed turbulent flows, the Schmidt numbers o^ 
and a, used in Equations 3.9 and 3.14, respectively, are generally 
assumed to be of the order of unity [2]. 



a^=1.0 

a^=^.0 (3.15) 



Indirect support of this assumption is given by Reynolds [34], who 
shows that the ratio of the eddy viscosity to the eddy diffusivities 
of enthalpy and of other properties such as fluid contaminants is close 
to unity for shear flows. 

In the analysis of turbulent flows using two-equations models, 
determination of the constants C , C^, C^ and C„ depends either on 
experiment and/or on an exact analysis of some simple flow. Launder and 
Spalding [3] proposed a wall region approach in which the turbulence 
convection and diffusion are negligible. They analytically determine 
the constant coefficients used in their two-equation model in this 
region. Frost, et al . [32] also used this approach in determining the 
coefficient used in their K-L model. Moreover, the Prandtl mixing 



21 



length, L = kZq, was specified as a boundary condition for the transport 
equation of length scale. Equation 3.14. Realistic results are obtained 
by applying this model to calculate atmospheric flow over forward-facing 
steps [32]. For flow over rectangular blocks, the wall region approach 
gives unrealistic results in both mean wind field and turbulence field 
[11]. Shi eh, et al. [11] investigated the influence of the coefficients 
C , Crj, Cr and C^, on the flow field computed and then proposed 

C = 0.416, Cn = 0.416, Cn = 1.44 and Cr = 0.6 for calculating flows 

H u R t ^ 

over a rectangular block. Realistic results in both mean wind field 
and turbulence field were obtained [11 ,26]. However, validation of 
the coefficients proposed by Shi eh, et al . [11] in other flow situa- 
tions requires further investigation. 

In the present study, the constant coefficients C , Cp, C^ and 
Cn described in the following are based on an understanding of the 
properties of turbulence as influenced by their transport process. 
Thus the coefficients are expected to have general applicability. The 
coefficients developed in this study are successfully applied to the 
cases of flow over semielliptical obstacles, fences, rectangular blocks 
and escarpments. Good agreement between the numerical solutions 
and the empirical data of both wind tunnel and full-scale field measure- 
ments are obtained. Comparisons of experiments and analyses are shown 
in Section 6.0. The coefficients are developed as follows. 

Taking the ratio of the production terms to the dissipation 
term in Equations 3.9 and 3.14, one obtains 



(3.16) 
22 



Pi/ 


f^' 1 


,? 


r9V. 9v.n 


9V. 


K _ 


U 


L 


T , J 


1 


^K 


N 


K 


9x. 9x. 
_ J '_ 


3X. 






['eI 


*/ 


l^M 



8X. 



ax. 



3X. 



(3.17) 



To establish some physical explanation of the relationship between the 
turbulence kinetic energy and the turbulence length scale, one can 
argue that the growth of the turbulence length scale Is a result of the 
dissipation of energy, which creates the biggest fatality rate for 
small eddies. Thus, when dissipation exceeds production, the length 
scale, which is taken to represent some effective mean of the collection 
of eddy sizes, will increase with the depletion of the smaller eddies. 
The reduction In length scale. In turn, may be thought of as being 
caused by the tendency of shear stress to rupture the larger eddies. 
Therefore, when the production of turbulence kinetic energy exceeds 
dissipation, the length scale is expected to decrease, and when dissipa- 
tion of turbulence kinetic energy exceeds production the length scale 
Is expected to grow. With this In mind it is proposed that 



— °^ -^ — 



which gives the relationship 



^L Vr 



= 1.0 



(3.18) 



This assumption is obviously valid for the local equilibrium case, I.e., 
P|. = Ew and P. E e. . However, it is not restricted to only equilibrium 
flows. 



23 



Due to the interaction between the mean motion and the turbulent 
motion, energy is extracted from the mean flow and converted into heat. 
The local production of turbulence energy, however, is not always 
balanced by the local viscous dissipation; thus the turbulence does not 
necessarily gain energy from the mean flow at the same rate as it loses 
energy through viscous dissipation. 

Conceptually, the ratio (C /Cq) appearing in Equation 3.16 can 
be determined from the ratio Piy/eiy and the local conditions of the 
turbulence and of the mean flow. Rodi [29] proposed an equation for 
evaluating the effect of the Piz/eiy ratio on the coefficients used in 
a K-£|, model. His coefficients vary significantly with the distortion 
rate of turbulence. However, an asymptotic constant value of C is 
achieved if Piy/siy > 1-0. The distributions of K and L are shown [11] 
to be strongly dependent on the ratios C /Cp, and C^/C Cn, respectively, 
rather than on the magnitude of the constants themselves. However, 
close behind a rectangular block, x = 1.0 H, where the mean flow is 
highly distorted, the distributions of K and L are rather insensitive 
to the value of these ratios (Figure 3.1a). On the other hand, in the 
relaxation flow region far downstream from the block, x = 20 H, the 
dependence of K and L on these ratios is seen to be much stronger 
(Figure 3.1b). 

Considering first the region directly behind the block where 
the computation is insensitive to the value of C /Cr. and Cr/C Cp,, it 
is reasonable to assume that the ratios C /Cr. and Cr/C Cn approach 
unity. Therefore, 



24 



r c 



• ■ 



^ = 0.416 

'D 

= 2.4 



C Cd 



7^= 1.0 



= 1.0 



C Cp 



z/H 




1 K/Kq 2 



(a) Profile at 1.0 H behind a rectangular block 



z/H 




1 K/K^ 2 



2 L/L. 



(b) Profile at 20.0 H behind a rectangular block 



Figure 3.1 Effect of the ratio of constants on the distribution of 
K/Kq and L/Lq behind a rectangular block [11] 



25 



C == Cn (3.19) 

y D 

Cp = C Cd (3.20) 

E y R 



Equations 3.19 and 3.20, of course, are consistent with Equation 3.18. 
Under the above assumptions, the ratio of turbulence production 
rate to turbulence dissipation rate is thus determined solely by the 
local turbulence properties, K and L, and the mean flow deformation; 
see Equations 3.16 and 3.17. 

As the turbulence reestablishes its fully developed state, the 
large eddies will break down through intertial interaction and transfer 
energy to the smaller eddies [28]. As the eddies become smaller, 
energy dissipation due to viscous effects becomes more and more impor- 
tant. Hinze [28] indicates that the turbulence should be characterized 
by the relative rate of change in turbulence energy (DK/Dt)/K. Because 
the rate of change in eddy size DL/Dt is dependent on the local turbu- 
lence properties, this rate of change can be assumed to reflect the 
condition of the local turbulence level, K, A characteristic parameter 
of turbulence can, therefore, be defined as: 

^ (DL/Dt)"^ 

It is noteworthy that the dimension of T^. appears to be the same as that 
of frequency, (time)" . An eddy frequency associated with vorticity 
fluctuation has long been recognized as an important parameter of 
turbulence [3]. 



26 



Now consider the range where the turbulence is characterized by 
P|, and e, . The factor 1^ can be expressed as 



K 



f 2 



1 




_L_ 


7 


"av. 9V." 
— - + — ^ 

3x. 9x. 


S'k 



3X. 



(3.22) 



Based on the assumption for deriving Equations 3.19 and 3.20, the ratio 
of coefficients is again assumed to approach unity, i.e., 



C Cd = 1.0 



(3.23) 



From Equations 3.18, 3.19, 3.20 and 3.23, the relationship between the 
constant coefficients can be written as follows: 



Cn = C 

D y 



Cn = 1/v^ 

R y 



C,- = v^ 



(3.24) 



It is interesting to note that the value of T^ for smaller eddies where 



£,, and P. are dominant becomes 



T^ =« 



2 L 



(3.25) 



which becomes /K/L when Equation 3.24 is introduced. This is consistent 



27 



with the assumption throughout this section that the relative turbu- 
lence properties; Equations 3.16, 3.17 and 3.22; are insensitive to 
the ratios of constant coefficients in regions of high rate of turbu- 
lence production. Equations 3.24 and 3.25 imply that this, assumption 
is also valid in regions of high rate of turbulence dissipation. 

The value of C can be determined from the concept of a constant 
shear layer where the length scale, L, is proportional to the distance 
from the solid wall, i.e., a mixing length hypothesis [11]. 



C = /FTpK = V^/ZK (3.26) 



where x is the wall shear stress and V^ is the friction velocity. 
Equation 3.26 indicates that the value of C depends on the flow being 
investigated. 

Gosman, et al . [2] propose that the value of C can be taken as 
an asymptotic universal constant in high turbulence flows where the 
turbulence Reynolds number R. = /KL/v is large. The value of C 
determined from Equation 3.26 is taken as the asymptotic value when 
applied to the high turbulence level flows considered in this study. 

For atmospheric flows, C = 0.416 is proposed [11]. Substituting 
this value into Equation 3.24, the following constants result. 



Cn = C = 0.416 
D y 



Cf^ = 1.55 



C^ = 0.645 (3.27) 

28 



These values are comparable to the previous work [11], where 
C = Cn = 0.416, Cn = 1.44 and Cc = 0.60 where used. 

y L) K t 

Launder and Spalding [3] indicate that the ratio V^ /K lies 
between 0.25 and 0.3 in various flows. A wider range, however, is 
reported by Zeman and Tennekes [35] who give 0.226 to 0.34. Field data 

for atmospheric boundary layers generally show smaller values of 

2 

V^ /K as compared with results from laboratory flows. This implies 

that the turbulence generated in the laboratory is smaller than that 

2 
generated in the real world. The ratio V^ /K for an atmospheric sur- 
face layer is considered to be between 0.17 [36] and 0.26 [35]. The 

numerical values of the constant coefficients of Equation 3.27 corre- 

2 
spond to the case of V^ /K = 0.17. No significant difference of 

2 
numerical solutions has been observed when V^ /K = 0.25 is employed. 

Consider the decay processes of turbulence behind grids. 

Launder and Spalding [3] propose that the value of Cn is half of the 

ratio Cnj/C if the K-L model is employed. With the aid of Equation 

3.19, one can easily show Cp is 0.5, which is comparable to the present 

predictions, Equation 3.27. In order to verify the capability of the 

present model for calculating turbulent flows, atmospheric flow 

over surface objects with different geometries, i.e., forward steps, 

rectangular blocks, fences, and semielliptical shaped bodies, is 

investigated. Comparison is made with wind tunnel data and full-scale 

field measurements. Good agreement betv/een the analytical solutions 

and the experimental data is obtained and is shown in Section 6.0. 



29 



4.0 FUNCTIONAL EQUATIONS IN GENERAL ORTHOGONAL 
CURVILINEAR COORDINATES 

The governing equations of motion for the turbulent atmospheric 
boundary layer under neutral stability conditions, Equations 3.1 , 3.2, 
3.9 and 3.14, are formulated in orthogonal curvilinear coordinates in 
this section. Considering a general orthogonal coordinates system; 
u.,,u«,u^, as illustrated in Figure 4.1; a differential line element 
in this system, dt, can be related to the Cartesian coordinates; 
x,y,z; by 

d? = e dx + e dy + e^ dz 

where e , e and e are unit vectors in the x, y and z directions, 
respectively. 



U 




Figure 4.1 An orthogonal curvilinear coordinate system 



30 



The magnitude of the vector ds can be expressed as 



(ds)^ = (h^du^)^ + (h2du2)^ + i^^^^^^^ 



where the metric coefficients, h.'s, are defined by 



^ = 



fax] 

LriJ 


2 

+ 




2 

+ 


1 3z 

3U. 


2' 



nl/2 



i = 1,2,3 



The gradient is defined by 



„ _ ^1 3 I ^2 3 ,^3 3 
- h^ 3u^ hg 3U2 hg 3U3 



where e-, , e^ and e^ are the unit vectors in the directions 1, 2 and 3, 
respectively. 

4.1 Continuity Equation 

The steady state conservation mass law for incompressible flow 
can be written in vector form as: 



V • ^ = 



(4.1) 



In the orthogonal curvilinear coordinates system. Equation 4.1 appears 
as 



V • V = 



1 



^^hS 



8U7(^2h3Vl)^3^(^3^V2)^4^'^'2V3^ 



= (4.2) 



31 



"Inertia' 
_ Forces _ 


+ 


" Body " 
_Forces_ 


■ + 


Surface 
_ Forces _ 



4.2 Conservation of Momentum Equations 

According to the D'Alembert principle, the forces acting on a 
control volume are balanced by an inertia force. The inertia force, 
in turn, represents the inflow of momentum to the control volume 
(Newton Law of Motion). 



= (4.3) 



The terms involving inertia forces and surface forces must be 
examined in the curvilinear coordinate system. In order to simplify 
the analysis, choose the curvilinear coordinates such that the constant 
Uo surfaces are flat planes (Figure 4.2). Therefore, the fluid motion 
in the directions of u, and u^ makes no contribution to the inertia 
forces in the direction of u^. The center of curvature, 0, for a 
'j„ coordinate, described by a position vector f measured from the 
original of the Cartesian coordinates system, lies on a constant u-. 
plane. The equations of flow written in the general orthogonal 
coordinates system are still very complex since the trace of the center, 
(Figure 4.2), may be described by a \/ery complex function in the 
x-y-z space. In order to reduce the complexity in the computational 
procedure, the coordinate system used in this study is taken as 
axisymmetric (Figure 4.2). This implies that the Uo coordinate repre- 
sents the angle of revolution, y, about the axis of symmetry from a 
given reference plane, e.g., the x-z plane. Referring to Figure 4.3, 
the inertia component of the rate of momentum change per unit volume, I, 
after neglecting the higher order terms can be expressed as follows: 

32 




Center- ^^ 

^■-^- 3 ^'a/ie 



^V 




^^^nter of 



ds 



'^tur^ 



e of 



'"""aPraW 



Constant u^ 
Plane 

Symmetry Axis 



Figure 4.2 An axisymmetric coordinate system 



33 



CO 

■1^ 



Center of Curvature of ds 




Symmetric Axis 



Center of Curvature of ds 



Figure 4.3 Velocity vectors of a fluid element shown in a constant 
Uo plane (reference Figure 4.2). 



I = -pen 



n)v 



: — Vn zrr + yn -J+ - ^o sin t> 7j+ 



Dt " '2 ^ "2 dt 



dt 



pe2 



'DV, 



Dt 



^lt^V,f-V3COS3;f 






(4.4) 



where 3 is the angle between the radius of curvature rg and the axis of 
symmetry and a - (71/2) - 3. Substituting 

V = r ^ 
n 1 dt 



V = r ^ 
h ^2 dt 



and 



V = r ^ 
^3 ^3 dt 



into Equation 4.4, the inertia forces per unit volume, I, become 



I = -pe. 



Dt rg r^ r^ 



- pe. 



DV, 



Dt 



f V V 

3 12 

:r- cos 3 + ^r^ 



- pe 



3 Dt 



(4.5) 



35 



Considering the steady state and substituting Equation 4.2 into 
Equation 4.5, the inertia terms in direction 1 are expressed as 



I. = 



br. fe (^2Vi^i) ' anr (^ V2V ' 3^7 (N^^zVi) 



h^h2..3 



— £l-J- + — t^ + 



1 



pV. 



sin B 



(4.6) 



\. J 



The terms in the first bracket on the right-hand side of Equation 4.6 
can be written as 



-p/h 



^•^2^3 



9h^ a 



3h ^'^1 

"13 2 1 3U2 Bu^ 1 2 3 1 r 1 2 3 I du^ 



and thus 



ii = -p 



J-V.(Vh^V^)- 



Vl^ Vi!!!i ¥i!!!i,Vi 



(4.7) 



2 2" 

^3 . ^2 

•^3 2 



The radii r-j and r2 are given by [2] 






(4.8) 



36 



1 ^\ 

h^h2 8u^ 



(4.9) 



Substituting Equation 4.8 into 4.7 and considering the coordinates 
system in which the metric coefficients h-, and h^ are invariant with 
respect to the direction 3 coordinate, Equation 4.7 becomes 



In=-P 



J-v.(?h^v^) 



V^^ 9h^ 

rr 3u7 



sin 3 



(4.10) 



Expressions for the momentum change in directions 2 and 3 are similar 
to that given by Equation 4.10 and are included later in the complete 
equations of motion. 

The surface forces in Equation 4.3 include contributions from 
the fluid pressure and from the shear stresses due both to molecular 
viscosity and to turbulent fluctuations. The contribution of fluid 



pressure to the force per unit volume, ?" , is 



F = -VP 
P 



(4.11) 



The surface stress T can be expressed in general form: 



T = T, + Tj + T3 



where 



Tj = «iL-i 



j = 1 ,2,3 and i = 1,2,3 



37 



The component T.. represents the shear stress acting In the ith direc- 
tion on a surface perpendicular to the jth direction. The components 
making up T.. are expressed as follows [2]: 



^n ' '^e 


[2 3^1 2V2- 
h^ 3U^ r^ 


^22 = ^e 


[2 ^^2 , ^v; 
^2 ^"2 ^2 



T33 = ~ [2(V^ sin 3 + V2 cos 3)] 



^12 ^ ^21 " ^e 



> 3" 2 


['^1 

H 



'2 3 



h, 3u, 



^2 



^23 " """32 " ^e 



"^3 3 
hg 3U2 


['3IJ 



^13 ' "^31 = ^e 



'''3 3 
h^ 3u^ 


f'3]l 

M 



(4.12) 



where y is the effective viscosity defined as 



Pg = y + Vt 



The contributions of the shear stresses to the force on a fluid 

element follow a similar derivation as that for the inertial terms given 

in Equation 4.10. Therefore, the general expression for the momentum 

balance given by Equation 4.3 reduces as follows (see also Reference [2]) 

38 



Di recti on 1 : 



T^v.[h,(7v, -t^)]- 



PV2 - T22 



p^a - ^33 



sin $ 



PV- - T^^ 



8h 



1 ^ 1 3P 



8u^ h^ 8u^ 



= 



(4.13) 



Direction 2: 






V- [h2(VV2 - t2)j - 



pV/ - T 



11 



^h - hs 



cos 3 



PV2 - T22 



8h 



2 ^ 1 ap 



it-^tH- = 



aup h2 9U2 



(4.14) 



Direction 3: 



^V.[h3(VV3-t3)].J-|L=o 



(4.15) 



where g-. , g« and g^ are the components of body force in directions 1 
2 and 3, respectively. 

4.3 Transport Equations for Turbulence Kinetic Energy and Length 
Scale 



The governing equations for the kinetic energy, K, and the 

turbulence length scale, L, can be stated as 

39 



[Convection] = [Diffusion] + [Sources] + [Dissipation] 

In vector notation, these equations are expressed as follows: (see 
also Equations 3.9 and 3.14) 



V • VK = 



V • VL 



Iv 

P 



Iv 

p 



(y^VK) + P|^ + e^ 



(y^vL) + Pl ^ ^L 



(4.16) 
(4.17) 



The expressions e^ and P, were defined previously in Equations 3.8 
and 3.12, respectively. However, the expressions Pj, and e in the 
curvilinear coordinates system will be discussed more thoroughly below. 
Gosman, et al . [2] derive the form of Pj. and e. in a curvilinear 
coordinates system from the following consideration: 



"Total Shear-Work" 
on Moving Fluid 



" Shear-Work Contribution 

to the Kinetic Energy of 

Mean Motion 



Shear-Work Contribution" 
to the Turbulence 
Kinetic Energy 



The total shear-work rate per unit volume can be written as 



f L. = V • [tiV, + %V« + t^Vj = u (W + W^) 
shear "-11 2 2 3 3-^ *^e^ m t^ 



(4.18) 



The shear-work rate is made up of W and of W4., which contribute to the 

m t 

kinetic energy of mean motion and turbulence kinetic energy, respec- 
tively. Details of the expressions W and W. are given in Reference [2]. 
Since the present turbulence model is directly concerned with the term 
Wx, it is necessary to write out W. in the following form [2]: 



40 



W, = 2 



'^(V^/h^)' 



LI 



3u. 



'3(V2/h2)^ 


2' 


■+ 


\ 3 
h2 9U2 







'2 9 



fv 1 

^2 



h, 8u, 



-12 



'^3 8 




-,2 

+ 


"^3 3 
^2 ^"2 


^3 



t2 



(4.19) 



Referring to Equations 3.9 and 3.14, the production rate of the turbu- 
lence kinetic energy per unit mass of the fluid is therefore 



K t t y t 



(4.20) 



and the dissipation rate of L is 



L 11 R ^ t 



(4.21) 



Substituting Equations 3.8, 3.12, 4.20 and 4.21 into the corresponding 
terms of Equations 4.16 and 4.17, respectively, the final form of 
Equations 4.16 and 4.17 appears as follows: 



/? 



V • VK = V • (v^VK) + C /K L W^ - Cj^^-^ 



(4.22) 



V . VL = V. (v,VL) - C^Cj^ ^ W^ + Cp/K 



(4.23) 



where W^ is given in Equation 4.19. 



41 



4.4 Vorticity and Stream Function Equations 

The procedure for obtaining the velocity components from the 
momentum equations. Equations 4.13 through 4.15, involves the solution 
of the Poisson equation for pressure [37] so that the continuity equa- 
tion, Equation 4.1, is satisfied. For two-dimensional flow situations, 
a stream function, ip, and a vorticity, w, are frequently used to reduce 
the number of variables from three (two velocity components and pres- 
sure) to two {^p and w) where no prediction of the pressure distribution 
is needed. The vorticity, w, and stream function, i|j, variables are 
used in this study and are formulated as follows. 

Defining the stream function and vorticity with the relationship 



h-w,^, (^-2^) 



'Z-^J^ (^-25) 



1 



^^2 



dh^V^ 9h^V^" 



3U.J du^ 



(4.26) 



the corresponding stream function and vorticity equations are [2]: 



y% = -pw (4.27) 



and 



42 



(P^") = h^ {3^7 [v • (-^2^2)] - 3^ [7 . (U,y] - ^ piL 



T„ ah^ 



h. au2 



3 


rT„ 3h^i 


3 


T22 31^2" 
hp 8U2 


9^2 


^22 ^^2 



(4.28) 



For two-dimensional flows, the vorticity vector, w, is always 
perpendicular to the plane of the flow and thus behaves as a scalar. 
Gosman, et al . [2] show that Equation 4.28 can be written as: 



V • (pVoi) = V • [V(y^aj)] + S 



(A) 



(4.29) 



where S is a source term to be discussed later. Using two vector 



to 



operators defined as follows [2]: 



v-xlr^^^ 



* X 3z z 9x 



(4.30) 



where the unconventional symbol V is used in this report for writing 
convenience; e and e are the unit vectors in the direction of the 
horizontal coordinate, x, and of the vertical coordinate, z, respectively, 



The term S is expressed as [2,11] 



(A) 



r-*- -^y 



S^ = 2[v(ej^ . V) • ?(e^ • V y^) + Vie^-V) • y(e^ • V y^)] (4.31) 

The expression of Equation of 4,31 in orthogonal curvilinear coor- 
dinates is (see Appendix): 



43 



0) 



^^2T= - 



8V^ 9y^ 9V^ 3y^ aVg 3y2 ^^2 ^^2 



8u^ 8U2 au^ au, 9u^ 9U2 Su^ 9u, 



+ iii 



39 3^2 



88 3^2 



3u-| au^ 



aup aun 



+ li. 



36 'h 



38 ^^1 



au^ au^ 



au^ au2 



+ V. 



ae ^^2 



ae ^^2 



aup au-j 



au^ au2 



+ V, 



ae ^^1 



ae ^^1 



au, au2 



au2 au. 



(4.32) 



where 



^1 ■ h^ au^ 



Po = 



ho auo 



and e, the angle between e^ and e^ , isdefined in the Appendix, Figure A.l 

The above equation implies that if the effective viscosity y 
is uniformly distributed, i.e., the derivatives of y with respect to 
both u, and u^ are zero or are negligible compared with the other 
terms in the vorticity transport equation. Equation 4.28, the S term 
can be wholly omitted from Equation 4.29. 

4.5 Analysis of Flow Parameters 

Since the problem considered here is atmospheric motion over 
bluff bodies, one can choose the characteristic body height, H, and the 
undisturbed upwind velocity at this height, Vm» as the characteristic 

rl 

length scale and velocity scale, respectively. 
Defining the dimensionless variables 



44 



^ = 



pHV, 



w == 



ojH 



K = 



L = 



V = 



H 



^e " pHV 



H 



(4.33) 



the resulting dimensionless form of the governing equations (Equations 
4.27, 4.29, 4.22 and 4.23) becomes: 



1 



h^h^ 



9 


^2 3$_ 


+ 3 


pi mJ 


8u^ 




aug 


h2 ^"2 



= -0) 



(4.34) 



3^(^2-Vi)-3U^(hi^V2)-3^ 



^2 



3U, 



^1 3 ,,. > 



- h,h„ S = 
1 2 w 



(4.35) 



3^(h2KV^).3^(h^KV2)-3^ 



'2 dK 



t z 9u 



9u. 



1 3K 



^ h2 ^"2 



- h^2 \ = 



(4.36) 



aU7(^2LVi)-9^(hiLV2)-3^ 



'2 8L 



^t ir 9u 



Ij 



9u- 



^t -- 
^ h. 



'1 9L 



du, 



- h^h^ \ = 



(4.37) 



45 



Considering that jig = (y^ + y)/pHVH5 one can rewrite Up as 



=CKl/2i:+T_ 
^e ^\i'^ "- Re 



(4.38) 



where 



Re = 



PHV, 



It is obvious that for turbulent flows where the Re number is yery 

'"1/2- 
large, y reduced to C K ' L and the effect of the molecular viscosity 

on the flow field can be ignored. The flow is thus insensitive to the 

Reynolds number based on molecular viscosity. 

In a neutrally stable atmospheric boundary layer, the wind 

profile over flat homogeneous terrain may be described by a logarithm 

law [4]: 



V = — In 

k: 



where 



1 +^ 



(4.39) 






z = 



'o = T 



The parameters V* and Zq are the friction velocity and surface roughness 
length scale, respectively, and von Karman's constant k is taken as 0.4. 
The friction velocity V* appearing in Equation 4.39 can be thought of 
as a measure of the turbulence level in the approaching flow (see 
Equation 3.26). 

46 



Examining Equations 4.35 through 4.39 and being cognizant of the 
above argument, the basic parameters involved in low-level atmospheric 

motion, where Coriolis effects are negligible, are found to be y, , V^, 

~l/2~ 
and Zq, where jl. = C K ' L. The first is a measure of turbulence trans- 
portation due to turbulence fluctuations, whereas the second and third 
characterize the turbulence level in the undisturbed upstream flow. 



47 



5.0 NUMERICAL ALGORITHM 

The governing equations. Equations 4.34 through 4.37, given 
in Section 4.0, are solved by utilizing the numerical procedure 
developed by Gosman, et al. [2], An outline of the procedure is 
given in Section 5.1. Theoretically, the procedure can be applied 
to any curvilinear coordinate system, however, a semi elliptical 
cylinder on a plane surface, discussed in Section 5.2, is chosen in 
this study. This geometry simulates a hill in natural terrain for 
which parametric variations of the aspect ratio can be carried out. 
The boundary conditions for the dependent variables oj, $, K and L 
are given in Section 5.3. In Section 5.4 the accuracy of the numerical 
results is discussed. 

5.1 Derivation of Finite Difference Equations 

Following the method of Gosman, et al . [2], Equations 4.34 
through 4.37 can be written in the general form: 

. r h 



9u 



* 



d±_ 



du, 



9 L -^ 



du 



, h^ ^^1 



9u, 



"2 2 



advection 



diffusion 



+ h-j h^ d = 



source 



(5.1) 



where ^ stands for ij), w, K or L. " The corresponding functional coeffi- 
cients a, b, c and d are listed in Table 5.1. 

48 



Table 5.1 Coefficient Functions of Equation 5.1 



1 
1 



-0) 

-S 



0) 



-S 



K 



-S. 



The general differential equation. Equation 5.1, is replaced by a finite- 
difference equation in the numerical computation procedure. The finite- 
difference equation is obtained by integrating Equation 5.1 over a 
typical cell appropriate to a central node p surrounded by nodes E, W, N 
and S, Figure 5.1. The resulting finite-difference equation is [2,38]: 



a< 



r \ 



8$ 



du. 



- ^, 



w 



9^ 
8Uo 



w 



Au« - 






du 



- 0. 



/■ N 



9$ 



3u- 



Au- 



advection 



-^(Ccf)) 



l3Si 



- b 



w 



^(C0) 



9S1 



w 



+ dp(A§i)^„(A§2)„3 = 



source 



where 



(Ai2)„3- 



diffusion 



[9S2 J 



- b. 



^(ccf)) 

l3S2 J 



(ASi) 



ew 



(5.2) 



49 




Figure 5.1 Grid arrangement with respect to the central node P 



50 



Sn and Sp are new independent variables which are connected to the 
u^-Up coordinates by the relations: 

ds-| = h.dUn ; dS2 = h^dUp 

The approximation for the advection terms in Equation 5.2 employs an 
"upwind differencing" technique. This is a one-sided rather than a 
centered-space differencing, where the scheme is backward differencing 
when the velocity is positive and forward differencing when it is 
negative. This formulation of the first-order terms gives greater 
numerical stability than can be obtained with central differences [32] 
Using the upwind differencing method, the first advection terms is 
approximated by [2]: 



0E -2 "^ *P 2 



Where Aip = i>^^ - i^^^ 



^se ^ KhE ^h-'h^ h^ 



V = l^^NE -^ ip -^ ^N ■*■ ^E^ 



Thus the value of <!>« in Equation, 5. 2 appears as follows 



({,g = (J)p if All; > 



51 



4)g = (J.£ if Ai|j < 

The term {d\p/^u^)^ has been approximated by 



'9$ ^ 






e 



The remaining advection terms are obtained similarly and the resulting 
finite form of the advection terms is: 

[Advection] - A A - I A.A. (5.4) 

P P i=NSEW ^ ^ 

where 

1=NSEW '^^■*^ ' '^N*N ^ Vs * '^E^E ' Vw 

'^E = fl^'i'SEN * I'i'SENl^ ' '''SEN = % * "''SE " '''NE " '''n 

^ = ||^'''nWS * l%sl] : %$ = '''N "" "''NW ■ '''SW - H 



'^N ^fl^^NW "" I'i'ENwl] ' hm'\''\^-^m-% 



S.r~. 



'^S = sf^E^ I%eI] • ""WSE = ""W + fsW - % - '''e 



Ap = A^ + Ay + A„ + As 



The above definition of the symbol ^-^Mcpu ts used throughout this study. 
Roache [37] indicates that the application of "upwind differencing" 



52 



technique may introduce an artificial viscosity. This pseudo viscosity, 
however, can be reduced by decreasing the grid size used in the numerical 
programming. The influence of the pseudo viscosity on the flow field is 
discussed in Section 5.3. 

A central -difference method is used to approximate the diffusion 
terms in Equation 5,2. It is assumed that the b's and c's can be 
approximated with a linear variation, i.e., 

bg « j[b^ + bp) 

^:r-(ccf))g ~ :: z vb.b) 

9s^ ^1,E " ^1,P 

The (AS]) and (As2)-,3 terms appearing in Equation 5.2 are approximated 
as 

(Ail)e„--^(Sl,E- St^„) (5.7) 

and 

respectively. Writing the remaining b's and c's in the forms of Equations 
5.5 through 5.8, the resulting finite-difference form becomes: 



[Diffusion] = I 

i=NSEW 



B,(c(J)),l-Bp(c4))p (5.9) 



1 ^-^'1 



where 

53 



^E " 4 ' ~ 

^1,E " 5i^p 



^W 4 ^^ ^~ 



^N "- ^ ^1,E ~ ^UW 
^N = 4 5 r 

^2,N " ^2,P 



^2,P " ^2,S 



Bp = B, + B^^ + B^ + B3 



Introducing Equations 5.7 and 5.8 into the source terms, the following 
finite-difference equation is obtained. 

[Source] = 1 dp(i^^^ - ii^„)(S2,N " ^2,5) (5.10) 

Summarizing, the finite-difference form of Equations 5.4, 5.9 and 5.10 
for integration of Equation 5.1 is 

<Pp = I C.<p. + D (5.11) 

^' i=NSEW ^ ^ 



where 



C. = (A- + B.c^OAAB 



D = -d^V /ZAB 
P P 



54 



1 -- ^ ^ ~ 



lAB = 



i=NSEW ^ P ^ 



Because of the nonlinear characteristics of the resulting finite- 
difference equation, Equation 5.2 is solved by an iterative 
successive substitution, technique. 

5.2 Elliptical Cylinder Coordinate System 

The orthogonal- curvilinear coordinates system used in the study 
is an elliptical cylinder system. This choice allows one of the 
coordinate lines to represent the solid boundary of a semielliptical 
body which in turn is assumed to simulate a hill. 

Mathematically, a two-dimensional semielliptical surface, A-B 
as shown in Figure 5.2(a), can be defined by the following equation: 



^ + ~ = 1.0 (5.12) 



where D is one-half the longitudinal length of the semielliptical body 
and H is the body height. In two-dimensional elliptical coordinates 
(u^. Up), Figure 5.2(b), the coordinates x and z can be expressed as 

X = -a cosh Ujp cos u, 

z = a sinh Up sin u-, (5.13) 

where a is one-half the length between the two focus points of the body; 

Up > and < u-. < 2-n. 

55 



Outer Boundary 




(1,411)^ wall (1,1) 



L D ^^^n? wall "(61', 41 ) 



(a) Schematic diagram of the numerical coordinates used 



u, =itt/2 



u^ = 2Tr/3 



u-| = 3Tr/4 




u. = 5tt/6 



{-a,o) (a,o) 



U-j = 27T 



(b) Two-dimensional semielliptic coordinate system 



Figure 5.2 The numerical and physical coordinate system for calculating 
the problem of flow over semiellipses 



56 



From Equations 5.12 and 5.13, it is evident that 

H = a sinh u^ 

D = a cosh Up (5.14) 

and therefore, a value [u2]u representing the solid ellipse boundary is 
determined by the following equation: 

[U2\ = coth"^ ^ (5.15) 

From Equations 5.14 and 5.15, the dimensionless form of Equation 5.13 is 
X = -cosh U2 cos u-j/sinh (u2)j^ 

z = sinh Up sin u-|/sinh (u2)|^ (5.16) 

The trace of the curvilinear coordinate surfaces of constant u-j and U2> 
respectively, on the x-z plane shown in Figure 5.2(b) are obtained from 
Equation 5.13. They are confocal ellipses and hyperbolas, respectively. 
The dimensionless metric coefficients h-i and h^ are as follows: 

h, = a/sinh2 u^ + sin2 u^ 



^2 = h 



where (see Equation 5.14) 
a = 1/sinh Lyi2\ 



57 



5.3 Numerical Coordinate System 

The origin of the numerical coordinates is situated at the up- 
stream corner of the semielliptical body, node A in Figure 5.2(a) with I 
indexing in the positive u, direction and J indexing in the positive u« 
direction. In the iteration process, the field was swept from left to 
right beginning at the outer boundary and proceeding in the decreasing 
J-direction. A total of 2,501 grid points (I x J = 61 x 41) are con- 
structed in the flow field. The variable mesh used is determined by 
the following equation: 



[u^lj = (I - l)Tr/60 ; I = 1, 2, 3, 



61 



[U2JJ = [U2]j=-, + I C^(C2)"'^ ; J = 2, 3, 4. ••• 41 (5.17) 
n=2 

where 

t"2]j=l = coth""" ^ 



and C-j and C^ are constants. A value assigned for C-, , then, represents 
the distance of the first grid point to the obstacle through Equation 
5.16. A value greater than unity assigned for Cp, however, is used to 
stretch the grid size and therefore changes the distance of the outer 
boundary from the obstacle. For D/H = 2.0, Equations 5.16 and 5.17 show 
that the outer boundary is a distance of approximately 30 H from the 
obstacle if C^ = 0.025 and C^^ = 1.05. The distance becomes 45 H if C^ 
is changed from 1.05 to 1.055. The distance from the first grid point 
to the obstacle is approximately 0.05 H in these cases. A distribution 



58 



of grid points in the x-z plane is shown in Figure 5.3 for the case 
D/H = 2.0, C^ = 0.025 and Z^^ = 1.05. 

5.4 Boundary Conditions 

The dimensionless form of the boundary conditions is discussed 
in this section. The outer boundary conditions are determined from the 
undisturbed flow conditions, whereas the obstacle and the wall boundaries 
are determined from the local equilibrium conditions. 

5.4.1 Outer boundary . It is assumed in the present study 
that the outer boundary. Figure 5.3, is located far away from the semi- 
elliptical body on the half plane so that the flow conditions at the 
outer boundary are not disturbed by the presence of the obstacle in 
the flow field. For atmospheric flow as considered in the study, the 
undisturbed wind velocity profile is the logarithmic velocity profile 
given by Equation 4.39. It is (Equation 4.39): 



V = — £n h + ^ 



where von Karman's constant k is taken as 0.4. In the Cartesian co- 
ordinates system, the undisturbed stream function $(z) is obtained by 
integrating this velocity profile over the height of the solution region, 



'^& = 



V^z. 



V dz = 



(1 + z^) £n(l + i^) - l^ 



(5.18) 



59 



en 

o 




Figure 5.3 Distribution of grid sizes for D/H = 2.0,. 
C^ = 1.05 



C^ = 0.025, 



where z 1s the ratio z/z . For a grid node, Q, on the outer boundary. 
Figure 5.3, the height, Zq, can be determined from its coordinates 
[(u-jJq, ("2^0^ ^y ^^^ following equations (see Equations 5.13 through 
5.15): 



sinh(u2)Qsin(uT)Q 
^ sinh[coth"^ (D/H)] 



(5.19) 



where (u])q and (u2)q are the values of u-, and u^ at the position of 
the node Q, respectively. At the node Q the boundary condition of the 
stream function ijj^ is therefore (see Equations 5.18 and 5.19): 



^Q " ^^^Q^ 



V^z 



L^l 


in 


L!al 


^Q 


[ ^oJ 




[ ^oJ 


^0 



(5.20) 



Similarly, all other boundary conditions of Tp for the nodes on 
the outer boundary are obtained. This procedure is also used to obtain 
the outer boundary- conditions for w, k and L, respectively. Thus, only 
the undisturbed conditions for w, k, and L, respectively, of the 
atmospheric boundary layer written in terms of the z are presented in 
the following. Equation 5.19 is used to convert (u-,,Ur,) to the (x,z) 
coordinate system. 

In the undisturbed atmospheric boundary layer, the dimensionless 
vorticity ai is written as follows: 



w 



(2)^-#=- 



3V 
9z 



-^oP 'I 



oJ 



(5.21) 



61 



The boundary condition for the turbulent length scale is based- 
on a mixing-length hypothesis. In the neighborhood of the wall, the 
mixing length is generally assumed to be linearly related to the normal 
distance from the surface. For a wery rough surface, the mixing length 
hypothesis does not go to zero but approaches a size on the order of the 
roughness length, z . The following relationship is assumed: 



L = k(z + z^) (5.22) 



The boundary condition for the turbulence kinetic energy, K, is derived 

from the assumption of constant shear layer. In terms of the Prandtl- 

Kolmogorov formula. Equation 3.5, the shear stress, x, of the undisturbed 
atmosphere is written as follows: 



pV*^ = C^p/K L 



faVl 



9z 



Substituting Equations 4.39 and 5.22 into the above equation, the 
following equation is obtained: 



K = 



V* 



C 



(5.23) 



5.4.2 Wall boundary . The wall boundary is assumed to be a 
solid wall; therefore, no flow passes through it. Thus, the wall boun- 
dary forms a streamline which is assigned the value zero. 

$ = (5.24) 



62 



Applying the nonslip conditions at the wall; i .e. , V, = 0; Vo = 0; the 
03 at the wall becomes (Equations 4.26 and 4.33): 



OJ 



wall 



9Vj 



1 ^h^ 






wall 



In the present study, the wall boundary is constructed by u, =0, 

u = TT and u« = {up)^ (see Figure 5.2 and Equation 5.15). Considering 

that 3V,/3U2 = along u, = and u, - tt, and that BV^/Bu-j = along 

Up - (^9)h> ^^^ dimensionless vorticity along the wall can be written as 

follows: 



03 



wall 



l!!2 



wall 



if u. = or U-. = TT 



and 



0) 



1 



9V 



wall 



h2'"2 



wall 



if Uj = (U2)|j 



Since the above equations can be written in a general form when the 
velocities V, and Vp are replaced by Equations 4.24 and 4.25, respectively, 
only the expression of wall vorticity along Up = (Up). , w^^, in terms of 
stream function (see Equation 4.24) is discussed in the following: 



- = 1 dU) 



U2=(u2)b 



(5.25) 



63 



To express Ja, in terms of known variables, the stream function around 
a nearby wall point, point NW in Figure 5.2a, is expanded in a Taylor 



series: 



^. 



= i,.^^ 



W 9u, 



^ 1 9^iJ 
Au« + ^ — ^ 



W 



? 9 



(Au^) + H.O.T. 



(5.26) 



W 



The second term on the right-hand side of Equation 5.26 is zero due to 
the nonslip condition. From Equations 5.25 and 5.26, the vorticity at 
node W on the wall boundary; w^; can be approximated by the following 
equation if the higher order terms of Equation 5.26 are neglected. 



2(^MU - $K|) 



0). 



^ m( 

'W h2(Au2)' 



W' 



(5.27) 



Assuming that the concept of a constant shear layer is applicable in the 
region very close to the wall, the kinetic energy of turbulence, K, and its 
length scale, L, take the following form, respectively [11]. 



I ICO) z 



K = 



'0 



^ J 



L = K z. 



(5.28) 
(5.29) 



5.5 Accuracy of Numerical Solutions 

The accuracy of numerical solutions is generally affected by the 
imposed outer boundary condition and the grid sizes used in the solution, 



64 



Theoretically, the limit of the outer boundary is at infinity. In a 

practical analysis, however, this cannot be imposed; therefore, in this 

study the conditions at infinity are applied at a large distance, z^gj.* 

from the elliptical body. The value of z„^^ is determined in such a way 

max 

that further increase in its value does not significantly effect the 

numerical solutions. The value of z„ is determined by numerical experi- 

max ^ 

ment. In the case D/H = 2.0, solutions obtained from z^^^ = 30 H and 

max 

45 H, respectively, are compared in Figures 5.4 through 5.7. This com- 
parison indicates that the difference between these two solutions is 
small (about 5% of the local value). The solutions for the case D/H = 2.0 

shown in Section 6.0 are therefore computed with z ,, = 45 H to assure 

^ max 

good compliance with the boundary conditions applied at the outer boun- 
dary. Equations 5.18 and 5.21 through 5.23. However, a value of 
z„3^, = 30 H would provide the necessary accuracy if desired. In the 

maX 

present study, the numerical computation for D/H = 1.02 is carried out 

for z^^^ = 30 H. 
max 

In addition to the outer boundary condition, the grid size used 
in a numerical procedure can significantly affect the accuracy of the 
solutions. Two cases are used to investigate the effect of the distri- 
bution of grid size on the numerical solutions. The first case is 
C-. = 0.05 and Cp = 1.02 in Equation 5.17, and the second case is 
C, = 0.05 and Co = 1.05. As mentioned in Section 5.3, the number 
of grid points constructed for numerical computation is 2,501 points 

(I X J = 61 X 41), which yields 2 ^^ = 30 for both cases mentioned 

max 

above. The number of grid points in the wall region of the second case 
is approximately twice the number of the first case. In the region 

65 







-3 



-2 



-1 






1 


x/H 




' ^max 


-- 45 H 


* ^max 


= 30 H 



Figure 5.4 Streamline, i|j, patterns of flow over a semiellipse for D/H = 2.0, 
(below i = 0.6 solutions coincide) 







z = 45 H 
max 



: z = 30 H 

max 



Figure 5.5 Distribution of turbulence kinetic energy, K, for D/H = 2.0 



en 

00 




: r = 45 H , 

max. ' 

=^max =30H 



Figure 5.6 Distributions of turbulence length scale, L, for D/H = 2.0 



en 




: z = 45 H, 
max 



: z = 30 H. 
max 



Figure 5.7 Vorticity, w, patterns for flow over a semiellipse for D/H =2.0 



near the outer boundary, more grid points are constructed for the first 
case than are constructed for the second case. Again, the effect of 
the grid size on the numerical solutions is small for the cases compared, 
Figure 5.8. Based on the above discussion, it was assumed that the grid 
distribution given by Equation 5.17 with C-j = 0.025 and C2 = 1.055 
would provide meaningful results. Therefore, this distribution is used 
throughout the study. The application of C-. = 0.025 and C2 = 1.055 
reduces the grid size to 0.05 H at the wall region and approximately 
7 H at the region near the outer boundary if D/H = 2.0. For D/H = 1.02, 
the grid size at the wall region and at the outer boundary region is 
0.025 H and 5 H, respectively. The distance from the outer boundary to 
the obstacle is approximately 45 H if D/H = 2.0 and is approximately 
30 H if D/H = 1.02. 

As mentioned in Section 5.1, the "upwind differencing" technique 
used in this study causes a pseudo viscosity, the value of which depends 
on the grid size [37]. Since the grid size near the obstacle is very 
small (not greater than 0.05 H), the effect of the pseudo viscosity is 
considered to be small in that region. An analytical method of 
determining the error caused by the pseudo viscosity effect is not 
available for nonlinear equations; the reader must therefore make his 
own judgement of this effect from the grid distribution shown in 
Figure 5.3. From previous experience, the author believes that this 
effect will not cause significant inaccuracies in the results. 

The convergence criteria for the numerical computation is 









< 10'^ (5.30) 

max 70 



20 



15- 



"max 



45 


0.025 


1.055 


30 


0.025 


1.050 


30 


0.050 


1.020 



z/H 10- 



5- 



0:5 




2.0 



Figure 5.8 Velocity profiles on top of a semiellipse for 
D/H = 2.0 and z^ = 0.005 



71 



M 



where 4) is the dependent variable at grid point P, and the superscript n 
denotes the nth iteration. Equation 5.30 insures good convergence in 
areas of small (f)p values. The present study was run on an IBM 370/3031 
computer. Execution of the present computer code requires approximately 
256 K bytes of main storage. Each iteration of the code requires approxi 
mately 9.5 sec on the central processing unit of the IBM 370/3031. The 
computational time required for convergence of the cases run in this 
study is approximately 180 minutes for each case if the undisturbed flow 
conditions, i.e., outer boundary conditions described in Section 5.4 
are the initial conditions. The computation time is reduced to 120 
minutes if the converged solutions of a case are used as the initial 
conditions of other cases. 



72 



6.0 RESULTS AND DISCUSSION 

The numerical results obtained in this study represent solutions 
of atmospheric flow over a two-dimensional semi elliptical body. The 
effects of aspect ratio, of distributions of surface roughness, and of 
molecular Reynolds number are discussed in this section. The numerical 
results are compared with the available experimental data. Since most 
currently available experimental data are for flow over rectangular- 
shaped bodies, i.e., fences, blocks and forward-facing steps, numerical 
results computed with the present model in a Cartesian coordinates system 
over these types of geometries are first tested. The agreement between 
experiment and theory is excellent. Since the computation algorithm used 
to compute flow over the semiellipse is the same as that for Cartesian 
coordinates, the present computed results for flow over semiellipses are 
believed to be good even though there is insufficient data available 
for a valid comparison. The flow separation created by a rectangular 
geometry, however, is quite different from separation caused by a 
semiellipse. The differences in the mechanism of flow separation are 
discussed in Section 6.5. 

6.1 Characteristics of Flow About a Two-Dimensional Semielliptical 
Obstacle 

The results discussed in this section are obtained for a semi- 
elliptical body with aspect ratio D/H = 1.02. The outer boundary is 
situated approximately 30 H away from the body. The dimensionless sur- 
face roughness scale, z« = 0,005, is uniformly distributed along the 

73 



wall and over the sem1 ell ipse. Substituting Zq = 0.005 into Equation 
4.39, the dimensionless friction velocity is obtained as V^ = 0.075. 
From Equations 3.27 and 5.23, the undisturbed upstream turbulence 
energy is therefore K = 0.0325. The molecular Reynolds number. Re, is 
assumed to be large enough so that its effect on the solutions is 

negligible, as discussed in Section 4.5. This assumption is reason- 

5 
able when Re is greater than 10 , as found from the present numerical 

study. The effect of Re is discussed later. 

The streamline pattern for Zq = 0.005 and D/H = 1.02 is shown 
in Figure 6.1. The flow is seen to converge as it passes over the 
obstacles. The convergence of the flow causes the wind to gain speed 
in the region above the obstacle. Figure 6.2. The constant velocity 
contours in Figure 6.2 show that the local maximum in velocity occurs 
at 0.1 H above the crest of the body and the local minimum occurs at 
1.0 H. Similar phenomena have been observed in full-scale field 
measurements [9]. More discussion of this phenomenon is given when 
the numerical results are compared with the experimental data in 
Section 6.4. 

Figure 6.2 shows very strong gradients occurring in the area 
close to the surface of the obstacle. This implies that a high produc- 
tion rate of turbulence kinetic energy, K, occurs there, as discussed 
in Section 3.0. The production rate of K in that region is stronger 
than its dissipation rate and, therefore, relatively high turbulence 
kinetic energy occurs in this region, as shown in Figure 6.3. Down- 
stream of the obstacle, the region of high turbulence expands grad- 
ually. Figure 6.3, due to the effects of convection and diffusion. A 

74 






2/H 




-.»L,,,,,L 



Figure 6.1 Streamline, ^, patterns over a semiellipse for D/H = 1.02 



and z = 0.005 





z/H 






1.14 



1,1 



1.05. 
1.0 ^, 

0.9 

0.8 

0,6 -- 



^>. 1.17 ^^" 



.1.14 
^1.10 



/ 



^1.05 






0.9 
-0.8 



0.4 



l^^-^'O, 



-3 



-2 




^-0.6 



(^-C::^^:::- 



0.4 

1_ 



Figure 6.2 V/V^ isotachs over a semiellipse for D/H = 1.02 
and Zq = 0.005 






z/H 



0.055' 



0.05 



i 
\ 
\ 

0.06 
s 



" .»--• 



\ 
I 

0.055 
/ 



/ 



0.06 



0.045-^^^ "'"^-^,''^ 

0.04 — rr""""---^ 



0.01. .---• 



^^-''^ 0.07.^- 
-^ ""0.06 



-3 



-2 




Figure 6.3 Distribution of turbulence kinetic energy, K, over a semiellipse 
for D/H = 1.02 and z^ = 0.005 



relatively low turbulence area is observed in the region approximately 
1 H to 2 H above the crest, as is expected due to the nearly uniform 
velocity profile. Figure 6.2. In the vicinity of both the upstream 
and the downstream corners of the obstacle, the turbulence kinetic 
energy is very small due to the motion of the fluid being restricted 
by the wall . 

The distribution of the turbulence length scale is shown in 
Figure 6.4. This figure indicates that the length scale L decreases as 
the flow approaches the obstacle. Then L increases rapidly as the flow 
passes over the obstacle. Downstream of the obstacle, the L decreases 
again with the exception that L increases in a region immediately 
behind the obstacle. The rapid increase of L on top of the obstacle 
implies that the length scale cannot be prescribed in that region by 
a mixing length hypothesis. At 0.5 H above the crest, the length scale 
L is equal to 0.4, which is the same value as L upwind at 1 H above the 
plane surface. The present results indicate that L would be under- 
predicted by a mixing length hypothesis in the region above the 
obstacle. 

Constant vorticity contours are shown in Figure 6.5. )lery 
strong vorticity is generated on the surface of the obstacle due to the 
nonslip condition at the wall. The distribution of vorticity in the 
flow field is determined by the velocity gradients. Convection and 
diffusion of vorticity are clearly shown in Figure 6.5. 



78 



to 



z/H 



■-^ 0.? 



^ 



"'"^-'-^i)^5 







\ 






Figure 6.4 Distribution of turbulence length scale, L, over a semiellipse 

for D/H = 1.02 and z = 0.005 
' 



z/H 




Figure 6.5 Vorticity, w, patterns over a semiellipse for D/H = 1.02 
and z^ = 0.005 



80 



6.2 Effect of Aspect Ratio 

The effect of variation in aspect ratio on the flow field for 
a nondimensional surface roughness Zq = 0.005 is shown in Figure 6.6. 
The V/Vu isotachs for the aspect ratio D/H = 1.02 and 2.0 are presented. 
One observes that the velocity above the crest of the obstacle increases 
with a decrease in aspect ratio. This effect is very strong in the 
lowest 1 H layer above the crest and is felt even up to 10 H, Figure 
6.7. Reference [19] explains that a stronger favorable pressure gradi- 
ent occurs near the leading edge of a smaller aspect ratio semi- 
elliptical body, thus causing higher wind speeds on top of the body. 

The constant contour lines of the dependent variables $i, co, 
K and L for D/H =2.0 were shown in Figures 5.4 through 5.7. Comparing 
these figures with Figures 6.1 and 6.3 through 6.5, one observes that 
the vorticity, oj, produced near the surface is significantly affected 
by the aspect ratio. However, the effect on the turbulence kinetic 
energy is small. The reason for this is that by definition the vorticity 
is a direct function of the velocity gradient; however, the local tur- 
bulence kinetic energy, K, is a balance between the diffusion, the pro- 
duction, the dissipation and the convection of K. Therefore, the local 
turbulence kinetic energy is dependent not only on the velocity gradient 
but also on the local turbulence level and the surface roughness. The 
effect of aspect ratio on K is small for the cases compared in the present 
study because the surface roughness has been held constant (?« = 0.005) 
along the wall. Moreover, the undisturbed turbulence kinetic energy also 
remains the same, i.e., K = 0.0325, for the same surface roughness value. 

81 



00 




: D/H 
: D/H 



2,0 
1.02 



Figure 6.6 V/V^ isotachs over semiell ipses for Zq = 0.005 



z/H 10- 




2.0 



Figure 6.7 Velocity profiles on the top of semiellipses of different 
aspect ratios for z = 0.005 





83 



6.3 Effect of Surface Roughness 

The effect of surface roughness on the flow field 1s described 
in two parts: (1) the effect of the local value of z, assigned to the 
obstacle and (2) the effect of the value of z^^ assigned to the upstream 



plane. 



The effect of the local z^. is studied for the aspect ratio 



D/H = 2.0. The value of Zq used to determine the undisturbed friction 
velocity, V^, Equation 5.14, is 0.005 in this case. The value of V^ is 
0.075, and therefore the undisturbed upstream K is 0.0325, Equation 5.23. 
Numerical results are obtained for z. = 0.01 and for z. = 0.005, where 
z. is the surface roughness scale of the obstacle. For increasing z. , 
the velocity decreases in the lower layer directly on top of the 
obstacle. Figure 6.8. The thickness of this layer is approximately 
1 H. The increase of surface roughness causes a larger drag on the 
flow field. This friction drag reduces the wind velocity in the lower 
layer. Figure 6.9, and thus displaces the streamline to a greater 
height, Figure 6.10. Since the kinetic energy lost in the mean wind 
field due to the effect of friction drag is converted into the turbulent 
fluctuations, the turbulence kinetic energy is increased with increasing 
z, , Figure 6.11 . 

The flow field about an obstacle is influenced not only by the 
roughness scale on the obstacle but also by its distribution along the 
wall. If the roughness scale is increased from z^ = 0.005 to z-. = 0.01, 
i.e., V^ increases from 0.075 to 0.087, the turbulence kinetic energy, 
K, in the upstream undisturbed flow increases from 0.035 to 0.044. 



84 



4- 












z 


Zu / 









/ 




0.005 


0.005 / / 


3- 




0.005 


0.010 / / 






0.010 


0.010 / / 
/ / 
/ / 
/ / 


z/H 2- 




J 


/ / 
/ / 
/ / 
/ / 
/ / 


1- 




/ 


/ 






/ 


/ 






4^ 


/ 






^ A 








-^ -J? 




- 











9 


r.o 1 


.'lo ^\^ 1. 



v/v, 



Figure 6.8 Velocity profiles on top of a semiellipse with different 
surface roughnesses for D/H = 2.0 



85 



00 



z/H 




z^ = 0.005 
z, = 0.01 



Figure 6.9 V/V|^ isotachs over a semiellipse for D/H = 



2.0 and z^ = 0.005 



00 
00 




-3 



-2 



-1 




x/H 




1 


. 


^b = 


■- 00 






: 


^b = 


0.01. 



Figure 6.11 Distribution of turbulence kinetic energy, K, over a semi- 



ellipse for D/H = 2.0 and z^ = 0.005 



The K in the flow field increases correspondingly, as one observes 
in Figures 6.11 and 6.12. From Equation 5.18, the maximum value of 
Ijj assigned at the outer boundary of the computed flow region, i.e., 
the value of ij) at z = 45, is 68.4 if Zq = 0.005 and V^ = 0.075, and 
is 76.2 if Zq = 0.01 and V^ = 0.087. Therefore, more mass flow is 
brought into the flow region computed for the case of larger z« and 
thus larger velocity in the flow field is obtained, Figure 6.13. 

The velocity profiles on top of the obstacles, shown in 
Figure 6.8, indicate a slower wind speed in the layer 0.2 H above the 
crest for 1^ = 0.01 than for Zq = 0.005. This is obviously a result 
of the larger drag force created by the rougher surface. In the layer 
between 0.2 H and 1.0 H, V increases for larger z„ due ta t'te enhanced 
turbulence mixing process which carries the upper layer fluid elements 
with high momentum into the lower layer. This interaction augments 
the momentum in the lower layer. 

Figure 6.14 illustrates that the velocity speed-up ratio, 
defined as the local speed at a height z above the crest of the obstacle 
divided by the undisturbed speed at the same height above the upstream 
surface, i.e., S = V(z)/Vq(z), increases with increasing z«. This 
agrees with the results of previous studies [32]. As explained pre- 
viously in Section 6.2, the aspect ratio has a significant influence 
on the speed-up ratio. From the limited cases discussed, the effect 
of the aspect ratio on the speed-up ratio can be seen to exceed the 
effect caused by the change of surface roughness. Figure 6.14. 



89 



'■JD 
O 



z/H 




Figure 6.12 Distribution of turbulence kinetic energy, K, over a semi- 
ellipse for D/H = 2.0, z^ = 0.01, and z^ = 0.01 



lO 



z/H 2- 




z = 0.005 and z. = 0.01 
b 

i = 0.01 and z. = 0.01 
b 



Figure 6.13 V/V^ isotachs over a semiellipse for D/H = 2.0 



z/H 




1.0 



1.25 



1.50 



1.75 



2.0 



V(z)/V^(z) 

Figure 6.14 Wind speed-up ratio on top of semiellipses with different 
surface roughnesses and different aspect ratios 



92 



6.4 Verification of Results 

Comparisons of the numerical results with experimental field 
data and with atmospheric boundary layer wind tunnel data are given 
in this section. Since the flow field is significantly affected by 
the geometry of the obstacle, the surface roughness scale and the 
upstream turbulence level (discussed in Sections 6.1 through 6.3), 
only experimental data obtained under well-defined conditions are 
used for comparison. 

Most data available for comparison are for flow over rectangular 
bluff bodies, i.e., steps, blocks and fences. Limited data are avail- 
able for flow over other geometries. In this section, therefore, 
comparison is first made for the case of flow over rectangular bluff 
bodies. For flow over semielliptical bodies, the numerical solutions 
are compared with the analytical results of Frost, et al . [19]. The 
present solutions are also compared with the results of flow over 
obstacles with geometries similar to semielliptical bodies. These 
cases are the analytical results of flow over a two-dimensional Gaus- 
sian hill [39], the wind tunnel data of flow over two-dimensional 
half- and full-sinusoidal bodies [7], the flow over a three-dimensional 
Gaussian hill [7], and the data of Bradley [9] obtained from field 
studies of flow over a full-scale hill. 

Numerical results for flow over a 90° escarpment computed with 
the current computational model in Cartesian coordinates [32] and the 
experimental wind tunnel data [6] are shown in Figure 6.15. The 
simulated undisturbed upstream flow in the wind tunnel corresponds to 

93 



x/H=0 x/H = 3.0 x/H=6.0 x/H=10.0 
x/H = -3.0 



z/H 












^imaiiim 1 ;o 'W/z/xwi ;q' w//^/(t?' i .0 'z'^^^^''/^*'^ 



V(z)/V^^(z) 
0.5 1.0 1.5 



: the present computed result 



: numerical results of Reference [32] 

• : wind tunnel data [6] 

Figure 6.15 Wind speed-up ratio over a 90° escarpment 



94 



atmospheric boundary flow over ground with ?« = 0.1 m If a model scale 
1 :'300 is assumed [6]. Since the height of the model in the wind 
tunnel study is 50 mm, at 1 : 300 this corresponds to a full-size scale 
body of 15 m. Accordingly, the dimensionless surface roughness 2^ is 
0.007. The numerical results shown in Figure 6.15, therefore, are 
computed with z« = 0.007 uniformly distributed along the wall. 
The numerical solutions obtained by using the coefficients C , Cp, Cj, 
and Cr- given in Reference [32] are also shown in Figure 6.15. The 
results are relatively insensitive to the coefficients used. 

Some differences between the numerical results and the 
experimental data [6] are shown in Figure 6.15. On top of the front 
upper corner of the step, the velocity speed-up ratio, V(z)/Vq(z), is 
slightly over-predicted. In other regions of the flow field computed, 
the numerical results are under-predicted in the region near the 
ground. These observations suggest that the surface roughness scale, 
z^, used in current computation is higher than the experimental value, 
as stated in Section 6.3. 

Figure 6.16 shows a comparison between the computed value of 
/K/V^ and the experimental data for a^ /V^ [6], where V^ is the upwind 
velocity at a height 10 H above the ground and Oy is the r.m.s. value 
of the turbulence fluctuation of the horizontal velocity components. 
Over level homogeneous terrains, many micrometeorological experiments 
suggest that 



Oy - 0.5 ay 



a„ = 0.8 aw (6.1) 

^3 ^1 

95 



x/H=0 x/H=3.0 x/H = 6.0 x/H = 10.0 




-. ■ - ^ - - . K/V (or aJM ) 



-^/^oo ^^ ^^^ present computed results 

• ^^v-j/V^ of experimental results [6] 

Figure 6.16 Turbulence intensity over a 90° escarpment 



96 



where o^ and a^ are the r.m.s. values of the turbulence fluctuations 
of the vertical and the lateral components, respectively [9]. By the 
definition of K, 



1 f 2 ^ 2 ^ 2 

2 hi + ^V2 + ^V3j 



(6.2) 



and the relationships of Equation 6.1, the following approximate result 
is obtained: 

av //K = 1.03 (6.3) 

for atmospheric flow over uniform terrain at a neutral stability con- 
dition. However, in actual fact the ratio a^^/K is affected by the 
terrain over which the flow passes [40]. As flow passes over obstacles, 
the ratio cr^ /a^ increases slightly, while the ratio Oy foy remains 
approximately equal to upstream value [7,9]. Limited data [9] shows 
that on top of a hill, the ratio a^ /a^ Increases from its upstream 
value 0.5 to a local value of 0.8. This reduces the value of the ratio 
cTy //!< from 1.03, Equation 6.3, to 0.93. Since the value of the ratio 
Oy //K depends on the terrain over which the flow passes, the comparison 
shown in Figure 6.16 should be made qualitatively rather than quanti- 
tatively. Figure 6.16 shows that the present results of /K/V^ agree 
well qualitatively with the experimental data of o^ /y^ [6]. The 
analytical results of v^K/V^ obtained by using the coefficients given 
in Reference [32] result in an over-prediction of v'K/V^ of approximately 
five times the present results, in the region next to the wall. It 
is concluded that the coefficients developed in the present study result 

97 



in an equivalent prediction of the mean wind field and a better pre- 
diction of the kinetic energy field, as compared with those developed 
in Reference [32]. 

Results from the present numerical model applied to a block 
geometry, when compared with experimental wind tunnel data and full- 
scale field data, show excellent agreement in both the mean velocity 
field and the turbulence kinetic energy field [11,26]. It is noteworthy 
that the upstream turbulence level causes significant effect on the 
flow field computed. An increase in upwind turbulence reduces the dis- 
tance, X , between the block and the flow reattachment point. Figure 
6.17. Reference [41] reports that reattachment of flow over a solid 
fence occurs at a distance 10 H to 15 H measured downstream from the 
fence. This is accurately predicted by the present numerical model in 
the range of V^/V,^ = 0.053 to 0.015, Figure 6.17. 

Figure 6.18 shows typical V/V^ isotachs computed and measured 
for a solid fence. Excellent agreement between the computational 
results and the experimental data [41] is observed. In the recirculat- 
ing flow region, the present results show a region of negative velocity, 
as one would expect. This region of negative velocity, however, is not 
shown in the experimental data. The measurement of velocity using hot 
wire [41] in the region close to the fence, where the velocity is low 
speed and highly fluctuating, can cause error in the measurements. 

Bradley [9] measured wind speed profiles on the top of a hill. 
The terrain surrounding the hill is roughly a uniform surface with 
Zq = 0.005. The hill, however, is covered by a forest. The shape of 
Bradley's hill, along with the shapes of other surface geometries for 

98 







0.05 



0.10 



0.15 



v.,/v 



H 



:the present computed results 

O :Raine and Stevenson [41] ; D/H ^ 0.0 
A :Woo, et al. [5] ; D/H = 0.75. 
• .-Frost and Shahabi [8]; D/H = 0.75 



Figure 6.17 Effect of the turbulence of the approaching wind on the size 
of the wake zone behind rectangular blocks 



o 
o 



Z/H 




16 



10 12 14 
x/H 

: The present computed results, 

: Raine and Stevenson [41]. 



20 22 24 26 



Figure 6.18 V/Vq isotachs .downstream of a solid fence, z = 0.01 



which either empirical or analytical data are available for comparison, 
is shown in Figure 6.19. Measured and computed vertical profiles of 
horizontal velocity on the top of the respective hill geometries shown 
in Figure 6.19 are presented in Figure 6.20. The velocity scale V^ 
used in Figure 6.20 is the upwind velocity at z^ = 0.2 H, where Zp is 
a reference height. 

Great differences between the various analytical models is 
clearly seen from Figure 6.20. These differences are due to the 
different assumptions used in the analyses. The length scale used 
in References [19] and [39] has been prescribed by the mixing-length 
hypothesis. This approach can cause significant differences in the 
analytical results (see Section 6.1). The boundary layer approxima- 
tion of Frost, et al . [19] and the different closure form used by 
Taylor [39] for the shear stress terms of the governing momentum 
equations can also create disagreement between the analytical models. 

The data of Bradley [9] show wery low values of V/Vp in the 
region z < z^. This is believed to be caused by the existence of a 
forest on the hill. The existence of the forest increases the effec- 
tive surface roughness and therefore decreases the velocity in the 
region next to the ground (see Section 6.3). An effort was carried 
out to simulate the effect of a forest on the flow field by increas- 
ing the surface roughness scale, z. , on the semiellipse. The value 
0.1 was assigned to Zl, while the ground surface roughness, z^., 
remained at 0.005. No converged numerical solutions are obtained 
due to the large magnitude of the abrupt change in surface roughness. 
The numerical solutions for z, = 0.01 and Zq = 0.005 show only slight 

101 



o 
A 

D 



- 1:2 Ellipse 
Bradley hill [9] 
Full sinusoidal shape [7] 
Gaussian high hill [39] 
Half sinusoidal shape [7] 



o 




-4H 



Figure 6.19 Various shapes of surface obstacles for which either 

analytical or experimental data are compared in Figure 6.20 



30.0 
20. Ol- 
io. 



z/z, 



1.0 - 



0.2 




▼ 3-D Gaussian hill [7] 
A2-D half sinusoidal hill [7] 
2-D full sinusoidal hill [7] 
• Bradley [9] 

— The Present Computed Results 

— Frost, et al . [19] 

— Taylor [39] 



0.6 



0.8 



1.0 



1.2 



1.4 



1.6 



V/V 



R 



Figure 6.20 Vertical velocity profiles on the top of the respective 
obstacles shown in Figure 6.19 



103 



Ill mil ■laiiiii ■■ ■iiiiiiBinminiiiniin 



reduction of the ratio V/V^, in the region z < Zj,. An accurate method 
of simulation of the effect of a forest has not yet been developed. 

Figure 6.20 shows that the analytical results of Frost, et al . 
[19] agree well with the data of Bradley [9]. It is noteworthy, how- 
ever, that the analytical results are obtained by a two-dimensional 
analysis with uniform z^ [19], while the experimental data are obtained 
for flow over a hill with a forest on its top [9]. The wind tunnel data 
for a three-dimensional Gaussian hill [7] are also consistent with the 
computed results of Frost, et al . [19] in the region z > Zn, Figure 
6.20, As shown in Figure 6.20, the present computed results agree 
excellently with the wind tunnel data [7] of flow over two-dimensional 
obstacles with geometries similar to a semiellipse. In the region 
z < Zp, Taylor's predictions [39] agree with the present computed 

- K 

results and then shift to agree with the computed results of Frost, 
et al . [19] if z > Zn, Figure 6.20. However, the value of Zq used by 
Taylor [39] is 0.0005 rather than 0.005 used by Frost, et al . [19] and 
by the present computation. 

Although the limited data compared in Figure 6.20 are obtained 
from flow over obstacles with similar geometries, one cannot form any 
conclusions from the agreement (or disagreement) between the analytical 
results and the experimental results until more experimental data are 
available for comparison. This is due to the fact that the upstream data 
were not available for Bradley [9] and the flow field can be significantly 
affected by the upstream turbulence level, the magnitude of surface 
roughness and its distribution, and the geometry of the obstacle. This 
is demonstrated by the present study. 

104 



Figure 6.21 indicates that the distribution of turbulence 
kinetic energy is significantly affected by the aspect ratio of the 
obstacle, as well as by the magnitude of the surface roughness scale. 
As mentioned earlier, the length scale, L, on top of the obstacle is 
under-predicted by a prescribed mixing length hypothesis. An under- 
prediction of L causes an under-prediction of the production rate of 
turbulence kinetic energy, K, and an over-prediction of the dissipation 
rate of K, Equation 3.9. Therefore, the results of Taylor [39] show 
relatively lower turbulence kinetic energy levels than the present 
computed results. However, his results agree well with the data of 
Bradley [9]. The two-dimensional model of Taylor [39] developed for 
flow over gentle terrain results in good agreement with the K/Kj^ profile 
for a set of three-dimensional, full-scale field data for a high hill. 
Figure 6.20, and in disagreement with the V/Vn profile, Figure 6.19. 
Therefore, more data are needed to verify the accuracy of Taylor's 
model when it is applied to the cases of flow over high hills. 

The present results show maximum values of turbulence kinetic 
energy occur at z = 0.3. Taylor's method [39], however, shows a 
constant turbulence layer due to his assumption of a constant shear 

layer in the wall region. His results [39] shown in Figure 6.20 were 

~ ~ 2 
obtained for the case z^ = 0.0005. Since K^ = (V*/C ) is used in the 

present study, K^ = 0.016 if C = 0.416 and V^ = 0.053, which corresponds 

to the case Zq = 0.0005. The present computed results show that the 

value of R can be as large as ten times its undisturbed value, Kq, 

if Zq = 0.0005 and D/H = 2.0, Figure 6.21. Therefore, the maximum 

value of K is 0.16 if K« = 0.016. The square root of K = 0.16 results 

105 



z/H 



8.0 
6.0 

4.0 



2.0 



1.0 

0.8 

0.6 
0,4 



0.2 



0.1 

0.08 

0.06 

0.04 



0.02 



0.01 



1.02 0.005 




2.00 0.005 




2.00 0.0005 




Taylor [39], z^ 


= 0.0005 


Bradley [9], z^ 


^ 0.005 




K/K, 



Figure 6.21 Distribution of dimensionless turbulence kinetic energy, 
K/Kq, on top of semiellipses with different z and D/H 
ratios 



106 



1n the value of the ratio /K/Vn = 0.4. Some evidence [5,7] exists 
that the value of cr^ /Vu is approximately 0.4 on top of surface 
obstacles. As discussed earlier, the value of Oy //K is 1.03 for flow 
over uniform, level terrain, Equation 6.3, and is 0.93 in the layer 
directly on top of Bradley's hill [9]. The value of /K/Vm =0.4 



obtained in the present computation is reasonable if the value of 
^Vt/^H = 0.4, which results in the ratio oy /v^ ~ 1.0. 



6.5 Phenomenon of Flow Separation 

As flows pass over a gradually sloping surface, fluid elements 
close to the surface gradually lose momentum due to surface friction. 
Flow separation occurs when the momentum of the fluid particle cannot 
overcome the adverse pressure gradient caused by the obstruction to the 
flow, i.e., the fluid elements flow in the direction of decreasing pres- 
sure. As contrasted to the gradually sloping surface, an abrupt or 
bluff-shaped obstacle creates a small separation on the front face due to 
the same phenomenon described above, i.e., the interaction of adverse 
pressure gradients and the viscous forces; however, a much larger 
downstream separation emanates from the sharp leading edge due to the 
inability of the flow to negotiate the change of surface configuration 
[1]. For flow over non-bluff bodies, e.g., flow over semiellipses, 
the mechanism of flow separation and reattachment is the same as that 
of flow over gradually sloping surface [1]. The flow parameters 
affecting the phenomenon of flow separation are therefore the scale 
of the viscous forces, the surface friction, and the geometry of the 

obstruction. In addition to these parameters, the magnitude of the 

107 



upwind turbulence level significantly influences the size of the 
separated flow regions (see Section 6,4), In the present study, 
the above mentioned parameters are scaled by the Reynolds number, 
Re, the surface roughness, Zq, the aspect ratio, D/H, and the 
dimensionless friction velocity, V^, respectively. The effect of 
these parameters on flow separation is discussed in this section. 
The obstacles considered are a rectangular block and different aspect 
ratio semielliptical cylinders. 

The characteristics of wind around a rectangular block have 
been investigated [11] by using the same numerical model as that used 
in this study. For all the cases studies in Reference [11], the 
presence of a block in the wind field always induces an upstream recir- 
culating flow region and a downstream separation bubble. In general, 
the upstream bubble is small and extends approximately 1 H upstream. 
The size of the downstream separation bubble, however, is significantly 
affected by the upstream turbulence level. In the present model, the 
undisturbed upstream turbulence kinetic energy, Kq, is proportional to 
(V^) . The increase of V^ causes a smaller downstream bubble. As 
mentioned in Section 6.3, an increase in z« of the obstacle tends 
to cause greater displacement of the flow approaching and passing 
over it. Therefore, a larger value of Zg induces a larger upstream 
separation bubble for flow over a semiellipse [19]. However, 
the computation of flow over a forward-facing step [32] shows a minimum 
upstream separation bubble at Zq = 0.05. In view of the fact that the 
quantities obtained in numerical computation show an approximate rather 
than an exact location of the flow separation and reattachment points, 

108 



the effect of Zq on the upstream flow separation should be resolved by 
experimental work. 

Summarily, the flow field upstream of the block is influenced 
by the surface friction, the pressure gradients, the viscous force, and 
the turbulence level of upstream wind. The downstream flow field, how- 
ever, is influenced by the upstream turbulence level and the shape of 
the obstacle. 

The cases of flow over semi elliptical configurations studied in 
this effort did not produce flow separation. The molecular viscosity, 
V, is assumed to be much smaller than the eddy viscosity, therefore, 

the influence of the molecular Reynolds number is not felt. For 

5 
D/H = 1.02, Zq = 0.001, the present computed results with Re = 10 have 

been compared with the case of Re = «>, and only insignificant effects 

2 
on the flow field occur. Assigning conditions Re = 10 , however, 

results in downstream flow separation. Figure 6.22. The flow reattach- 
ment distance, x , where x is measured from the lower corner of the 
ellipse (Figure 6.22) is approximately 2 H, which is much smaller than 
that for laminar flow over a circular cylinder (x = 10 H has been 
reported [42]). In view of the fact that the turbulence model described 
in Section 3.0 is developed under the assumption of high Reynolds num- 
bers, reliability of the results obtained by applying the model to low 
Re flows is uncertain. However, assuming an eddy viscosity of zero and 
solving the system of equations described herein, a reattachment distance, 
x^> which extends to 7 H for D/H = 1.0001, is found. This value may 
still be an under-prediction. A coarse grid size used in the computa- 
tion procedure can result in under-prediction of the downstream separa- 
tion bubble of a circular cylinder [42]. 

109 



z/H 



1.8 

1.5 



-1.2 

1.0 

0.8 

-0.6 

"""^"^^0.4 



-3 



-2 




•0.2 



0.01 



-0.08"~| ' Q'Q '^^^ 



x/H 



Figure 6.22 Streamline,^ , patterns over a semiellipse for D/H = 1.02, 



z^ = 0.001, and Re = 100 



Figure 6.23 shows the co contours corresponding to the case 
shown in Figure 6.22. The distribution of the oi is wery similar to 
the results of a laminar flow calculation for flow over a circular 
cylinder [42]. This implies that the turbulence model developed under 
the assumption of high Reynolds number can be applied to low Reynolds 
number flows by assuming a \/ery small value for turbulence kinetic 
energy, which leads to a small value of v. . It is noteworthy, however, 
that K should not be zero; otherwise a singularity occurs at the term 
E, of Equation 3.14. 

For laminar flow, increasing Re increases the size of the flow 
separation bubble [42]. An increase in Re can be generated by increas- 
ing the upstream velocity and therefore increasing the momentum of the 
approaching flow. This increased inertia force convects the separation 
bubble downstream. In the case of turbulent flow, an increase in Vn 
causes a decrease in V^/Vm» if V^ is kept constant. The downstream 
separation bubble is therefore increased. Figure 6.17. The mechanism 
governing the size of the flow separation bubble, however, is different 
in the laminar flow case and the turbulent flow case. In low Re flows, 
the size of the reversed flow region is controlled by viscous diffusion 
of momentum rather than by turbulent diffusion. The shear layer shed 
behind the body spreads more rapidly and reattaches sooner in turbulent 
flow. 



Ill 



z/H 



-0.3 



^^^'^ZT- — -^,. -oTs"^^^-- 



-0.5 




-0.5 



-0.4 



-3 



x/H 



Figure 6.23 Vorticity, w, patterns over a semiellipse for D/H = 1.02, 



Zq = 0.001, and Re = 100 



7.0 CONCLUSIONS 

The values of the constant coefficients developed in the present 
study for use in the K-L turbulence model are applied to the solution of 
the problems of atmospheric flow over two-dimensional obstacles. Good 
agreement between the numerically computed results and the experimental 
results from both atmospheric boundary layer wind tunnels and full-scale 
field tests are obtained. 

The results of studying the influence of the flow parameters 
on the flow field for the range of parameters investigated in this 
study lead to the following conclusions: 

1. The character of the flow field is significantly affected 
by the surface roughness length scale, z^/H; the upwind 
turbulence level, V^/Vn; and the geometries of the obstacles. 

2. The drag force exerted on the flow field by the solid sur- 
face is larger for increased z^/H and thus more effectively 
retards the flow near the surface, causing the streamlines 
to be further displaced from the obstacle. 

3. The turbulence intensity of the fluid motion is increased in 
the wall region for increased Zq/H. In turn, the momentum 
of the fluid elements in the wall region is reduced due to 
the energy being extracted from the mean flow field and 
converted to turbulence energy. 

4. The thickness of the layer of atmosphere next to the ground 
where the flow is significantly affected by the local surface 

113 



roughness depends on the distribution of the roughness and 
its magnitude. This thickness is generally less than one 
characteristic body height. 

5. Flow separation generally commences from the upstream corner 
of the bluff body. For flow over a convex surface, the flow 
separates at a position downstream of the crest of the 
obstacle if Re is low. For high Re, the flow is not likely 
to separate, particularly if the upstream turbulence is high. 

6, An increase in upstream turbulence level reduces the size 
of the downstream flow separation bubble. This effect is 
attributed to enhanced turbulent mixing which transports 
momentum much more quickly from the upper layer to the 
lower layer. 

The results of the present study indicate that the turbulent 
flow field upstream of an obstacle is mainly affected by the surface 
friction and the pressure gradient. Downstream of the obstacle, however, 
the flow is mainly affected by the upstream turbulence level. Addition- 
ally, the mean velocity field and the turbulence field are significantly 
affected by the geometry of the obstacle. For low Reynolds number 
flows, the size of the downstream reversed flow region is controlled 
by viscous diffusion of momentum rather than by turbulent diffusion. 
Therefore, the shear layer shed behind the obstacle spreads more 
rapidly in turbulent flow, resulting in a smaller downstream separation 
bubble. 



114 



BIBLIOGRAPHY 



1. Frost, W., and C. F. Shieh. "Guidelines for Siting WECS Relative 
to Small -Scale Terrain Features." Report prepared under Department 
of Energy Contract No. EY-76-C-062443 by FWG Associates, Inc., 
Tullahoma, Tennessee, September, 1979. 

2. Gosman, A. D., W. M. Pun, A. K. Runchal , D. B. Spalding, and 
M. Wolfshtein. Heat and Mass Transfer in Recirculating Flows . 
London and New York: Academic Press, 1969. 

3. Launder, B. D., and D. B. Spalding. "Two-Equation Models of Tur- 
bulence," Mathematical Models of Turbulence . London and New York: 
Academic Press, 1972. 

4. Frost, W., B. H. Long, and R. E. Turner. "Engineering Handbook on 
the Atmospheric Environmental Guidelines for Use in Wind Turbine 
Generator Development," National Aeronautics and Space Administra- 
tion Technical Report 1359, George C. Marshall Space Flight Center, 
Alabama, December, 1978. 

5. Woo, H. G. C, J. A. Peterka, and J. E. Cermak. "Wind-Tunnel 
Measurements in the Wakes of Structures," National Aeronautics 
and Space Administration CR 2806, George C. Marshall Space Flight 
Center, Alabama, March, 1977. 

6. Bowen, A. J., and D. Lindley. "A Wind Tunnel Investigation of the 
Wind Speed and Turbulence Characteristics Close to the Ground over 
Various Escarpment Shapes," Boundary-Layer M eteorology , 

12(No. 3):259-271, October, 1977. 

7. Meroney, R. N., V. A. Sandborn, R. J. B. Bouwmeester, H. C. Chien, 
and M. Rider. "Sites for Wind Power Installations: Physical 
Modeling of the Influence of Hills, Ridges and Complex Terrain on 
Wind Speed and Turbulence," Final report prepared under Department 
of Energy Contract No. EY-76-S-06-2438, AOOl , by Colorado State 
University, Fort Collins, Colorado, June, 1978.^ 

8. Frost, W., and A. M. Shahabi. "A Field Study of Wind over a 
Simulated Block Building," National Aeronautics and Space Adminis- 
tration CR 2804, George C. Marshall Space Flight Center, March, 
1977. 

9. Bradley, E. G. "An Experimental Study of the Profiles of Wind 
Speed, Shearing Stress and Turbulence at the Crest of a Large 
Hill," Unpublished report by Division of Environmental Mechanics, 
Commonwealth Scientific and Industrial Research Organization, 
Canberra City, A.C.T., Australia, 1978. 



115 



10. Frost, W. "Review of Data and Prediction Techniques for Wind 
Profiles Around Manmade Surface Obstructions," AGARD Conference 
Proceedings No. 140 on Flight in Turbule nce. London: Technical 
Editing and Reproduction, Ltd., 1973. Pp.' 4.1-4.18. 

11. Shieh, C. F., W. Frost, and J, Bitte. "Neutrally Stable Atmo- 
spheric Flow over a Two-Dimensional Rectangular Block," National 
Aeronautics and Space Administration CR 2926, George C. Marshall 
Space Flight Center, Alabama, November, 1977. 

12. Tani, N. "On the Wind Tunnel Test of the Model Shelter-hedge," 
Bulletin of the National Institute of Agricultural Sciences , 
A(No. 6):l-80, March, 1958. 

13. Plate, E. J. "Aerodynamics of Shelterbel ts," Agricultural 
Meteorology , 8(No. 3):203-222, May, 1971. 

14. Chang, S. C. "Velocity Distributions in the Separated Flow 
behind a Wedge Shaped Model Hill," Report prepared under United 
States Army Research Grant No. DA-AMC-28-043-G30 by Colorado 
State University, Fort Collins, Colorado, March, 1966. 

15. Counihan, J., J. C. R. Hunt, and P. S. Jackson. "Wake Behind 
Two-Dimensional Surface Obstacles in Turbulent Boundary Layers," 
Journal of Fluid Mechanics , 64(Part 3):529-563, July, 1974. 

16. Taylor, P. A. "Numerical Studies of Neutrally Stratified Planetary 
Boundary-Layer Flow Above Gentle Topography," Boundary-Layer 
Meteorology , 12(No. l):37-60, August, 1977. 

17. Jackson, P. S., and J. C. R. Hunt. "Turbulent Wind Flow over a 
Low Hill," Quarterly Journal of the Royal Meteorological Society , 
101:929-955, 1975. 

18. Deaves, D. M. "Wind Over Hills: A Numerical Approach," Journal 
of Industrial Aerodynamics , l(No. 4): 371 -391, August, 1976. 

19. Frost, W., J. R. Maus, and G. H. Fichtl. "A Boundary Analysis of 
Atmospheric Motion over a Semi-Elliptical Surface Obstruction," 
Boundary-Layer Meteorology , 7{No. 2):165-184, October, 1974. 

20. Lin, C. L., and S. C. Lee. "Transient State Analysis on Applica- 
tion of Computers to Fluid Dynamic Analysis and Design," 
Unpublished report prepared by Polytechnic Institute of Brooklyn, 
Brooklyn, New York, January, 1973. 

21. Haussling, H. J. "Viscous Flows of Stably Stratified Fluids over 
Barriers," Journal of the Atmospheric Sciences , 34:581-602, April, 
1977. 



116 



22. Guven, 0,, V. C. Patel , and C. Farell . "A Model for High- 
Reynolds-Number Flow Past Rough-Walled Circular Cylinders," 
Transactions of the ASME Journal of Fluid Engineering , 

99 (No. 3): 486-494. September. 1977. 

23. Bradshaw, P., and F. Y. F. Wong. "The Reattachment and Relaxation 
of a Turbulent Shear Layer," Journal of Fluid Mechanics , 
52(Part 1):113-135, March, 1972. 

24. Towsend, A. A. "Self-Preserving Flow Inside a Turbulent Boundary 
Layer," Journal of Fluid Mechanics , 22(Part 4):773-797, August, 
1965. 

25. Vincent, J. H. "Model Experiments on the Nature of Air Pollution 
Transport Near Buildings," Atmospheric Environment , 11 (No. 8): 
765-774, August, 1977. 

26. Shieh, C. F., and W. Frost. "Applications of a Numerical Model 
to WECS Siting Relative to Two-Dimensional Terrain Features." 
Paper presented at the Fifth International Conference on Wind 
Energy, Fort Collins, Colorado, July 8-14, 1979. 

27. Mulhearn, P. J., and E. F. Bradley. "Secondary Flows in the Lee 
of Porous Shel terbelts," Boundary-Layer Meteorology , 

12(No. l):75-92, August, 1977. 

28. Hinze, J. 0. Turbulence . Second edition. New York: McGraw- 
Hill, Inc., 1975. 

29. Rodi , W., and U. Karlsruhe. "Turbulence Models for Environmental 
Problems," Prediction Methods for Turbulent Flows , von Karman 
Institute for Fluid Dynamics Lecture Series 1979-2, Rhode Saint 
Genese, Belgium, 1979. Pp. 1-89. 

30. Mellor, 6. L., and H. J. Herring. "A Survey of the Mean Turbulent 
Field Closure Models," AIAA Journal , n(No. 5):590-599, May, 1973. 

31. Bradshaw, P. "The Understanding and Prediction of Turbulent 
Flows," Aeronautical Journal , 76:403-418, July, 1972. 

32. Frost, W., J. Bitte, and C. F. Shieh. "Analysis of Neutrally 
Stable Atmospheric Flow over a Two-Dimensional Forward-Facing 
Step," AIAA Journal , 18(No. l):32-38, January, 1980, 

33. Lewellen, W. S. "Use of Invariant Modeling," Handbook of 
Turbulence , W. Frost and T. H. Moulden, editors. Vol. 1. 
New York: Plenum Press, 1977. Pp. 237-280. 

34. Reynolds, A. J. "The Variation of Turbulent Prandtl and Schmidt 
Numbers in Wakes and Jets," International Journal of Heat Mass 
Transfer , 19:757-764, July, T9^76. 

117 



35. Zeman, 0., and H. Tennekes. "A Self-Contained Model for the 
Pressure Terms in the Turbulent Stress Equations of the Neutral 
Atmospheric Boundary Layer," Journal of the Atmospheric Sciences , 
32(No. 9):1808-1813, September, 1975. 

36. Frost, W., W. L. Harper, and G. H. Fichtl. "Analysis of Atmo- 
spheric Flow over a Surface Protrusion Using the Turbulence 
Kinetic Energy Equation," Boundary-Layer Meteorology , 

8(Nos. 3/4):401-417, April/June, 1975. 

37. Roche, P. J. Computational Fluid Dynamics . Revised printing. 
Albuquerque: Hermosa Publishers, 1976. 

38. LeFeuvre, R. F. "The Prediction of 2-Dimensional Recirculating 
Flows Using a Simple Finite-Difference Grid for Non-rectangular 
Flow Fields," Computers and Fluids , 6 (No. 4): 203-21 8, December, 
1978. 

39. Taylor, P. A. "Some Numerical Studies of Surface Boundary-Layer 
Flow above Gentle Topography," Boundary-Layer Meteorology , 
n(No. 4):440-465, July, 1977. 

40. Panofsky, H., M. Vilardo, and R. Lipschutg. "Terrain Effects on 
Wind Fluctuations," Fifth International Conference on Wind 
Engineering . Vol. 1. Fort Collins: Colorado State University, 
1979. 

41. Raine, J. K. , and D. C. Stevenson. "Wind Protection by Model 
Fences in a Simulated Atmospheric Boundary Layer," Journal of 
Industrial Aerodyamics , 2CNo. 2):159-180, June, 1977. 

42. Tuann, S., and M. D. Olson. "Numerical Studies of the Flow 
around a Circular Cylinder by a Finite Element Method," Computers 
and Fluids , 6(No. 4):219-240, December, 1978. 



118 



APPENDIX 

EXPRESSION OF S^ IN ORTHOGONAL CURVILINEAR 
COORDINATE SYSTEMS 

Consider an orthogonal curvilinear system (un,Up) rotated an 
angle e with respect to the Cartesian system (x,z) as shown in 
Figure A. 1(a), the unit vectors e and e in (x,z) system can be 

X z 

expressed as: 



e = cos e e. - sin e e 2 



e = sin e e-| + co.- 9 e^ 



(A.l) 



where e. and e^ are the unit vectors in the (u-jjUp) system. 

The distances between two nearby points in the (u^jU,,) system, 
i.e., 6i and dn in Figure A. 1(b), can be written as 



d£ = h-jdu-i 



dn = h^dup 



(A. 2) 



where h, and h^ are metric coefficients in (u-|,U2) system. The deriva- 
tive of a quantity (^ with respect to u, can be expressed as: 



3A_ 



3U. 



3xJ dZ dU-, 



3j) 

dz 



dz d_l 
95, 3U 



= h. 



M. 



3cJ) 



cos e T^ + sin e .^ 

dX dZ 



(A. 3) 



Similarly 



119 




(a) Two-dimensional orthogonal curvilinear coordinate 
systems 




(b) A small surface element in (u,,Up) system 



Figure A.l Schematic diagram of the relationship between general 
orthogonal curvilinear coordinate system and Cartesian 
coordinate system 



120 



34) 



3U 



2Ju. 



= ho -sin 6^4+ cos e ^^ 



{A. 4) 



Multiplying Equations A. 3 and A. 4 with h^ cos 6 and h^ sin 6, respec- 
tively, and then subtracting, one obtains 



9(t) ^ cos 9 9(j) _ sin 9 9(t) 



8X h-, 9Un 



hp aup 



(A.5) 



In turn, multiply Equations A. 3 and A. 4 with h^ sin 9 and h^ cos e, 
respectively, and adding gives 



54) ^ si_n_e 34) _|_ cos 9 34) 
3z h-, sun hp 9u- 



(A.6) 



Substituting Equations A.l, A.5 and A. 6 into Equation 4.31, the expres- 
sion for V and V in the {Ui,U2) system becomes: 



V4) 



_]- 94) + _2 94) 
h-, 3u-| hp su^ 



^* h2 3U2 h^ 3u^ 



(A. 7) 



The velocity vector V appearing in Equation 4.30 can be written 



as: 



-> . .. -> 



V -~ V^e^ + V^e^ 



(A. 8) 



Expanding Equation 4.30 with the aid of Equations A.l, A. 7 and A. 8 



121 






^X X ^ ^x X 



aup 9U-] 



au-i 3u2 



3p^ 3V, ^ Sy, 3V^ 



SU^ 9U-| 



8U^ 8U2 



(A.9) 



where 



8u . 9u 

COS e e sin e e 



^x " ~h 



1 ^^1 



h2 3U2 



^ 3y r. 9y 
Sin 9 e cos 9 e 

^z h-, 9U-, h2 9U2 



V = Vt cos 6 - V^ sin 
X 1 2 



V = Vt sin 9 + V^ cos 
z I 2 



(A. 10) 



Substituting Equation A. 10 into A.9 and arranging gives Equation 4.32 in 
the text: 



h^h^f 



9V, 9y, 9V. 3y, 3V2 9^2 ^^2 9y2 



9Ui 9U2 



9U2 9U-J 



9u-| 9U2 



9U2 9u-[ 



39 ^^2 



38 ^^2 



9U-. dU^ 



9U2 9Ui 



+ y^ 



36 ^^1 



39 ^^1 



3U2 9Ui 



9Ui 9U2 



+ V. 



99 ^^2 



99 '^2 



9U2 9u-| 



9Ui 9U2 



+ V, 



96 ^'^l 



39 ^^1 



9Ui 9U2 



9U2 9U. 



(A. 11) 



where 



122 



Pi - 









(A. 12) 



123 



1. REPORT NO. 

NASA CR-3430 



2, GOVERNMENT ACCESSION NO. 



4. TITLE AND SUBTITLE 

NumeJ'ical Solutions of Atmospheric Flow Over 
Semielliptical Simulated Hills 



7. AUTHOR(S) 

Chih Fang Shieh and Walter Frost 



9. PERFORMING ORGANIZATION NAME AND ADDRESS 

Atmospheric Science Division 

The University of Tennessee Space Institute 

Tullahoma, Tennessee 



12. SPONSORING AGENCY NAME AND ADDRESS 



National Aeronautics and Space Administration 
Washington, D.C. 20546 



3. RECIPIENT'S CATALOG NO. 



5. REPORT DATE 

June 1981 



6. PERFORMING ORGANIZATION CODE 



8. PERFORMING ORGANIZATION REPORT # 



10. WORK UNIT NO. 

M-348 



1 1 CONTRACT OR GRANT NO. 

NAS8-32692 



13. TYPE OF REPOR'i' & PERIOD COVERED 

Contractor 
{Interim Report) 



l-l. SPONSORING AGENCY CODE 

989 13 ZZ 



15. SUPPLEMENTARY NOTES 

Marshall Technical Monitor: Dennis W. Camp 



16, 



ABSTRACT 

Analysis of atmospheric motion over obstacles on plane surfaces is carried out 
in the present study to compute simulated wind fields over terrain features. 
Emphasis is on a semielliptical, two-dimensional geometry. Numerical simulation of 
flow over rectangular geometries is also discussed. The numerical procedure utilizes 
a two-equation turbulence model, and the selection of the necessary constant 
coefficients in the model is considered an important part of the present study. 

In the present approach the partial differential equations for the vorticity, 
stream function, turbulence kinetic energy, and turbulence length scale are solved 
by a finite-difference technique. The numerical solutions have been compared with 
available experimental data; agreement is good. 

It is found that the mechanism of flow separation induced by a semiellipse is 
the same as in the case of flow over a gradually sloping surface for which the flow 
separation is caused by the interaction between the viscous force, the pressure force, 
and the turbulence level. For flow over bluff bodies, e.g., solid fences, a downstrean 
recirculation bubble is created due to the inability of the flow to negotiate with 
the abrupt change of the surface shape. Increasing the aspect ratio and/or increasing 
the turbulence level results in flow reattachment close behind the obstacle. 



17. KEY WORDS 

Low-level flow 
Turbul ence 
Turbul ence models 
Vorticity 



19. SECURITY CLASSIF. (of thli roportl 

Unclassified 



18. DISTRIBUTION STATEMENT 

Unclassified - Unlimited 



Subject Category 47 



20. SECURITY CLASSIF. (of this page) 

Unclassified 




For sale by National Technical Information Service, Springfield, Virginia 221«1 



NASA-Langley, 1981