Numerical calculation of turbulent flow around two-dimensional hills

Journal o f Wind Engineering and Industrial Aerodynamics, 21 (1985) 307--321

307

Elsevier Science Publishers B.V., Amsterdam -- Printed in The Netherlands

NUMERICAL CALCULATION OF T U R B U L E N T FLOW A R O U N D TWO-DIMENSIONAL HILLS

G.C. BERGELES

Laboratory o f Aerodynamics, National Technical University o f Athens, Athens (Greece) (Received January 22, 1985; accepted in revised form July 5, 1985)

Summary Predictions are presented of the two-dimensional turbulent flow around a hill; the two-dimensional Reynolds equations are written in an orthogonal curvilinear coordinate system and transformed to finite difference form after discretisation in physical space. Turbulence is simulated by the k--e model of turbulence; the first part of the paper compares predictions with available laboratory data; in the second part the results of the developed method are compared with results of other numerical methods; therefore, the validity of the method, its flexibility to treat any hill shape and its reliability compared to other methods is assessed.

Notation

he

coefficient in finite difference equation linking the central node, P, with adjacent nodes, i (E, W, N and S), for the dependent variable, ¢ (¢ = u, v, k, e) ratio of the metric coefficients (l~/l~) of the curvilinear orthogonal coordinate system height of the hill

k

turbulent kinetic energy

l l~ /l~

length scale of turbulence for flat plate metric coefficients for the ~ and ~ directions of the curvilinear orthogonal coordinate system eharacteristic length of the hill pressure source term for the dependent variable, ¢ (¢ = u, v, k, e) mean velocity c o m p o n e n t in the ~ direction u velocity fluctuations mean velocity c o m p o n e n t in the r~ direction Reynolds shear stress in ~ direction Cartesian coordinates

A~ (A~, A~V , etc)

L P

S¢ u u' u - - p U'V' X, y

0167-6105/85/$03.30

u z + v'2 + w'Z ) 2

© 1985 Elsevier Science Publishers B.V.

308

Greek letters

77 ÷

P

Indices E,W, N, S ¢ P P t oo

rate of dissipation of turbulent kinetic energy coordinate distance in the orthogonal coordinate system (~,17) dimensionless distance from the ground (77÷ = pu~/~) viscosity coordinate distance in the orthogonal coordinate system (~ ,77) density

grid nodes adjacent to the central node P dependent variable (~ = u, v, k, e) nearest to the ground grid point central grid node turbulent flow free stream conditions

Introduction The knowledge of the flow development around hills, or generally in complex terrains, is of significant importance to various branches of engineering such as wind and environmental engineering, pollution, etc. The numerical study of flows, particularly around hills, has progressed in the last few years. This study presents a numerical model for the calculation of such flows. Various studies have been reported which differ in the level of assumptions regarding the flow field or the physical domain. A basic contribution to the subject is the work of Jackson and Hunt [1] (henceforth called JH). Their model applies an analytical perturbation theory valid for low hills without flow reversal. Under these limitations it can give quantitative estimates of the effects of the hill on the wind. Other studies have applied the JH theory in predicting the flow around various shapes of escarpment. Jensen [2] has applied the JH model to the flow past a ramp-like escarpment for which wind-tunnel measurements were available. He reported that the maximum velocity perturbation was predicted successfully, but the vertical and downwind variations of the velocity perturbations were not calculated accurately. Sykes [3] has also developed an analytical theory based on matched asymptotic expansion (along the lines of the triple-deck theory) for a turbulent boundary layer encountering a short two-dimensional hump w i t h o u t flow reversal. The theory concludes that the leading-order velocity perturbations are those of the linearised potential flow solution. Both studies [ 1] and [ 3 ] , despite their limitations, shed light on the structure of the interaction of the viscous layer near the hill with the inviscid region. Numerical solutions of the governing equations have much greater poten-

