# 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