- Email: [email protected]
Wind over hills: A numerical approach
Journal of Industrial Aerodynamics, 1 (1975[1976) 371--391 © Elsevier Scientific Publishing Company, Amsterdam -- Printed in The Netherlands
371
WIND O V E R H I L L S : A N U M E R I C A L A P P R O A C H
D.M. DEAVES Environmental Sciences Research Unit, Cranfield Institute of Technology, Cranfield, Bedford, MK43 OAL (Gt. Britain) (Received June 19, 1975; in revised form December 1, 1975)
Summary The equations of motion of atmospheric flow over a surface obstruction are formulated in terms of two flow variables, then written in finite difference form and solved numerically subject to specified closure assumptions and boundary conditions. Results for particular hill shapes are presented and are found to agree fairly well with experimental and theoretical published results. The method is valid for any shape of slender obstruction, whereas the methods of various other workers are less easy to apply to a given hill, requiring either a transformation of the lower boundary, or knowledge of an integral involving the slope of that boundary. The present approach is also shown theoretically to be insensitive to the closure and boundary conditions, provided these are realistic and consistent.
1. I n t r o d u c t i o n T h e e f f e c t s o f t o p o g r a p h y , and o f hills in particular, on profiles b o t h o f m e a n wind speeds and o f gusts, have always b e e n o f interest t o engineers conc e r n e d with the design and siting o f tall structures. T h e c u r r e n t r e c o m m e n d a tions o f the Building R e s e a r c h S t a t i o n [1] suggest t h e use o f a speed-up f a c t o r o f 1.1 f o r m o s t hill shapes; i.e. the m e a n w i n d speed at a given h e i g h t above the hill t o p is assumed t o be 1 0 ~ h i g h e r t h a n t h a t o f t h e i n c i d e n t profile at the same h e i g h t a b o v e level ground. T h e r e has b e e n a fair a m o u n t o f r e c e n t w o r k o n this p r o b l e m , m o s t o f which w o u l d suggest t h a t B.R.E.'s value is unlikely t o p r e d i c t a high e n o u g h speed-up f o r m o s t hill shapes. F r o s t , Maus and S i m p s o n [2] have solved n u m e r i c a l l y t h e t u r b u l e n t b o u n d a r y - l a y e r e q u a t i o n s f o r f l o w over relatively low ellipses, w h i c h were i n t e n d e d t o r e p r e s e n t such o b s t r u c t i o n s as a i r p o r t buildings, etc. T h e s e results, w h i c h also r e q u i r e d an estimate o f t h e pressure d i s t r i b u t i o n f r o m the inviscid flow solution, are t h e r e f o r e o f limited value t o the p r e s e n t p r o b l e m , since m o s t n a t u r a l hills w o u l d cause larger d i s t u r b a n c e s t o the flow. T h e m e t h o d o f J a c k s o n and H u n t [3] a m o u n t s t o analytical p e r t u r b a t i o n t h e o r y , taking w i n d profiles w h i c h have m e r e l y b e e n displaced over t h e h u m p as a first a p p r o x i m a t i o n t o t h e solution. Calculation o f speed-up e f f e c t s t h e n
372 involved taking an integral function of the shape of the hill, a procedure which evidently limits the applicability of the m e t h o d to analytically well-behaved hill shapes, even though the method may in principle be applied to any shape. Taylor and Gent [4] have produced a numerical method based upon similar assumptions to those of the present work. The method is limited, however, by their use of a conformal mapping to produce a straight-line lower b o u n d a r y to the grid system, with a resultant increase in the complexity of the system of equations. There is therefore a need for a simple method which may readily be applied to any shape of hill. The w o r k of Alexander and Coles [ 5], which provided a numerical m e t h o d for the solution of the Navier S t o k e s flow equations over an arbitrarily shaped hill, was subsequently amended by the author [6] and forms the basis of the present analysis. This numerical approach has been followed in preference to some of the more analytical attempts, since it has the advantage that the flow pattern over any shape of hill m a y readily be examined. 2. Derivation of equations of m o t i o n Although most hills will evidently induce 3D effects into the flow pattern, it would appear that such effects are likely to be quite small (see Deaves [6] ). It is also evident that consideration of the flow as taking place in t w o dimensions only (i.e. along-wind and vertically), will yield velocity profiles over the top of the hill which are upper limits to the full 3D solutions. This can be seen since, in the 3D case, flow around the sides of the hill will tend to reduce the speed-up experienced over the top. The equations of motion are therefore derived from the time-averaged 2D Navier-Stokes equations under the following additional assumptions : (i) In a relatively high speed flow, equivalent to the windiest conditions, the turbulent stresses, which vary as the square of velocity, will greatly exceed both viscous and Coriolis forces. These latter are therefore neglected, as are any density fluctuations. (ii) It is assumed that the turbulent stresses can be related to the mean-flow field by introducing a scalar eddy viscosity, e, where, for example, - u ' w ; = e
+
. The primes refer to the differences between instantaneous values
and the means of the appropriate quantities, which are defined after eqn. (2.3). (iii) The ratio of height (h) to "length" (L) of the hill shall be taken to be small enough that various simplifications can be made to the final forms of the equations. In practice, any such constraint imposed will apply to the maximum slope to be encountered. However, the streamlines will only follow that slope very close to the hill, having a slope whose absolute value is smaller away from the hill. Thus, it is to be expected that the procedure may still provide reasonable estimates of the flow field even if this condition is n o t satisfied (see results in section 7).
373 (iv) Separation, if it takes place, is n o t likely to be dealt with adequately. The m e t h o d is therefore only likely to provide accurate results for nonseparated flow, although, even for separated flow, the results may be valid sufficiently far upwind of the separation point. The basic equations of m o t i o n which may then be derived are those of continuity: ~u --+ 3x
3w
=0
3z
(2.1)
and those of m o m e n t u m (Navier--Stokes): +
-
3x -
-
-
-
~ - -
3z +
-
3x
-
p 3x -
3z
-
3x +
p 3z
(2.2)
3z
- -
(2.3)
3x
3z
where u = along wind (x) c o m p o n e n t of velocity w = vertical (z) c o m p o n e n t of velocity P = pressure p = density (constant) Equations (2.2) and (2.3) can then be combined by differentiating to eliminate P, and the following completely general form obtained:
3x
3z
3x
~z
3z +
_ 32
3x 32 [_u-~-Tw,] 3x 2
32
3. Relation of turbulent stresses to mean flow The equations which are used to relate these stresses to the mean flow velocities are those of Hinze [7]: -u'2
3
u'2 + v'2 + w'2
=2e3-x
-v ' 2 + 3
u '2 + v '2 + w '2
=0
(3.2)
-w '2+--
u ' 2 + v ' 2 + w '2
=2e-3z
(3.3)
r/p = -u'w'
= -w'u'
+
(3.4)
= e
~x
374 where v = cross wind (y) c o m p o n e n t of velocity. T = shear stress. The additional energy terms in eqns. (3.1) to (3.3) are included so that addition of these three equations does n o t imply that the turbulent energy (1/2p(u'2 + v'2 + w':)} is zero (as it does in Boussinesq's original formulation). Notice that, for a 2D flow, eqn. (3.2) implies that v '2 = l~(u'2 + w':), and ~U
that, for a steady flow with no x-dependence, ~x and
~}W
~z
are zero. Solution of
eqns. {3.1) to (3.3) then implies u': = v'2 = w '2
{3.5)
It is generally found that, for non-isotropic wall turbulence {where the turbulent flow has one solid boundary), eqn. {3.5) does n o t hold, but that u '2 > v '2 > w '2, although they are of the same order of magnitude for the greater part of the boundary layer. It will be noted from eqn. {2.4) that the only turbulence quantities required are ( w '2 - u '2 ) x z and derivatives of (r/p), and that any slight inconsistencies due to this formulation will only be present in the first of these terms. Since this term involves the x derivative of the difference of two turbulent stresses, it is assumed that eqns. {3.1) to {3.4) provide an adequate representation (see assumptions in section 5). Equation {3.4) is widely used in turbulent boundary~w au layer work, generally omitting the Ox term, since -az- p r e d o m i n a t e s for a boundary layer. If eqns. {3.1) to {3.4) are then inserted into eqn. {2.4), the following equation of motion is obtained: U--
~x = 4
ax
+W--
~z
- -
~z
ax
+
axaz
\ ~x
e
~z ~
~x 2
- -
~z
+
{3.6)
~x
In order to solve this, e must be specified in terms of u and w. Before looking at possible forms of e, the m e t h o d of solution, which is independent of t h a t form, is examined. Alexander and Coles [5] solved a similar equation by introducing the streamfunction $ and vorticity ~2 defined by: a~ ~z ' and ~2 -
a~ w = - -ax ~w
~u
ax
~z
(= -V2~)
However, if the upstream velocity is assumed to behave logarithmically,
{3.7) (3.8)
375
au -
-
1 au 1 ~---* ¢¢ as z ~ 0, implying t h a t ~2 -* - ~ as z -~ 0. At b e s t , ~ ~ -- which is
az
z
~z
z0
usually very large. Since ~2 is given finite values at discrete points in a grid, it is obvious t h a t this m e t h o d will introduce large errors at points near the lower au boundary. As z -* 0, although - - - * o% it is expected that e -~ 0 in such a way az au that r = p e - - tends to a finite value, rs, the surface shear stress. It was thereaz fore decided to formulate eqn. (3.6) in terms of ~, and the product: Q = -ea = e
ax
= eV
(3.9)
which has no singularities. Notice that Q will be a good estimate for aw au at points where - - is significant compared with - - . ax az Substitution of $ and Q into eqn. (3.6) then yields: V:$
= Q/e
V2 Q
= - -1
e
+2
rip except
(3.10) ( aU Q -ax
aQ) az
+ W ~
a2~ a2e ax 2 az 2
2--
--
a~ ~
axaz
Q ( U -be ~ ax
+W
a~z )
a2e a2~a~e -+ axaz
(3.11)
--az 2 a-xx2.J
These equations are independent of the form of e, which remains to be specified before solutions m a y be obtained. 4. Closure of the equations The original work of Alexander and Coles [5] employed an eddy viscosity e which was assumed to be constant with height. Now, near the ground surface, eqn. (3.4) becomes: au rip = e--
az
(4.1)
In the layer adjacent to the surface, the stress r is constant, as is p (see asau sumption 2(i)). Thus, if e is also constant, so is - - , which therefore implies az t h a t u = a z + b near z = 0. This assumption of constant e is evidently unrealistic, and results in a forced "linearisation" of the velocity profile over the hill. An excellent account of results of research on the nature of e is given by Launder and Spalding [8] where various classes of flow are examined. The closure of eqns.
376 (3.10) and (3.11) used in the present work is based upon the well known mixing-length hypothesis of Prandtl. This assumes that (4.2)
e = lm 2
where Im is a length, known as the mixing length. The local height is now defined by (4.3)
= z - hf(x)
where h f ( x ) is the height of the hill at horizontal position x, and it is assumed that
(4.4)
lm = k (~ + Zo)
where k = yon Karman's constant, z0 = roughness length of underlying terrain. Equation (4.4) is valid near the surface, and up to a height equal to about 1A of the boundary-layer thickness (around 300m, say). The implications of such a restriction to the applicability of the m e t h o d are discussed in the n e x t section. Equations (4.1), (4.2) and (4.4) may now be combined to yield: Ou k (~ + Zo) -~z = x / ( r / p ) = u .
in the constant-stress layer, where u , is the friction velocity. Integration of this immediately yields the familiar logarithmic velocity profile: u. u=
k
In I \
÷zo Zo
(4.5)
! /
which is incorporated as the upstream and downstream boundary conditions. The general form which is used for e is therefore obtained from (4.2) and (4.4) as: e = (k(~ +z0)~ 2
~-zUl
(4.6)
In order to effect a simplification before inserting eqn. (4.6) into eqns. (3.10) and (3.11), it is noted that, from eqn. (3.9), a u / ~ z ~ Q / e for most of the boundary-layer thickness. Substituting this into eqn. (4.6), the following is obtained: e = k ( ~ + Z o ) JQ ~
(4.7)
e has therefore been related quite simply to one of the flow variables, Q, which will be specified at each point of a finite grid.
377 Recently, Bradshaw [9] has suggested that the form of e over an obstacle will also depend upon additional strain-rate terms, so that: e= (k(5 + z0)) 2
°wl
~u +q~2 ~ x
(4.8)
where 4~2 is of order 10. Even if this is so, the form of eqn. (4.7) will still be valid, although there would be a more severe restriction on the ratio h/L than would otherwise have to be imposed (see comments in section 5). If eqn. (4.'/) is then substituted into eqns. (3.10) and (3.11), the following final forms of the flow equations are obtained :
Q V2qJ k(~" + Zo) IQ 11/2
(4.9)
aQ ~Q 2 Q ( - u h f ' + w) 2eV2Q = u - - + w - - +F ~x ~z z + z0
(4.10)
where
F
4e Lax 2 Oz 2
OxazSx~z
~z: ax:.J
The terms involving second derivatives of e are grouped together in the term F, and left in their unexpanded form, since it is found in practice that l~ is very small compared with the other terms involved in the calculation of Q. A number of calculations for identical hills, firstly including, then excluding, the F term, gave results which were identical to within the accuracy considered. Apart from these empirical results, there are analytical reasons for neglecting F, as each term is equivalent to a second derivative with respect to x, and these are small compared with the second z derivative present in the V2 term. 5. Assumptions of the analysis It was mentioned in section 2(iii) that there would be a restriction upon h/L in order that various terms become negligible. Before examining these restrictions, it is necessary to define L more explicitly. If h is the vertical, and L is the horizontal length scale, then it would seem sensible to define h/L as the m a x i m u m slope on the windward side of the hill. (If there is a very large, but localised, maximum, it may be better to estimate L from an averaged slope taken near the steepest part of the hill). The derivation of eqns. (3.10) and (3.11) have involved no assumption on the slenderness of the hill. Omitting the last term of eqn. (3.11) requires that: a2qJ < <
so that
< < 1.
(5.1)
378
The simplification made between eqns. (4.6) and (4.7) also requires the same condition, although, if: ~x 2 ~ ~u ~ Q ~z
e
-~
-~z2 ,ther~
[1 + [_
, so that:
e = k ( ~ , + e o ) IQ 11~
+-~
(5.2) (5.3)
(5.4)
Thus, it is required that (hJL) 2 < < 2, which is evidently satisfied if (h/L) 2 < < 1. If, however, eqn. (4.8) is valid, then: --- + ¢2 - - -~-~z ~x e
1+(¢2 +1)
(5.5)
so that the error in e is: 1/~(¢2 + 1)(h/L) 2 . Thus, (h/L) 2 << 2/(¢2 + 1) = 0.18 in the case where ¢2 ~- 10, and this is the strictest condition that may have to be applied. Since the relationship (4.8) is only expected to hold near the ground, three regions may be defined for the ratio h/L: (a)
(h/L) 2 << 0.18
(b)
(h/L) 2 << 1 b u t (h/L) ~ < / < 0.18
(c)
(h/L) 2 < / < 1
The m e t h o d should hold w i t h o u t reservation in the first region (h/L) < 1/6 (say), and may be expected to give reasonable estimates in the region (b), (h/L) < 1/3 (say), being less reliable above this limit. As has already been stated, in addition to this limitation on h/ L, eqn. (4.4) is only strictly correct up to around 300m. As the numerical grid generally extends to a height of around 10 times that of the hill (see comments in section 6), this implies that h <-- 30m in order that the equations e m p l o y e d remain valid throughout the computation region. However, it is shown in the Discussion (section 8) that the results are n o t likely to be sensitive to the exact form of closure, and therefore of lm, used, provided the upstream u, r and Im profiles are consistent and realistic. Since corrections to the logarithmic wind profile are likely to become significant at around 100m, it may be strictly necessary to amend all the boundary conditions when dealing with all b u t the most slender of hills. However, in view of the theoretical insensitivity of the results to such changes, it will be assumed that the present method is accurate for h <: 30m and reasonably so for h <-- 100m.
379
6. Boundary conditions Equation (4.9) is solved for ~b after setting Q to its previously c o m p u t e d (or initial) value. Equation (4.10) is solved for Q after linearising it by evaluating u, w and e from the previously computed values of $ and Q, together with the appropriate value of z. These two elliptic partial-differential equations can then be solved by writing them in finite difference form, and applying relaxation techniques to update the values of ~ and Q until some preset criterion for convergence is satisfied. The m e t h o d is basically the same as that used by Alexander and Coles [5], and is found to work well with these modified equations. The solutions obtained will evidently depend upon a sensible specification of the appropriate boundary conditions, which are given on the 4 boundaries shown in Fig. 1: G
(z : ~G)
T ..........
...... H (~=o)
Xo
x,u
~
........................... <
L
(x-xo)
~
2 I
(x:x~)
Fig. 1. Flow region, showing boundaries, uo(z) is the incident logarithmic velocity profile; L is a characteristic length scale of the hill, as defined in section 5.
X0 is where the incident logarithmic profile is specified. XI is where the velocity profile is assumed to have recovered its original logarithmic form. H is the ground and hill surface where the no-slip condition applies. This does n o t include all the line z -- 0. G is a level above which the flow is assumed to be undisturbed by the presence of the hill. The values of ~ and Q may therefore easily be specified on X0, XI, and G, being obtained by integrating eqn. (4.5). On X0, XI : =
u, L(Z r + z0)
In
\----~o z + z° l - z]
(6.1)
380 Q = u, 2
(6.2)
On G: =-~
(zo + z 0 ) I n
Z0
-ZG
Q = u, 2
(6.3) (6.4)
where za is the height of the level G above the ground. The lower boundary, H, does not coincide with a grid line, and also involves a rather less trivial specification of Q. Thus, on H, (6.5)
= 0
and Q is obtained from eqn. (4.9) as the limit as 2 -* 0: }-~0 [(z + z°)V~ lira (~2~/__1_1 ~ O'/~ =knm ~l = k ~--.o ~z21~,+Zo ]
(6.6)
~2¢ 02~ since ~z2 >> --Sx2 near ~ = 0. Now, consider:
lirn I
~-~0
~(z)
~ + z0 (~ + z0) l n - - -
Z0
] =A(say)
(6.7)
z
Since both numerator and denominator in this limit tend to zero as ~-~0, L'Hbpital's rule may be applied to obtair~ A=lim I
~0
~/az ]
In ( z +z°t
\
(6.8)
Z0 /
As the same still applies, further application of L'Hbpital's rule gives:
A = lim [ ~ / _ _ _ 1 _ 1
~
~0 L~z2/£, +Zoj
(6.9)
Combining eqns. (6.6), (6.7) and (6.8), the value of Q on H is :
o =
~-*o
(,~ +Zo) In
z + Zo
(6.10)
Zo
At any position over the hill, it is assumed that Q approaches its surface value smoothly, as it will upstream. Thus eqn. (6.10) is generally used with ~ as the
381 value appropriate to the nearest grid point, in the vertical direction, at which is non-zero. Setting a " l i d " at the level G represents a compromise, but is necessary for the numerical method. It is generally found necessary to take zG ~- a h + b L before the effects of the hill are sufficiently small, where a ~-- 10, b ~-- 0.1 Thus, for most hills, where h / L ~-- 0.1, the height of this lid m a y be taken to be approximately 10 times that of the hill. It is also found t h a t the results at the hill top are insensitive to the position of the downstream boundary XI, provided it is sufficiently far from the peak. Typical minimal distances are f o u n d to be around 5L. If a m a x i m u m of 50 grid points evenly spaced in the vertical direction is allowed in order to produce a workable computer program, only the lowest 5 will intersect the hill. Thus, for a hill of height 100m, the finest grid available would give 20m spacing. Since engineers concerned with wind flow over hills generally have the erection of some structure in mind, they are rarely interested in more than about the lowest 100m above the hill top. Evidently, such a grid with constant spacing would be much too coarse. To overcome this problem, it was decided to use a two-layer grid system, in which the lower region, up to a height of, say, twice that of the hill, used a refined grid, while the upper region, where gradients of ~ and Q are small, used a rather coarser spacing. However, such a scheme then poses various problems at the interface. For a grid with regular spacing, the Taylor series expansions of any function about a point in the grid may be combined in straightforward ways to provide first and second derivatives to reasonable accuracies. At the interface, where the grid spacing depends on the direction taken, more complicated expressions are required in order to achieve the same order of accuracy. Provided these expressions are computed carefully, the results obtained are reasonably continuous across the interface. The convergence criterion used is t h a t the values of both $ and Q at each point of the grid should differ from their previous values by less than a certain fraction of those values; i.e., (fn - fo)/fo < ~ where 5 ~ 0.003, and the suffices n and o refer to new and old values respectively. In these circumstances, it is found that convergence is obtained after between 30 and 150 iterations, which, for a 40 x 40 computation grid, takes about 5--25 minutes of ICL 1903T computer time. The wide range of expected times results from various values of h / L ; the higher this value, the greater the time required to obtain convergence of the solution. 7. Comparison o f results Computations have been made, using eqns. (4.9) and (4.10), for the flow over various shapes and sizes of hills. Comparisons with other results -- theoret~ ical, experimental and observed -- are shown in Figs. 2 to 9, notes on which are given below.
382 In all except Figs. 2 and 3, the "fractional speed-up ratio" As is used, since it gives a more sensitive indication of the agreement between the sets of results than do actual velocities. As is defined (as in Jackson and H u n t [3] ) by: As =
u(x,z) U(Xo, z - h f ( x ) )
1
(7.1)
where u (x0,z) is the upstream velocity profile (on boundary X0 ) h f ( x ) is the local hill height at position x. Note that eqn. (6.10) may be used in conjunction with eqn. (6.8) to express the surface value of As in terms of rs, the surface shear stress: Upstream surface shear stress = ro = p u , 2
(7.2)
Over the hill, rs--= Qs = lira ~
~-*0 L
a
(k ~ (z))
_~ 2
(~ +Zo) In z- - + z0 Zo
k
_ lim Ii -~-*0
n
.......
O~ 0Z ~ +Zo
2
(7.3)
Zo
Combining eqns. (7.2) and (7.3) then yields:
1;m - u (x,z) Ts/ro = z-*O -U~ In (z + z 0 t / -~
= l~m ~o
~ Zo
2
Ij
u (x,z)
] 2
(7.4)
_u [Xo. Z U - h f ( x ) l
Finally, eqns. (7.1) and (7.4) give: /1 +AS) 2 rs/r° =lim ~ 0 ~*
(7.5)
Equation (7.5) therefore provides a cross-check between the surface shear stress (PQs = Ts) and the value of As calculated from the $ values adjacent to the boundary. It should be noted from eqn. (2.4) that the pressure has been eliminated from the flow equations at an early stage. In fact, the only quantities used in the computation are $ and Q (see eqns. (3.10) and (3.11)), and, from these, local velocity values may be obtained. Also, from the Q values, and via eqns. (6.10) and (7.4), crude estimates may be made of the surface shear stress. In none of the examples presented below was separation predicted; w h e t h e r this
383
is a property of these flows, or of the particular numerical m e t h o d employed, is n o t certain• The sensitivity of As to changes in the mean wind speed has n o t been investigated, u , being given values appropriate to high windspeed, neutral atmospheric conditions ( u , = 1 m/s in most cases). For appreciably lower values of u,, it is quite likely that non-neutrality may a n y w a y induce departures from the incident logarithmic profile. In each case, values of h and L are given, L being chosen in accord with the definition of section 5. Notice that the values of L used by Jackson and Hunt [3] do n o t agree with this definition, but, as a result of the particular hill shapes they studied, are only 65% of those used in the present work (e.g. 50m in the examples of Figs. 7 to 9).
u(a) lu(100)
u (~) / u (100) 0"5
3"0
3"0
05
~6 14
I •o I • o/I 8 6
6
3
6 METRES
9
12
15
18
Fig. 2. W e s t c o t t e m b a n k m e n t v e l o c i t y profiles. - c o m p u t e d curve; . . . . . c o r r e c t e d for o b l i q u e i n c i d e n c e ; • • o b s e r v e d s p e e d s (4 runs).
c o m p u t e d curve
In Fig. 2, h = 5.4m, L = 13.5m and h / L = 0.4. Westcott is a disused railway e m b a n k m e n t in South Buckinghamshire and runs almost due N--S, situated in open c o u n t r y with Z0 ~ 5cm. The observations represent 4 different recordings, of average windspeed, 7.5 m/s at 10m, for each of which the incident wind direction was about 20 ° from normal. There were 2 anemometer masts, each with 3 sets of anemometers; one mast was at the f o o t of the slope, while the other was on top of the e m b a n k m e n t as shown. The theory over-estimates the velocities over the embankment, although, when corrected for oblique incidence (dashed curve), the agreement is rather better. In Fig. 3, h = 26m, L = 173m and h / L = 0.15. These measurements were taken on the BBC's transmission mast at Sutton Coldfield, which is set on a hillock about 30m above the surrounding terrain. Since there were no measurements upstream of the hill, the computed results correspond to an upstream roughness length of l m , which is reasonable considering the nature of the surrounding countryside. Again, there is a slight over-estimation of the hill-top velocities.
384
WfNDSPEED
(METRES/SEC) 10
WINDSPEED 20
(METRES/SEC) 10
I
20
i
I
100
o 40
20
~lOOm ~-
c o m p u t e d curve; • • observed speeds.
Fig. 3. S u t t o n Coldfield v e l o c i t y profiles. - -
In Fig. 4, h = 13m, L = 26m and h/L = 0.5. Comparison is shown with fullscale field tests [10] made 40 km N of Christchurch, N.Z. The escarpment is 13m high, with a slope of 26 ° ~-- tan -1 (1/~). The computed results actually correspond to an extended symmetrical embankment of total length 312m, so that the downstream boundary conditions do not quite correspond to thefield conditions. The theory again over-estimates the speed-up, except near the 10
DOS N A
posN
Q:
!
\
0.0
AS
0.5
00 i
~
B
AS
G5
•
A
B
65M ~ 6 5 M ~
Fig. 4. Sloping-escarpment fractional speed-up ratios. ( B o w e n and L i n d l e y [ 10 ]).
c o m p u t e d curve; . . . .
field data
385 ground at position A, where the relative coarseness of the numerical grid may affect the results. In Fig. 5, h = 7.5m, L = 15.0 and h / L = 0.5. The results presented here were obtained from wind-tunnel tests which were designed to model flow over an
I
~ i.-
P0S N A
2
00
l-
AS
0"5
2
0"0
POS N
B
05
3M
Fig. 5. M o t o r w a y - e m b a n k m e n t f r a c t i o n a l s p e e d - u p ratios, t u n n e l data.
c o m p u t e d curves; . . . .
wind
e m b a n k m e n t of height 7.5m, slopes of tan -1 (1/~) and total length 60m. These results are very similar to those shown in Fig. 4. In Fig. 6, h = 100m, L = 160m and h / L = 0.625. Computations are here compared with wind-tunnel results over a model hill [11], scaled up to a height of 100m. In spite of the relative "peakiness" of this hill, and the fact that the height is strictly outside the range of applicability of the method, a remarkably good agreement is obtained. Both Jackson and H u n t [3] and Taylor and Gent [4] give results for "bellshaped" hills, for which: f ( x ) = 1/(1 + ( x / L o ) : )
The remainder of the results refer to such hills. In Fig. 7, h = 20m, L = 77m and h / L = 0.26. Comparison is given with b o t h sets of results, for a low, b u t relatively peaked hill. Considering the sensitivity of As, the agreement of all three sets of results is quite good, especially near the surface where the greatest speed-up effects are felt. In Fig. 8, h = 10m, L = 77m and h / L = 0.13. Further comparisons are given of As for a hill of half the height of that in Fig. 7, for which Jackson and H u n t [3] therefore predict that As will be exactly halved. The present m e t h o d predicts values of As some 30% higher than theirs near the surface, although this corresponds to only a 10% increase in actual windspeed.
386
200
II
150
m I00
i.-
i
o m
50
i
I
.......
O5
i10
&5
Fig. 6. Sinusoidal-hill fractional speed-up ratios. Hill height = 1 0 0 m ; - - . . . . wind tunnel data (Plate and Lin [ 11 ]).
c o m p u t e d curve;
In Fig. 9, h = 10m, L = 77m and h/L = 0.13. Surface shear-stress profiles are shown, and compared with those given by Taylor and Gent [4]. The m a x i m u m value at the top of the hill, and those immediately upwind of this point, correspond very well, although the present method fails to predict the minima on the upwind and downwind slopes. Table 1 (L = 77m, L0 = 50m) compares the prediction of rs/ro at the hill top using the results of Jackson and Hunt [3], Taylor and Gent [4] and the present method. (Notice that, from (7.5), rs/To is even more sensitive than As). Since there is no clear correspondence between the different sets of results, evaluation of the relative merits of each method is difficult. TABLE 1 Comparison o f p r e d i c t i o n s o f normalised surface shear-stress m a x i m u m (rs/r o ) given b y different authors for bell-shaped hills w i t h L = 7 7 m (L 0 = 5Ore). h (metres)
5
10
20
T a y l o r and G e n t [ 4 ] J a c k s o n and H u n t [ 3 ] present m e t h o d
1.8 1.6 2.2
3.2 2.2 3.0
10.0 3,4 4.3
387
30
20
u~ tu
v
t-
w ca
o
O.S
00
10
AS
Fig. 7. Bell-shaped hill fractional speed-up ratios. computed curve; . . . . [4 ] computation; . . . . Jackson and Hunt [3 ] theory.
Taylor and Gent
8. Discussion It may be thought that the solutions obtained by the method given above will depend critically upon the boundary conditions and closure assumption used. It is acknowledged in particular that the constant boundary value of Q will be unrealistic for all but the m o s t slender of hills. The object of this discussion is to show the theoretical insensitivity of the results to these factors. The specification, eqn. (4.6), of the variation of e over the hill, is based on the mixing-length hypothesis. An alternative, which depends on this hypothesis only for the specification of the upstream and downstream viscosities, is to assume that, given these boundary values, e then remains constant along a streamline. Thus, in the undisturbed upstream flow where the velocity profile is assumed logarithmic, e = u , k ( z + Zo ). Since, for every value of z there corresponds a value of $, z may be written as a function of ~, and the eddy viscosity throughout the flow region may then be written: e = u.k(z(~)
+ zo)
(8.1)
388
1
1 I
\ \ I 0"5
0"0
\
I 10 AS
Fig. 8. Bell-shaped hill fractional speed-up ratios. Hunt [3] theory.
c o m p u t e d curve; .... Jackson and
3
Tsti °
E == u~
I -60
! -30 Df ST ANCE
FROM HILL
I 0 TOP
I 30
I 60
(METRES)
Fig. 9. Bell-shaped hill surface shear-stress profiles, Gent [4] computations.
computed curve; . . . .
Taylor and
389 where z ( $ ) is the solution of: =-~
(z(~) + z o ) l n
-z($)
z
(8.2)
Z0
The equations of motion then take the simplified forms : ey2~ = Q eV2Q = u
(8.3)
aQ
aQ r + w--+-ax az 2
(8.4)
where P is defined after eqn. (4.10). At any given stage in the computation, there is a value of $ at the grid point under consideration. Using this value, z ($) may be found from eqn. (8.2), and e, which is then used in eqns. (8.3) and (8.4), is given b y eqn. (8.1). Since eqn. (8.2) cannot be re-written to give an explicit formula for z ( $ ) , an iterative process must be used in its solution, thereby adding considerably to the computation time (by approximately doubling it). This formulation is seen to be superior to the previous assumption, eqn. (4.6), when conditions are examined at the top "lid" G of the finite grid. When e is constant on a streamline, because ~ is held constant on G (see eqn. (6.3)), so is e. However, with e related to ~, the local height, by eqn. (4.6), it is evident that e will vary as a function of x along this "lid". Since G was introduced as a level at which the incident flow remained undisturbed, such a variation of eddy viscosity there is an unwelcome feature of the model. Equations (8.1) to {8.4) have been used for the calculation of the flow over a few hill shapes, and were found to give almost identical results to those obtained from eqns. (4.7), (4.9) and (4.10). This close agreement would seem to result from the fact that, even using eqn. (4.7), the lines of constant eddy viscosity remain very close to the streamlines near the hill, although they tend to diverge higher up. As has been noted above, b o t h the forms of eddy viscosity used depend on the assumed upstream profiles of velocity and shear stress. Suppose n o w that these profiles are taken as any functions of height
(8.5)
u = Uo(Z), r/p = to(z) = Qo(Z)
Then, from eqn. (4.1) and (4.2),
Im
takes the form (8.6)
Im = toah/ ~ Uo / az
Inserting this into eqn. (4.2) and approximating a u / a z ~-Q/e yields
e = (to i Q I )~/~/auo /
az
Defining A as the fractional change in Q,
(8.7)
390
Q - Q0
(=
Q - to
)
eqn. (8.7) may be substituted into the equations of motion (3.10) to give:
V~-
~Uo Q ~u0 (1 + ~)1~ az (to iQ L)' ~ - ~z
(8.8)
The boundary conditions eqns. (6.2) and (6.4), which imply that the whole of the computation grid is within the constant-stress layer, were chosen for their ease of application. It is acknowledged that they are unlikely to be realistic for high hills, for which large values of h necessitate setting a high "lid" G on the computation grid. It is apparent that the results, in the form of velocities obtained from the values of ~, will be most sensitive to eqn. (4.9), which in turn depends upon I Q I '~. It has been shown in the derivation of eqn. (8.8) for generalised incident profiles that V2~b is independent of the upstream stress profile, varying only with z~, the fractional variation of Q from its incident form. The derivation of eqn. (8.8) would therefore seem to imply that the values of such quantities as As over the hill are n o t sensitive to the particular form of the upstream boundary conditions used. 9. Conclusions The present analysis provides a method of calculating the effects of arbitrarily-shaped surface obstacles -- in particular of hills -- on atmospheric wind flow. The results for which comparisons have been made show that the prediction of speed-up effects over hills are generally of the correct order of magnitude. Examination of the assumptions of the analysis would imply that results obtained for particularly large or steep hills should be treated with caution. It has been shown, theoretically, that the m e t h o d used is unlikely to be very sensitive to the assumed upstream conditions of stress and velocity profile. Since these quantities have been given variations appropriate to the lowest layers of the atmosphere, it is expected that the results will be most accurate in that layer, where the speed-up will be greatest. References 1 Building Research Station, The assessment of wind loads, Digest No. 119, H.M. Stationery Office, London, 1970. 2 W. Frost, J.R. Maus and W.R. Simpson, A boundary layer approach to the analysis of atmospheric motion over a surface obstruction, N.A.S.A. Rept. C.R. 2182, 1973. 3 P.S. Jackson and J.C.R. Hunt, Turbulent wind flow over a low hill, Q. J. R. Meteorol. Soc., 101 (1975) 929---955. 4 P.A. Taylor and P.R. Gent, A model of atmospheric boundary layer flow above an isolated two dimensional "hill": an example of flow above gentle topography, Boundary Layer Meteorol., 7 (1974) 349--362.
391 -
5 A.J. Alexander and C.F. Coles, A theoretical study of wind flow over hills, 3rd Int. Congr. on Wind Loads on Buildings and Structures, Tokyo, 1971. 6 D.M. Deaves, A theoretical study of three dimensional effects on wind flow over hills, M. Sc. Thesis, Loughborough University of Technology, England, 1973. 7 J.O. Hinze, Turbulence: An introduction to its Mechanism and Theory, McGraw-Hill, New York, 1959, p. 21. 8 B.E. Launder and D.B. Spalding, Mathematical Models of Turbulence, Academic Press, London, 1972. 9 P. Bradshaw, Variations on a theme of Prandtl, Imperial College Aero. Report 71--13, AGARD Meeting, 1971. 10 A.J. Bowen and D. Lindley, Measurements of the mean wind flow over various escarpment shapes, 5th Australasian Conference on Hydraulics and Fluid Mechanics, Christchurch, N.Z., 1974. 11 E.J. Plate and C.Y. Lin, The velocity field downstream from a two-dimensional model hill, Rept. No. CER65 EJP14, Colorado State University, 1965.
371
WIND O V E R H I L L S : A N U M E R I C A L A P P R O A C H
D.M. DEAVES Environmental Sciences Research Unit, Cranfield Institute of Technology, Cranfield, Bedford, MK43 OAL (Gt. Britain) (Received June 19, 1975; in revised form December 1, 1975)
Summary The equations of motion of atmospheric flow over a surface obstruction are formulated in terms of two flow variables, then written in finite difference form and solved numerically subject to specified closure assumptions and boundary conditions. Results for particular hill shapes are presented and are found to agree fairly well with experimental and theoretical published results. The method is valid for any shape of slender obstruction, whereas the methods of various other workers are less easy to apply to a given hill, requiring either a transformation of the lower boundary, or knowledge of an integral involving the slope of that boundary. The present approach is also shown theoretically to be insensitive to the closure and boundary conditions, provided these are realistic and consistent.
1. I n t r o d u c t i o n T h e e f f e c t s o f t o p o g r a p h y , and o f hills in particular, on profiles b o t h o f m e a n wind speeds and o f gusts, have always b e e n o f interest t o engineers conc e r n e d with the design and siting o f tall structures. T h e c u r r e n t r e c o m m e n d a tions o f the Building R e s e a r c h S t a t i o n [1] suggest t h e use o f a speed-up f a c t o r o f 1.1 f o r m o s t hill shapes; i.e. the m e a n w i n d speed at a given h e i g h t above the hill t o p is assumed t o be 1 0 ~ h i g h e r t h a n t h a t o f t h e i n c i d e n t profile at the same h e i g h t a b o v e level ground. T h e r e has b e e n a fair a m o u n t o f r e c e n t w o r k o n this p r o b l e m , m o s t o f which w o u l d suggest t h a t B.R.E.'s value is unlikely t o p r e d i c t a high e n o u g h speed-up f o r m o s t hill shapes. F r o s t , Maus and S i m p s o n [2] have solved n u m e r i c a l l y t h e t u r b u l e n t b o u n d a r y - l a y e r e q u a t i o n s f o r f l o w over relatively low ellipses, w h i c h were i n t e n d e d t o r e p r e s e n t such o b s t r u c t i o n s as a i r p o r t buildings, etc. T h e s e results, w h i c h also r e q u i r e d an estimate o f t h e pressure d i s t r i b u t i o n f r o m the inviscid flow solution, are t h e r e f o r e o f limited value t o the p r e s e n t p r o b l e m , since m o s t n a t u r a l hills w o u l d cause larger d i s t u r b a n c e s t o the flow. T h e m e t h o d o f J a c k s o n and H u n t [3] a m o u n t s t o analytical p e r t u r b a t i o n t h e o r y , taking w i n d profiles w h i c h have m e r e l y b e e n displaced over t h e h u m p as a first a p p r o x i m a t i o n t o t h e solution. Calculation o f speed-up e f f e c t s t h e n
372 involved taking an integral function of the shape of the hill, a procedure which evidently limits the applicability of the m e t h o d to analytically well-behaved hill shapes, even though the method may in principle be applied to any shape. Taylor and Gent [4] have produced a numerical method based upon similar assumptions to those of the present work. The method is limited, however, by their use of a conformal mapping to produce a straight-line lower b o u n d a r y to the grid system, with a resultant increase in the complexity of the system of equations. There is therefore a need for a simple method which may readily be applied to any shape of hill. The w o r k of Alexander and Coles [ 5], which provided a numerical m e t h o d for the solution of the Navier S t o k e s flow equations over an arbitrarily shaped hill, was subsequently amended by the author [6] and forms the basis of the present analysis. This numerical approach has been followed in preference to some of the more analytical attempts, since it has the advantage that the flow pattern over any shape of hill m a y readily be examined. 2. Derivation of equations of m o t i o n Although most hills will evidently induce 3D effects into the flow pattern, it would appear that such effects are likely to be quite small (see Deaves [6] ). It is also evident that consideration of the flow as taking place in t w o dimensions only (i.e. along-wind and vertically), will yield velocity profiles over the top of the hill which are upper limits to the full 3D solutions. This can be seen since, in the 3D case, flow around the sides of the hill will tend to reduce the speed-up experienced over the top. The equations of motion are therefore derived from the time-averaged 2D Navier-Stokes equations under the following additional assumptions : (i) In a relatively high speed flow, equivalent to the windiest conditions, the turbulent stresses, which vary as the square of velocity, will greatly exceed both viscous and Coriolis forces. These latter are therefore neglected, as are any density fluctuations. (ii) It is assumed that the turbulent stresses can be related to the mean-flow field by introducing a scalar eddy viscosity, e, where, for example, - u ' w ; = e
+
. The primes refer to the differences between instantaneous values
and the means of the appropriate quantities, which are defined after eqn. (2.3). (iii) The ratio of height (h) to "length" (L) of the hill shall be taken to be small enough that various simplifications can be made to the final forms of the equations. In practice, any such constraint imposed will apply to the maximum slope to be encountered. However, the streamlines will only follow that slope very close to the hill, having a slope whose absolute value is smaller away from the hill. Thus, it is to be expected that the procedure may still provide reasonable estimates of the flow field even if this condition is n o t satisfied (see results in section 7).
373 (iv) Separation, if it takes place, is n o t likely to be dealt with adequately. The m e t h o d is therefore only likely to provide accurate results for nonseparated flow, although, even for separated flow, the results may be valid sufficiently far upwind of the separation point. The basic equations of m o t i o n which may then be derived are those of continuity: ~u --+ 3x
3w
=0
3z
(2.1)
and those of m o m e n t u m (Navier--Stokes): +
-
3x -
-
-
-
~ - -
3z +
-
3x
-
p 3x -
3z
-
3x +
p 3z
(2.2)
3z
- -
(2.3)
3x
3z
where u = along wind (x) c o m p o n e n t of velocity w = vertical (z) c o m p o n e n t of velocity P = pressure p = density (constant) Equations (2.2) and (2.3) can then be combined by differentiating to eliminate P, and the following completely general form obtained:
3x
3z
3x
~z
3z +
_ 32
3x 32 [_u-~-Tw,] 3x 2
32
3. Relation of turbulent stresses to mean flow The equations which are used to relate these stresses to the mean flow velocities are those of Hinze [7]: -u'2
3
u'2 + v'2 + w'2
=2e3-x
-v ' 2 + 3
u '2 + v '2 + w '2
=0
(3.2)
-w '2+--
u ' 2 + v ' 2 + w '2
=2e-3z
(3.3)
r/p = -u'w'
= -w'u'
+
(3.4)
= e
~x
374 where v = cross wind (y) c o m p o n e n t of velocity. T = shear stress. The additional energy terms in eqns. (3.1) to (3.3) are included so that addition of these three equations does n o t imply that the turbulent energy (1/2p(u'2 + v'2 + w':)} is zero (as it does in Boussinesq's original formulation). Notice that, for a 2D flow, eqn. (3.2) implies that v '2 = l~(u'2 + w':), and ~U
that, for a steady flow with no x-dependence, ~x and
~}W
~z
are zero. Solution of
eqns. {3.1) to (3.3) then implies u': = v'2 = w '2
{3.5)
It is generally found that, for non-isotropic wall turbulence {where the turbulent flow has one solid boundary), eqn. {3.5) does n o t hold, but that u '2 > v '2 > w '2, although they are of the same order of magnitude for the greater part of the boundary layer. It will be noted from eqn. {2.4) that the only turbulence quantities required are ( w '2 - u '2 ) x z and derivatives of (r/p), and that any slight inconsistencies due to this formulation will only be present in the first of these terms. Since this term involves the x derivative of the difference of two turbulent stresses, it is assumed that eqns. {3.1) to {3.4) provide an adequate representation (see assumptions in section 5). Equation {3.4) is widely used in turbulent boundary~w au layer work, generally omitting the Ox term, since -az- p r e d o m i n a t e s for a boundary layer. If eqns. {3.1) to {3.4) are then inserted into eqn. {2.4), the following equation of motion is obtained: U--
~x = 4
ax
+W--
~z
- -
~z
ax
+
axaz
\ ~x
e
~z ~
~x 2
- -
~z
+
{3.6)
~x
In order to solve this, e must be specified in terms of u and w. Before looking at possible forms of e, the m e t h o d of solution, which is independent of t h a t form, is examined. Alexander and Coles [5] solved a similar equation by introducing the streamfunction $ and vorticity ~2 defined by: a~ ~z ' and ~2 -
a~ w = - -ax ~w
~u
ax
~z
(= -V2~)
However, if the upstream velocity is assumed to behave logarithmically,
{3.7) (3.8)
375
au -
-
1 au 1 ~---* ¢¢ as z ~ 0, implying t h a t ~2 -* - ~ as z -~ 0. At b e s t , ~ ~ -- which is
az
z
~z
z0
usually very large. Since ~2 is given finite values at discrete points in a grid, it is obvious t h a t this m e t h o d will introduce large errors at points near the lower au boundary. As z -* 0, although - - - * o% it is expected that e -~ 0 in such a way az au that r = p e - - tends to a finite value, rs, the surface shear stress. It was thereaz fore decided to formulate eqn. (3.6) in terms of ~, and the product: Q = -ea = e
ax
= eV
(3.9)
which has no singularities. Notice that Q will be a good estimate for aw au at points where - - is significant compared with - - . ax az Substitution of $ and Q into eqn. (3.6) then yields: V:$
= Q/e
V2 Q
= - -1
e
+2
rip except
(3.10) ( aU Q -ax
aQ) az
+ W ~
a2~ a2e ax 2 az 2
2--
--
a~ ~
axaz
Q ( U -be ~ ax
+W
a~z )
a2e a2~a~e -+ axaz
(3.11)
--az 2 a-xx2.J
These equations are independent of the form of e, which remains to be specified before solutions m a y be obtained. 4. Closure of the equations The original work of Alexander and Coles [5] employed an eddy viscosity e which was assumed to be constant with height. Now, near the ground surface, eqn. (3.4) becomes: au rip = e--
az
(4.1)
In the layer adjacent to the surface, the stress r is constant, as is p (see asau sumption 2(i)). Thus, if e is also constant, so is - - , which therefore implies az t h a t u = a z + b near z = 0. This assumption of constant e is evidently unrealistic, and results in a forced "linearisation" of the velocity profile over the hill. An excellent account of results of research on the nature of e is given by Launder and Spalding [8] where various classes of flow are examined. The closure of eqns.
376 (3.10) and (3.11) used in the present work is based upon the well known mixing-length hypothesis of Prandtl. This assumes that (4.2)
e = lm 2
where Im is a length, known as the mixing length. The local height is now defined by (4.3)
= z - hf(x)
where h f ( x ) is the height of the hill at horizontal position x, and it is assumed that
(4.4)
lm = k (~ + Zo)
where k = yon Karman's constant, z0 = roughness length of underlying terrain. Equation (4.4) is valid near the surface, and up to a height equal to about 1A of the boundary-layer thickness (around 300m, say). The implications of such a restriction to the applicability of the m e t h o d are discussed in the n e x t section. Equations (4.1), (4.2) and (4.4) may now be combined to yield: Ou k (~ + Zo) -~z = x / ( r / p ) = u .
in the constant-stress layer, where u , is the friction velocity. Integration of this immediately yields the familiar logarithmic velocity profile: u. u=
k
In I \
÷zo Zo
(4.5)
! /
which is incorporated as the upstream and downstream boundary conditions. The general form which is used for e is therefore obtained from (4.2) and (4.4) as: e = (k(~ +z0)~ 2
~-zUl
(4.6)
In order to effect a simplification before inserting eqn. (4.6) into eqns. (3.10) and (3.11), it is noted that, from eqn. (3.9), a u / ~ z ~ Q / e for most of the boundary-layer thickness. Substituting this into eqn. (4.6), the following is obtained: e = k ( ~ + Z o ) JQ ~
(4.7)
e has therefore been related quite simply to one of the flow variables, Q, which will be specified at each point of a finite grid.
377 Recently, Bradshaw [9] has suggested that the form of e over an obstacle will also depend upon additional strain-rate terms, so that: e= (k(5 + z0)) 2
°wl
~u +q~2 ~ x
(4.8)
where 4~2 is of order 10. Even if this is so, the form of eqn. (4.7) will still be valid, although there would be a more severe restriction on the ratio h/L than would otherwise have to be imposed (see comments in section 5). If eqn. (4.'/) is then substituted into eqns. (3.10) and (3.11), the following final forms of the flow equations are obtained :
Q V2qJ k(~" + Zo) IQ 11/2
(4.9)
aQ ~Q 2 Q ( - u h f ' + w) 2eV2Q = u - - + w - - +F ~x ~z z + z0
(4.10)
where
F
4e Lax 2 Oz 2
OxazSx~z
~z: ax:.J
The terms involving second derivatives of e are grouped together in the term F, and left in their unexpanded form, since it is found in practice that l~ is very small compared with the other terms involved in the calculation of Q. A number of calculations for identical hills, firstly including, then excluding, the F term, gave results which were identical to within the accuracy considered. Apart from these empirical results, there are analytical reasons for neglecting F, as each term is equivalent to a second derivative with respect to x, and these are small compared with the second z derivative present in the V2 term. 5. Assumptions of the analysis It was mentioned in section 2(iii) that there would be a restriction upon h/L in order that various terms become negligible. Before examining these restrictions, it is necessary to define L more explicitly. If h is the vertical, and L is the horizontal length scale, then it would seem sensible to define h/L as the m a x i m u m slope on the windward side of the hill. (If there is a very large, but localised, maximum, it may be better to estimate L from an averaged slope taken near the steepest part of the hill). The derivation of eqns. (3.10) and (3.11) have involved no assumption on the slenderness of the hill. Omitting the last term of eqn. (3.11) requires that: a2qJ < <
so that
< < 1.
(5.1)
378
The simplification made between eqns. (4.6) and (4.7) also requires the same condition, although, if: ~x 2 ~ ~u ~ Q ~z
e
-~
-~z2 ,ther~
[1 + [_
, so that:
e = k ( ~ , + e o ) IQ 11~
+-~
(5.2) (5.3)
(5.4)
Thus, it is required that (hJL) 2 < < 2, which is evidently satisfied if (h/L) 2 < < 1. If, however, eqn. (4.8) is valid, then: --- + ¢2 - - -~-~z ~x e
1+(¢2 +1)
(5.5)
so that the error in e is: 1/~(¢2 + 1)(h/L) 2 . Thus, (h/L) 2 << 2/(¢2 + 1) = 0.18 in the case where ¢2 ~- 10, and this is the strictest condition that may have to be applied. Since the relationship (4.8) is only expected to hold near the ground, three regions may be defined for the ratio h/L: (a)
(h/L) 2 << 0.18
(b)
(h/L) 2 << 1 b u t (h/L) ~ < / < 0.18
(c)
(h/L) 2 < / < 1
The m e t h o d should hold w i t h o u t reservation in the first region (h/L) < 1/6 (say), and may be expected to give reasonable estimates in the region (b), (h/L) < 1/3 (say), being less reliable above this limit. As has already been stated, in addition to this limitation on h/ L, eqn. (4.4) is only strictly correct up to around 300m. As the numerical grid generally extends to a height of around 10 times that of the hill (see comments in section 6), this implies that h <-- 30m in order that the equations e m p l o y e d remain valid throughout the computation region. However, it is shown in the Discussion (section 8) that the results are n o t likely to be sensitive to the exact form of closure, and therefore of lm, used, provided the upstream u, r and Im profiles are consistent and realistic. Since corrections to the logarithmic wind profile are likely to become significant at around 100m, it may be strictly necessary to amend all the boundary conditions when dealing with all b u t the most slender of hills. However, in view of the theoretical insensitivity of the results to such changes, it will be assumed that the present method is accurate for h <: 30m and reasonably so for h <-- 100m.
379
6. Boundary conditions Equation (4.9) is solved for ~b after setting Q to its previously c o m p u t e d (or initial) value. Equation (4.10) is solved for Q after linearising it by evaluating u, w and e from the previously computed values of $ and Q, together with the appropriate value of z. These two elliptic partial-differential equations can then be solved by writing them in finite difference form, and applying relaxation techniques to update the values of ~ and Q until some preset criterion for convergence is satisfied. The m e t h o d is basically the same as that used by Alexander and Coles [5], and is found to work well with these modified equations. The solutions obtained will evidently depend upon a sensible specification of the appropriate boundary conditions, which are given on the 4 boundaries shown in Fig. 1: G
(z : ~G)
T ..........
...... H (~=o)
Xo
x,u
~
........................... <
L
(x-xo)
~
2 I
(x:x~)
Fig. 1. Flow region, showing boundaries, uo(z) is the incident logarithmic velocity profile; L is a characteristic length scale of the hill, as defined in section 5.
X0 is where the incident logarithmic profile is specified. XI is where the velocity profile is assumed to have recovered its original logarithmic form. H is the ground and hill surface where the no-slip condition applies. This does n o t include all the line z -- 0. G is a level above which the flow is assumed to be undisturbed by the presence of the hill. The values of ~ and Q may therefore easily be specified on X0, XI, and G, being obtained by integrating eqn. (4.5). On X0, XI : =
u, L(Z r + z0)
In
\----~o z + z° l - z]
(6.1)
380 Q = u, 2
(6.2)
On G: =-~
(zo + z 0 ) I n
Z0
-ZG
Q = u, 2
(6.3) (6.4)
where za is the height of the level G above the ground. The lower boundary, H, does not coincide with a grid line, and also involves a rather less trivial specification of Q. Thus, on H, (6.5)
= 0
and Q is obtained from eqn. (4.9) as the limit as 2 -* 0: }-~0 [(z + z°)V~ lira (~2~/__1_1 ~ O'/~ =knm ~l = k ~--.o ~z21~,+Zo ]
(6.6)
~2¢ 02~ since ~z2 >> --Sx2 near ~ = 0. Now, consider:
lirn I
~-~0
~(z)
~ + z0 (~ + z0) l n - - -
Z0
] =A(say)
(6.7)
z
Since both numerator and denominator in this limit tend to zero as ~-~0, L'Hbpital's rule may be applied to obtair~ A=lim I
~0
~/az ]
In ( z +z°t
\
(6.8)
Z0 /
As the same still applies, further application of L'Hbpital's rule gives:
A = lim [ ~ / _ _ _ 1 _ 1
~
~0 L~z2/£, +Zoj
(6.9)
Combining eqns. (6.6), (6.7) and (6.8), the value of Q on H is :
o =
~-*o
(,~ +Zo) In
z + Zo
(6.10)
Zo
At any position over the hill, it is assumed that Q approaches its surface value smoothly, as it will upstream. Thus eqn. (6.10) is generally used with ~ as the
381 value appropriate to the nearest grid point, in the vertical direction, at which is non-zero. Setting a " l i d " at the level G represents a compromise, but is necessary for the numerical method. It is generally found necessary to take zG ~- a h + b L before the effects of the hill are sufficiently small, where a ~-- 10, b ~-- 0.1 Thus, for most hills, where h / L ~-- 0.1, the height of this lid m a y be taken to be approximately 10 times that of the hill. It is also found t h a t the results at the hill top are insensitive to the position of the downstream boundary XI, provided it is sufficiently far from the peak. Typical minimal distances are f o u n d to be around 5L. If a m a x i m u m of 50 grid points evenly spaced in the vertical direction is allowed in order to produce a workable computer program, only the lowest 5 will intersect the hill. Thus, for a hill of height 100m, the finest grid available would give 20m spacing. Since engineers concerned with wind flow over hills generally have the erection of some structure in mind, they are rarely interested in more than about the lowest 100m above the hill top. Evidently, such a grid with constant spacing would be much too coarse. To overcome this problem, it was decided to use a two-layer grid system, in which the lower region, up to a height of, say, twice that of the hill, used a refined grid, while the upper region, where gradients of ~ and Q are small, used a rather coarser spacing. However, such a scheme then poses various problems at the interface. For a grid with regular spacing, the Taylor series expansions of any function about a point in the grid may be combined in straightforward ways to provide first and second derivatives to reasonable accuracies. At the interface, where the grid spacing depends on the direction taken, more complicated expressions are required in order to achieve the same order of accuracy. Provided these expressions are computed carefully, the results obtained are reasonably continuous across the interface. The convergence criterion used is t h a t the values of both $ and Q at each point of the grid should differ from their previous values by less than a certain fraction of those values; i.e., (fn - fo)/fo < ~ where 5 ~ 0.003, and the suffices n and o refer to new and old values respectively. In these circumstances, it is found that convergence is obtained after between 30 and 150 iterations, which, for a 40 x 40 computation grid, takes about 5--25 minutes of ICL 1903T computer time. The wide range of expected times results from various values of h / L ; the higher this value, the greater the time required to obtain convergence of the solution. 7. Comparison o f results Computations have been made, using eqns. (4.9) and (4.10), for the flow over various shapes and sizes of hills. Comparisons with other results -- theoret~ ical, experimental and observed -- are shown in Figs. 2 to 9, notes on which are given below.
382 In all except Figs. 2 and 3, the "fractional speed-up ratio" As is used, since it gives a more sensitive indication of the agreement between the sets of results than do actual velocities. As is defined (as in Jackson and H u n t [3] ) by: As =
u(x,z) U(Xo, z - h f ( x ) )
1
(7.1)
where u (x0,z) is the upstream velocity profile (on boundary X0 ) h f ( x ) is the local hill height at position x. Note that eqn. (6.10) may be used in conjunction with eqn. (6.8) to express the surface value of As in terms of rs, the surface shear stress: Upstream surface shear stress = ro = p u , 2
(7.2)
Over the hill, rs--= Qs = lira ~
~-*0 L
a
(k ~ (z))
_~ 2
(~ +Zo) In z- - + z0 Zo
k
_ lim Ii -~-*0
n
.......
O~ 0Z ~ +Zo
2
(7.3)
Zo
Combining eqns. (7.2) and (7.3) then yields:
1;m - u (x,z) Ts/ro = z-*O -U~ In (z + z 0 t / -~
= l~m ~o
~ Zo
2
Ij
u (x,z)
] 2
(7.4)
_u [Xo. Z U - h f ( x ) l
Finally, eqns. (7.1) and (7.4) give: /1 +AS) 2 rs/r° =lim ~ 0 ~*
(7.5)
Equation (7.5) therefore provides a cross-check between the surface shear stress (PQs = Ts) and the value of As calculated from the $ values adjacent to the boundary. It should be noted from eqn. (2.4) that the pressure has been eliminated from the flow equations at an early stage. In fact, the only quantities used in the computation are $ and Q (see eqns. (3.10) and (3.11)), and, from these, local velocity values may be obtained. Also, from the Q values, and via eqns. (6.10) and (7.4), crude estimates may be made of the surface shear stress. In none of the examples presented below was separation predicted; w h e t h e r this
383
is a property of these flows, or of the particular numerical m e t h o d employed, is n o t certain• The sensitivity of As to changes in the mean wind speed has n o t been investigated, u , being given values appropriate to high windspeed, neutral atmospheric conditions ( u , = 1 m/s in most cases). For appreciably lower values of u,, it is quite likely that non-neutrality may a n y w a y induce departures from the incident logarithmic profile. In each case, values of h and L are given, L being chosen in accord with the definition of section 5. Notice that the values of L used by Jackson and Hunt [3] do n o t agree with this definition, but, as a result of the particular hill shapes they studied, are only 65% of those used in the present work (e.g. 50m in the examples of Figs. 7 to 9).
u(a) lu(100)
u (~) / u (100) 0"5
3"0
3"0
05
~6 14
I •o I • o/I 8 6
6
3
6 METRES
9
12
15
18
Fig. 2. W e s t c o t t e m b a n k m e n t v e l o c i t y profiles. - c o m p u t e d curve; . . . . . c o r r e c t e d for o b l i q u e i n c i d e n c e ; • • o b s e r v e d s p e e d s (4 runs).
c o m p u t e d curve
In Fig. 2, h = 5.4m, L = 13.5m and h / L = 0.4. Westcott is a disused railway e m b a n k m e n t in South Buckinghamshire and runs almost due N--S, situated in open c o u n t r y with Z0 ~ 5cm. The observations represent 4 different recordings, of average windspeed, 7.5 m/s at 10m, for each of which the incident wind direction was about 20 ° from normal. There were 2 anemometer masts, each with 3 sets of anemometers; one mast was at the f o o t of the slope, while the other was on top of the e m b a n k m e n t as shown. The theory over-estimates the velocities over the embankment, although, when corrected for oblique incidence (dashed curve), the agreement is rather better. In Fig. 3, h = 26m, L = 173m and h / L = 0.15. These measurements were taken on the BBC's transmission mast at Sutton Coldfield, which is set on a hillock about 30m above the surrounding terrain. Since there were no measurements upstream of the hill, the computed results correspond to an upstream roughness length of l m , which is reasonable considering the nature of the surrounding countryside. Again, there is a slight over-estimation of the hill-top velocities.
384
WfNDSPEED
(METRES/SEC) 10
WINDSPEED 20
(METRES/SEC) 10
I
20
i
I
100
o 40
20
~lOOm ~-
c o m p u t e d curve; • • observed speeds.
Fig. 3. S u t t o n Coldfield v e l o c i t y profiles. - -
In Fig. 4, h = 13m, L = 26m and h/L = 0.5. Comparison is shown with fullscale field tests [10] made 40 km N of Christchurch, N.Z. The escarpment is 13m high, with a slope of 26 ° ~-- tan -1 (1/~). The computed results actually correspond to an extended symmetrical embankment of total length 312m, so that the downstream boundary conditions do not quite correspond to thefield conditions. The theory again over-estimates the speed-up, except near the 10
DOS N A
posN
Q:
!
\
0.0
AS
0.5
00 i
~
B
AS
G5
•
A
B
65M ~ 6 5 M ~
Fig. 4. Sloping-escarpment fractional speed-up ratios. ( B o w e n and L i n d l e y [ 10 ]).
c o m p u t e d curve; . . . .
field data
385 ground at position A, where the relative coarseness of the numerical grid may affect the results. In Fig. 5, h = 7.5m, L = 15.0 and h / L = 0.5. The results presented here were obtained from wind-tunnel tests which were designed to model flow over an
I
~ i.-
P0S N A
2
00
l-
AS
0"5
2
0"0
POS N
B
05
3M
Fig. 5. M o t o r w a y - e m b a n k m e n t f r a c t i o n a l s p e e d - u p ratios, t u n n e l data.
c o m p u t e d curves; . . . .
wind
e m b a n k m e n t of height 7.5m, slopes of tan -1 (1/~) and total length 60m. These results are very similar to those shown in Fig. 4. In Fig. 6, h = 100m, L = 160m and h / L = 0.625. Computations are here compared with wind-tunnel results over a model hill [11], scaled up to a height of 100m. In spite of the relative "peakiness" of this hill, and the fact that the height is strictly outside the range of applicability of the method, a remarkably good agreement is obtained. Both Jackson and H u n t [3] and Taylor and Gent [4] give results for "bellshaped" hills, for which: f ( x ) = 1/(1 + ( x / L o ) : )
The remainder of the results refer to such hills. In Fig. 7, h = 20m, L = 77m and h / L = 0.26. Comparison is given with b o t h sets of results, for a low, b u t relatively peaked hill. Considering the sensitivity of As, the agreement of all three sets of results is quite good, especially near the surface where the greatest speed-up effects are felt. In Fig. 8, h = 10m, L = 77m and h / L = 0.13. Further comparisons are given of As for a hill of half the height of that in Fig. 7, for which Jackson and H u n t [3] therefore predict that As will be exactly halved. The present m e t h o d predicts values of As some 30% higher than theirs near the surface, although this corresponds to only a 10% increase in actual windspeed.
386
200
II
150
m I00
i.-
i
o m
50
i
I
.......
O5
i10
&5
Fig. 6. Sinusoidal-hill fractional speed-up ratios. Hill height = 1 0 0 m ; - - . . . . wind tunnel data (Plate and Lin [ 11 ]).
c o m p u t e d curve;
In Fig. 9, h = 10m, L = 77m and h/L = 0.13. Surface shear-stress profiles are shown, and compared with those given by Taylor and Gent [4]. The m a x i m u m value at the top of the hill, and those immediately upwind of this point, correspond very well, although the present method fails to predict the minima on the upwind and downwind slopes. Table 1 (L = 77m, L0 = 50m) compares the prediction of rs/ro at the hill top using the results of Jackson and Hunt [3], Taylor and Gent [4] and the present method. (Notice that, from (7.5), rs/To is even more sensitive than As). Since there is no clear correspondence between the different sets of results, evaluation of the relative merits of each method is difficult. TABLE 1 Comparison o f p r e d i c t i o n s o f normalised surface shear-stress m a x i m u m (rs/r o ) given b y different authors for bell-shaped hills w i t h L = 7 7 m (L 0 = 5Ore). h (metres)
5
10
20
T a y l o r and G e n t [ 4 ] J a c k s o n and H u n t [ 3 ] present m e t h o d
1.8 1.6 2.2
3.2 2.2 3.0
10.0 3,4 4.3
387
30
20
u~ tu
v
t-
w ca
o
O.S
00
10
AS
Fig. 7. Bell-shaped hill fractional speed-up ratios. computed curve; . . . . [4 ] computation; . . . . Jackson and Hunt [3 ] theory.
Taylor and Gent
8. Discussion It may be thought that the solutions obtained by the method given above will depend critically upon the boundary conditions and closure assumption used. It is acknowledged in particular that the constant boundary value of Q will be unrealistic for all but the m o s t slender of hills. The object of this discussion is to show the theoretical insensitivity of the results to these factors. The specification, eqn. (4.6), of the variation of e over the hill, is based on the mixing-length hypothesis. An alternative, which depends on this hypothesis only for the specification of the upstream and downstream viscosities, is to assume that, given these boundary values, e then remains constant along a streamline. Thus, in the undisturbed upstream flow where the velocity profile is assumed logarithmic, e = u , k ( z + Zo ). Since, for every value of z there corresponds a value of $, z may be written as a function of ~, and the eddy viscosity throughout the flow region may then be written: e = u.k(z(~)
+ zo)
(8.1)
388
1
1 I
\ \ I 0"5
0"0
\
I 10 AS
Fig. 8. Bell-shaped hill fractional speed-up ratios. Hunt [3] theory.
c o m p u t e d curve; .... Jackson and
3
Tsti °
E == u~
I -60
! -30 Df ST ANCE
FROM HILL
I 0 TOP
I 30
I 60
(METRES)
Fig. 9. Bell-shaped hill surface shear-stress profiles, Gent [4] computations.
computed curve; . . . .
Taylor and
389 where z ( $ ) is the solution of: =-~
(z(~) + z o ) l n
-z($)
z
(8.2)
Z0
The equations of motion then take the simplified forms : ey2~ = Q eV2Q = u
(8.3)
aQ
aQ r + w--+-ax az 2
(8.4)
where P is defined after eqn. (4.10). At any given stage in the computation, there is a value of $ at the grid point under consideration. Using this value, z ($) may be found from eqn. (8.2), and e, which is then used in eqns. (8.3) and (8.4), is given b y eqn. (8.1). Since eqn. (8.2) cannot be re-written to give an explicit formula for z ( $ ) , an iterative process must be used in its solution, thereby adding considerably to the computation time (by approximately doubling it). This formulation is seen to be superior to the previous assumption, eqn. (4.6), when conditions are examined at the top "lid" G of the finite grid. When e is constant on a streamline, because ~ is held constant on G (see eqn. (6.3)), so is e. However, with e related to ~, the local height, by eqn. (4.6), it is evident that e will vary as a function of x along this "lid". Since G was introduced as a level at which the incident flow remained undisturbed, such a variation of eddy viscosity there is an unwelcome feature of the model. Equations (8.1) to {8.4) have been used for the calculation of the flow over a few hill shapes, and were found to give almost identical results to those obtained from eqns. (4.7), (4.9) and (4.10). This close agreement would seem to result from the fact that, even using eqn. (4.7), the lines of constant eddy viscosity remain very close to the streamlines near the hill, although they tend to diverge higher up. As has been noted above, b o t h the forms of eddy viscosity used depend on the assumed upstream profiles of velocity and shear stress. Suppose n o w that these profiles are taken as any functions of height
(8.5)
u = Uo(Z), r/p = to(z) = Qo(Z)
Then, from eqn. (4.1) and (4.2),
Im
takes the form (8.6)
Im = toah/ ~ Uo / az
Inserting this into eqn. (4.2) and approximating a u / a z ~-Q/e yields
e = (to i Q I )~/~/auo /
az
Defining A as the fractional change in Q,
(8.7)
390
Q - Q0
(=
Q - to
)
eqn. (8.7) may be substituted into the equations of motion (3.10) to give:
V~-
~Uo Q ~u0 (1 + ~)1~ az (to iQ L)' ~ - ~z
(8.8)
The boundary conditions eqns. (6.2) and (6.4), which imply that the whole of the computation grid is within the constant-stress layer, were chosen for their ease of application. It is acknowledged that they are unlikely to be realistic for high hills, for which large values of h necessitate setting a high "lid" G on the computation grid. It is apparent that the results, in the form of velocities obtained from the values of ~, will be most sensitive to eqn. (4.9), which in turn depends upon I Q I '~. It has been shown in the derivation of eqn. (8.8) for generalised incident profiles that V2~b is independent of the upstream stress profile, varying only with z~, the fractional variation of Q from its incident form. The derivation of eqn. (8.8) would therefore seem to imply that the values of such quantities as As over the hill are n o t sensitive to the particular form of the upstream boundary conditions used. 9. Conclusions The present analysis provides a method of calculating the effects of arbitrarily-shaped surface obstacles -- in particular of hills -- on atmospheric wind flow. The results for which comparisons have been made show that the prediction of speed-up effects over hills are generally of the correct order of magnitude. Examination of the assumptions of the analysis would imply that results obtained for particularly large or steep hills should be treated with caution. It has been shown, theoretically, that the m e t h o d used is unlikely to be very sensitive to the assumed upstream conditions of stress and velocity profile. Since these quantities have been given variations appropriate to the lowest layers of the atmosphere, it is expected that the results will be most accurate in that layer, where the speed-up will be greatest. References 1 Building Research Station, The assessment of wind loads, Digest No. 119, H.M. Stationery Office, London, 1970. 2 W. Frost, J.R. Maus and W.R. Simpson, A boundary layer approach to the analysis of atmospheric motion over a surface obstruction, N.A.S.A. Rept. C.R. 2182, 1973. 3 P.S. Jackson and J.C.R. Hunt, Turbulent wind flow over a low hill, Q. J. R. Meteorol. Soc., 101 (1975) 929---955. 4 P.A. Taylor and P.R. Gent, A model of atmospheric boundary layer flow above an isolated two dimensional "hill": an example of flow above gentle topography, Boundary Layer Meteorol., 7 (1974) 349--362.
391 -
5 A.J. Alexander and C.F. Coles, A theoretical study of wind flow over hills, 3rd Int. Congr. on Wind Loads on Buildings and Structures, Tokyo, 1971. 6 D.M. Deaves, A theoretical study of three dimensional effects on wind flow over hills, M. Sc. Thesis, Loughborough University of Technology, England, 1973. 7 J.O. Hinze, Turbulence: An introduction to its Mechanism and Theory, McGraw-Hill, New York, 1959, p. 21. 8 B.E. Launder and D.B. Spalding, Mathematical Models of Turbulence, Academic Press, London, 1972. 9 P. Bradshaw, Variations on a theme of Prandtl, Imperial College Aero. Report 71--13, AGARD Meeting, 1971. 10 A.J. Bowen and D. Lindley, Measurements of the mean wind flow over various escarpment shapes, 5th Australasian Conference on Hydraulics and Fluid Mechanics, Christchurch, N.Z., 1974. 11 E.J. Plate and C.Y. Lin, The velocity field downstream from a two-dimensional model hill, Rept. No. CER65 EJP14, Colorado State University, 1965.