309 tial than analytical methods. Relevant studies here are the works of Taylor and Gent [4] and Taylor [5] for two-dimensional flow above gentle topography. Again, in these models flow reversal in the flow field is not allowed, and boundary layer t y p e approximations are assumed. In the work of Taylor [5], the governing equations are expressed in Cartesian form and a coordinate transformation is used to transform the hill shape to a straight line, the equations being solved on the transformed plane. The method uses empirical modifications to account for the non-orthogonality of the vertical coordinate to the ground surface. Deaves [6] presents a numerical solution of the governing equations in Cartesian coordinates and in primitive form. The method, however, cannot properly treat the shape of the hill (within the frame of the Cartesian coordinates), and has problems in imposing the correct boundary conditions on the hill surface. Problems also arise in the refinement of the Cartesian grid, particularly near the hill contour. The method gives quantitatively good results for low hills, b u t as Deaves [6] points out, examination of the assumptions of the analysis imply that the results obtained for particularly large or steep hills should be treated with caution. In a later publication, Deaves [7] solved the governing equations in stream function and vorticity formulation, using a transformed grid which overcomes most of the problems associated with the Cartesian grid. Again, flow reversal cannot be accounted for by the model. In a recent review article by Walmsley [8] it is stated that the numerical models of Taylor and Gent and of Deaves are probably among the best available for two-dimensional, neutrally stratified, turbulent, unseparated flows above gentle topography. Effects which lead to flow reversal, however, may be difficult to handle without significant modifications to the original models. This paper presents a numerical model for environmental flows around two-dimensional hills which has the generality previous models lack, i.e., it can be applied for any hill shape, for large or steep hills with or without flow reversal. The problems of introducing into the numerical solution of the governing equations the hill shape, correct boundary conditions along the contour of the hill making allowance for the appearance of reverse flow caused by the presence of the hill, can be alleviated by applying numerical solution techniques to the governing elliptical equations expressed in a curvflinear orthogonal coordinate system, in which the contour of the hill coincides with a boundary coordinate line. Also, the use of such a coordinate system, in which the velocity components are locally aligned to the grid, minimises errors due to numerical diffusion. Previous work along these lines has not been conducted for environmental flows. However, for internal flows there has been progress in the last few years (Antonopoulos [9], Gosman and Rapley [10]). In these two studies the SIMPLE numerical algorithm of Patankar and Spalding [11] is used for the solution of Reynolds equations expressed in a curvilinear orthogonal coordinate system.

310 Its reliability and success has been assessed, along with the generality of the method. Particular to each problem is the proper generation of the orthogonal curvilinear coordinate system fitted to the solution domain of the equations. The next paragraphs present the mathematical foundation of the problem, the m e t h o d of the grid generation, and the solution algorithm. The study is concluded with the presentation and discussion of the results.

Mathematical foundation and numerical solution Equations and boundary conditions The time-averaged Reynolds equations which express the turbulent flow field around a two-dimensional hill can be written in an orthogonal curvilinear coordinate system as a

1 --

a

1

--

(puln

1~l,~ a~

) +--

~

ltl~ 37?

) = 0

(pvl~

(1)

- - n -an - (pvl~¢) -l~l - -n a~ ( p u l n ~ ) + ltl

#~ ~

ltln

where ¢ = u, v and Sv source terms, and It, In are the metric coefficients of the coordinate system. B

:

'

:

~

'

' .

.

.

• • • '

~ " .

I

i

"

I ~ - F - - 7

:

:

.

i

'

i

~-

~ - - ~ - F - q ~ [ " '

"

~

~

:

"

~ .......

i

=

]

'

'

I

T

~

,

'-~,~

~

,

,

•

......

C

i

'

:

'

+

' ; I ,

! '

,

"

'

;

I

', ~ .

' ' '

" '

I

:,

i

/

/ I • | . [ R T / P / - h ~ k ~ H h ~ t l ; [ ~4.---ff-~" , ,

• : "

.....

, '

, "

, : 'l

~

:

• [ '

'

/

~ "

~

'

,

•

i ,

i

:

I

~

: '

,

,

:

'

I t

.

/

,

.

!

,

',

'

,

' '

.

~

i

, /

' Ii

I i

!

I,

I

[

' ,

" [

:

, •

'

!

[

'

[

,

!

i

~l

,

i

', t

~,

m I-

,

,

'

!

~

!

:

+-.~-_

; .~

= i ,

' I

I

'

,

'

'

I

3h

!

[

,

!

;

'

,

Fig. 1. Orthogonal curvilinear coordinates for the RTP hill.

[

l ,

:, ;

l '

'

,

,4-4- ~

A 5h

"

i

i

, , • i, :~ L~_~_:_~--~ ' ~ - r 4,1 -I:~, I LT ~t-~, ~ 4I -I~ - ~ ;

' - - 1 - - - ~ - F :

: ......

,

~-4-4-~t$ 4-~-, * ~--~-4 ....... ,-~-4--+-~ t +

'

I

- 7 - T /

•

/ :

;

.....

.

4..... ~-~-4~J ~ - 1 - 4 - 4

/

-~-~-'T: ,

.

•

'

" r " - ' - 7 " ~ T .

t

,~ r I

'

~_L~

I i

I

' ,

i ' / I

:

'. ,

L / ,

r-l i

I I

311

The orthogonal coordinate system fitted to the shape of the hill is shown in Fig. 1 where the 77 = 0 line coincides with the lower boundary of the integration domain. The system of equations retains the m o m e n t u m fluxes in the streamwise ~ direction and renders the equation elliptical. Therefore, in principle, flow reversal in the flow field can be properly treated. Equation (2), due to its ellipticity, requires boundary conditions around the boundary as shown in Fig. 1. Along the lower line, the no-slip boundary condition on velocity is applied, whilst at the inlet plane a logarithmic velocity profile is applied. The top plane is taken far away, so that the flow is not influenced by the presence of the hill and free stream conditions can be imposed. The downstream plane is also taken far downstream from the hill, so that the flow there is parabolic. In fact, the locations of the last two planes are selected after numerical experimentation so that their particular location has no influence on the results. Turbulence

closure

The Reynolds stresses ( - p u'v') which appear in the time-averaged Reynolds equations have to be modelled with the help of known or calculable quantities, in order to render the set of equations (1) and (2) complete. The introduction of the Boussinesq hypothesis for turbulent viscosity leads to the two-equation turbulence model of Jones and Launder [12] for the turbulent kinetic energy (k) and its rate of dissipation (e). Viscosity is calculated from Pt = p k 2 / e where k and e are obtained from the solution of two transport equations of the form of eqn. (2) with appropriate values of the source term S , . For details see ref. 9. Equations for k and e are also elliptical, and, therefore, boundary conditions have to be specified on the contour of the integration domain. They are shown analytically in Table 1. Particular attention is given to the near wall region where steep gradients of the dependent variables exist and, therefore, a fine numerical grid is needed for adequate space resolution in this area. This need is avoided by the use of a special wall treatment, TABLE1 Boundary conditions u

v

k

e

l n ( E y *) - - u~ 0.41

v = 0

Kolrnogorof's k distribution f o r flat p l a t e

k3/2/1

BC

u = u~

v = O

CD

Ou/~x

AD

u = O

u

AB

= 0

~v/~x v = O

= 0

~k/~y

= O

~k/~x

= 0

~k/by

=0

w h e r e l = l e n g t h scale f o r flat p l a t e

ae/ax = 0 ep = hp312/lp

p is t h e grid n o d e nearest to the surface

312

using the so-called wall functions, with the help of which the viscous sublayer is bridged to the fully turbulent region. Details of the philosophy and applications of the wall functions can be found in the work of Launder and Spalding [ 13]. The effect of surface roughness on the near-wall local equilibrium of the flow is introduced via the constant in the logarithmic law of the wall. The use of an orthogonal coordinate system bypasses the need for empirical modifications which previous workers [3] were obliged to introduce into the mixing length distribution. Also, the effects on turbulence of stream line curvature are included in the k--e model, thus giving it a higher degree of universality compared to the zero equation model of previous workers.

The numerical grid For the numerical solution of the governing partial differential equations, as expressed by eqns. (1) and (2), an orthogonal curvilinear coordinate system fitted to the solution domain must be generated. The boundaryfitted orthogonal curvilinear coordinate system is obtained by solving a set of Laplace equations (3 and 4) with Neumann boundary conditions, so as to ensure orthogonality of the coordinate lines to the boundaries. ~2x 1 a2x + =0 a~2 h 2 8 ~

(3)

b2y

1 ~2y + ~)~2 h2 ~C~2 = 0

(4)

where h is the ratio of the scale factors in the two orthogonal directions and 7. Figure 2 shows typical orthogonal curvilinear coordinate systems obtained after the solution of the system of Laplace equations for various hills.

3

Fig. 2. Orthogonal curvilinear coordinates for the JH hills.

~r

,"

313

Numerical solution o f the equations All the equations can be cast into the general form (eqn. 2) with appropriate values of the source term depending on the variable under consideration. Following the practice of the SIMPLE algorithm of Patankar and Spalding [11], the dependent scalar variables are stored at the centre of the cells whilst the corresponding velocities are specified at the cell interfaces. The equation is integrated over the finite volume surrounding the control point, and by application of the Gauss and mean value theorems, the equation is transformed into finite difference form (A~ - S2~p) ~p = A~b E + Aw~b ~ w + A ~ ¢ N + A ~ s + S~p

(5)

where A~ = A~ + A~v + A~ + A ~ and coefficients A 7 link the dependent variables on adjacent grid nodes and contain convection and diffusion fluxes through corresponding interfaces. Equations like eqn. (5) are written for all the dependent variables, u, v, k, e. The method of solution is based on the observation that if the pressure field were known then the m o m e n t u m equations could be solved. Consequently, with the known velocity field the rest of the equations (k, e) could be solved in an iterative mode in order to get the turbulent viscosity. Now since the velocity field is calculated with a guessed pressure field the continuity equation is not satisfied and mass imbalances occur. These mass imbalances are used to correct the pressures via a pressure correction equation. Finally, with the corrected pressure field, the iteration cycle is repeated and a new velocity field, which is more accurate, is obtained from the m o m e n t u m equations. This m e t h o d has been extensively used by various workers [ 14, 15].

Grid independency tests The results which are presented here have been obtained with a numerical grid of 60 X 40 in the x,y directions, respectively. Preliminary results have been also obtained with coarser grids of 30 X 15 and 45 × 25, which show that the results presented here are grid independent. Special care was taken to achieve grid-independent and converged results since previous studies [16, 17] have shown that discrepancies observed in predictions with measurements by previous workers were due to the dependency of the results on the grid, and not to the inadequacy of the turbulence model. Figure 3 shows the rate of convergence of the flow around the hill of Fig. 2. The fineness of the grid affects the rate of convergence, and the mass imbalances in the field are reduced to 10 -4 of the mass inflow to the solution domain in approximately 400 iterations, depending, of course, on the degree of ellipticity of the flow field. The calculations have been performed on a PRIME 450 minicomputer.

314 10 -1

lo-2 _m

E_

\ lo-3

\

\}_ i h~lO

1(~4

1 , 100

h- 20

_ ~'.~

I Iterah°ns 200

.. 1 300

\\

,

I ' ~ 400

5OO

Fig. 3. Rate of convergence for the JH hills. Presentation and discussion o f t h e results

In order to test the accuracy of the code developed there is a need for reliable experimental data for flows around two-dimensional hills. Such data can only comprise laboratory measurements in wind tunnels where controlled and well-defined boundary conditions can be imposed. A systematic study of such a flow has been conducted by the environmental laboratory at Research Triangle Park during the joint Soviet--American hill study program, the catalogued results of which have been published [18]. The study has examined hills of low, m e d i u m and high slope, the last meaning that there is flow reversal on the downstream side of the hill. The prediction of the last flow is a serious test for the present code, since one of the features of the method is its capability to predict flow reversal. Figure I shows the hill-fitted orthogonal curvilinear coordinate for a hill height-to-length ratio of 3 (only a few of the grid lines are plotted). The predicted velocity distribution at various locations upstream and downstream of the hill and its comparison with the measurements of ref. 18 is shown in Fig. 4a. O n the top of the hillthe speed increase is predicted accurately within few percent, as is the far-downstream recovery of the flow. The measurements at the downstream base corner of the hillin ref. 18 were taken with various velocity probes, such as hot-wires and pitot tubes. The presence of the flow reversal there makes hot-wire measurements unreliable, whilst the magnitude of the velocity (0.5 m s-l), and the high degree of turbulence make the experimental error of the pitot measurements

315

rather significant. Figure 4b shows a comparison between the predicted and measured velocity distributions. The general level of agreement could initially be regarded as satisfactory considering the uncertainty of the measurements in the region. A further test of the code is shown in Fig. 5, where predictions are compared with the wind-tunnel measurements of Bowen and Lindley [19] over an escarpment of slope 2:1. The flow field near the escarpment which is the most difficult to predict compares favorably with the measurements. Deaves [6, 7] has solved the parabolised Naviet--Stokes equations for flow around a hill described by Jackson and Hunt [1] numerically. Apart from the problems which the use of Cartesian coordinates creates in imposing correct boundary conditions, and also the errors due to numerical diffusion, the code developed by Deaves cannot

15(

Measurements Ref. 18

A,I t

•

.

y(mm',

" /~jd

10o

/

I

/

" // / //

/4

a

/i

i/ I a/ • i I A/"

01

0.2

0.4 U/Uco

i

."

•I

50

o -0

:

/oi

I

o

o .I

0.6

Fig. 4. Velocity predictions for the RTP hill.

0.~

10

316 1000

D

100

o D

y(mm)

[]

0

9

-

A

0

°f

:

DO 10

D0

/- -

A

\~ D0

A

\

Measurements

Ref. 18

INSTRUMENT

LI~ m/s

A

HOT --WIRE

B

PITOT

4 4

0

PITOT

8

\ - 0 1.5

.

0.

.

015 .

0 30 045 U/U~

0 '6 0

0 ' 75

090

1,~5 '~

Fig. 4 (continued).

treat separated flow properly. Jackson and Hunt have applied their theory for two hills of a shape described by the equation y = hi[1 + (x/L) 2]

The curvilinear orthogonal coordinate system created for the study of hills of this family (one low and one high with h = 10, h = 20 and L = 50) are shown in Fig. 2. The predicted overspeed factor (defined as the percentage increase of velocity compared to the undisturbed velocity at the same distance from the surface), As, on the top of the hill is shown in Fig. 6

I

1

.11

o

~

o•\

1~ x/h---- / \

\

\

~

\

i

1.2

i

13

- ~ °I"~

Overspeed f a c t o r

..Q.-~

~'"

i

1.4

~

~--~×

i

1.5

Fh-'- 2h

e s c a r p m e n t slope 2:1

O • M e a s u r e m e n t s [19]

o

o~"

\

-x/h=3 \

Fig. 5. Velocity predictions for the escarpment of 2:1 slope.

1.0

2_o!

3.0

y/h

4.0

5.0

9

16

I

L----SOre

0.5 Overspeed f a c t o r

h---- lore

\

Jackson and Hunt hill

Fig. 6. Predicted overspeed factor for the low JH hill.

y(m)

10

15

1.0

-q

¢o

318

,~i~o ...... L-:-I:I- ........................... /VTol X\

I

I

- 60 m

-- 30 rn

h:

. . . . . . . . . . . .

0

lOr'n

L:50m

,

3Ore--

T o p Of hill

Fig. 7. Predicted surface shear stress distribution for the low JH hill. 30 Jackson and Hunt hill h:20m

L=5Om

20

\ y(m)

\

10

\

o's Overspeed f a c t o r

Fig. 8. Predicted overspeed factor for the high JH hill.

10

I

.

.

.

6Ore

.

.

319

with the predictions obtained by Jackson and Hunt [1] and Deaves [6]. The near-surface overspeed factor predicted by the present m e t h o d is in excellent agreement with the JH theory, and in disagreement by about 30% with the predicted values of Deaves. A parameter more sensitive to the hill characteristics is the surface shear stress distribution. The predicted values over the hill are shown in Fig. 7. The maximum shear stress distribution over the hill predicted by the present method is in agreement with that predicted by the JH model and almost 50% in disagreement with the predictions of Deaves [6] and Taylor and Gent [4]. The method also predicts the upstream and downstream minima of shear stress which are theoretically expected. For the higher hill (h = 20 m), the predicted overspeed factor on the top of the hill (Fig. 8) is again in excellent agreement with the JH values and in good agreement with the predictions of Deaves, even though one should expect the differences in the latter to be more pronounced. For the shear stress distribution over the hill the predictions of the present method are shown in Fig. 9. Detailed predictions of the surface shear stress are not given in refs. 4, 6 and 7, however, their maximum values on the top of the hill are summarised in Table 2. 40 /

A --90

-60

-- 30

~

•

0

JH

[1]

Jackson and Hunt hill m

30

TOp of hill

Fig. 9. Predicted surface shear stress d i s t r i b u t i o n for t h e high J H hill. TABLE 2 M a x i m u m values o f the shear stress for t h e t w o hills Reference

h = 10 m

h -- 20 m

T a y l o r a n d G e n t [4 ] Deaves [6 ] Deaves [7 ] Jackson and Hunt [1] Present method

3.2 3.0 2.2 2.2 2.2

10.0 4.3 3.8 3.4 3.6

60

I 90 (m)

320 Again, the maximum predicted shear stress values are in very good agreement with those of the JH theory and the improved model of Deaves [7] whilst the corresponding values of the numerical approaches of Taylor and Gent [4] and Deaves [6] differ substantially. The present method is, in principle, free of the limitations of the numerical methods of previous workers. Also, numerical diffusion errors are minimal in the present method since, as the predictions show, the velocity vector is very well aligned in the majority of the flow field with the ~ (= constant) coordinate line. Conclusions

This paper has presented a numerical approach to solving the two-dimensional time-averaged Navier--Stokes equations over hills. The numerical solution is obtained in a " b o d y - f i t t e d " orthogonal curvilinear coordinate system, in which errors due to numerical diffusion are minimal, and boundary conditions and the shape of the boundary are properly taken into account. All terms in the governing differential equations are retained, hence, flow reversal in the flow domain can be predicted. An advanced, twoequation k--e turbulence model has been used, and any errors due to numerical convergence of the method and grid dependency have been minimised by grid dependency tests and convergence monitoring. Previous numerical methods do not have the generality of the present method, as was mentioned during the presentation of the results, and the discrepancies which have been noticed particularly for high hills can safely be attributed to their limitations. Comparison of the present predictions with those of the JH theory for slender hills under their limiting assumption, reveals the success of their theory. However, the present method for two-dimensional hills does not have any limitations and can be used for predicting flow fields over irregular two-dimensional hill shapes with a high degree of reliability. References 1 P.S. Jackson and J.C.R. Hunt, Turbulent wind flow over a low hill, Q.J.R. Meteorol. Soc., 101 (1975) 92~ ~55. 2 N.O. Jensen, Escarpment induced flow perturbations, a comparison of measurements and theory, J. Wind Eng. Ind. Aerodyn., 15 (1983) 243--251. 3 R.I. Sykes, An asymptotic theory of incompressible turbulent boundary layer flow over a small hump, J. Fluid Mech., 101 (1980) 647--670. 4 P.A. Taylor and P.R. Gent, A model of atmoa~heric boundary-layer flow above an isolated two-dimensional 'hill'; an e-Rmple of flow above 'gentle topography', Boundary Layer Meteorol., 7 (1974) 349--362. 5 P.A. Taylor, Numerical studies of neutrally stratified planetary boundary-layer flow above gentle topography. I. Two-dimensional c u e s , Boundary Layer M~taorol.,

12 (1977) 37--60. 6

D.M. Deaves, Wind over hills: a numerical approach, J. Ind, Aerodyn., 1 (1975/ 1976) 371--391.

321 7 D.M. Deaves, Computations of wind flow over two-dimensional hills and embankments, J. Wind Eng. Ind. Aerodyn., 6 (1980) 89--111. 8 J.L. Walmsley, Modelling of boundary-layer wind flow above complex terrain: a review, Environment Canada, Atmospheric Environment Service, Rep. AQRB83-005-L, 1983. 9 K. Antonopoulos, Predictions of flow and heat transfer in rod bundles, Ph.D. Thesis, Imperial College, 1979. 10 A.D. Gosman and C.W. Rapley, A prediction method for fully developed flow through non-circular passages, Proc. 1st Int. Conf. on Numerical Methods in Laminar and Turbulent Flows, University of Swansea, 1978. 11 S.V. Patankar and D.B. Spalding, A calculation procedure for heat, mass and moment u m transfer in three-dimensional parabolic flows, Int. J. Heat Mass Transfer, 15 (1972) 1787. 12 W.P. Jones and B.E. Launder, The prediction of laminarisation with a 2-equation model of turbulence, Int. J. Heat Mass Transfer, 15 (1972) 301. 13 B.E. Launder and D.B. Spalding, Mathematical Models of Turbulence, Pergamon Press, Oxford, 1972. 14 G. Bergeles, A.D. Gosman and B.E. Launder, Double row discrete-hole cooling: an experimental and numerical study, J. Eng. Power, 102 (1980) 498--503. 15 S.V. Patankar, D.K. Basu and S.A. Alpay, Prediction of the three-dimensional velocity field of a turbulent deflected jet, J. Fluids Eng., 99 (1977) 758. 16 G. Bergeles and N. Athanassiadis, Numerical study of the flow around a surface mounted prism, Proc. Symp. on Refined Modelling of Flows, Paris, 1982, Vol. 1, pp. 49--57. 17 L.P. Hackman, G.D. Raithby and A.B. Strong, Numerical predictions of flows over backward-facing steps, Int. J. Numerical Methods Fluids, 4 (1984) 711--724. 18 I.V. Nekrasov, L.H. Khurshudyan and W.H. Snyder, Soviet--American hill study project, U.S. Environmental Protection Agency, Research Triangle Park, EPA600/4-81-067, 1979. 19 A.J. Bowen and D. Lindley, A wind tunnel investigation of the wind speed and turbulence characteristics close to the ground over various escarpment shapes, Boundary-layer Meteorol., 12 (1977) 259--271.