01 Cwe - Murakami 1990.pdf

517

Journal of Wind Engineering and Industrial Aerodynamics, 36 (1990) 517-538 Elsevier Science Publishers B.V., Amsterdam-Printed in The Netherlands

COMPUTATIONAL WIND ENGINEERING S. Murakami”

ABSTRACT The state of the art of computational wind engineering is reviewed. Firstly,- the numerical method for simulating turbulent airflow in wind engineering is described briefly. The diagnostic system for assessing the results of numerical simulation is described and error estimation and mesh resolution are discussed. Next the time-dependent flowfield is predicted by 3D Large Eddy Simulation and the results are illustrated using the techniques of animated graphics. Lastly an example, of practical application is described using the flowfield prediction around an actual building complex. INTRODUCTION A computational method for predicting the velocity field and pressure field around a building is investigated in this paper. In the field of engineering, new approaches based on numerical method have been developed rapidly in recent years. These approaches were made possible by the rapid advance of computer technology. The airflow around a building placed within an atmospheric surface boundary layer is fully turbulent and very complicated. It is composed of separations at the windward corners of the building, the circulations behind it, etc. There are many difficulties in clarifying these complex flowfields using usual wind tunnel techniques. Therefore, the establishment of a three-dimensional numerical method for simulating turbulent flow, namely computational wind engineering, is required for use in the field of wind engineering. RELATION BETWEEN COMPUTATIONAL AND EXPERIMENTAL WIND ENGINEERING The results of numerical simulation cannot be free from various types of numerical errors. There are many factors which produce numerical errors, for example, unfitness of the turbulence modelling, inaccuracy of the finite difference scheme, the macroscopic treatment of the viscous sublayer at the solid wall boundary, etc. Therefore it ‘Professor, Institute of Industrial Science, University 22-l 7-chome Roppongi Minato-ku, Tokyo, Japan

0167-6105/90/$03.50

0 1990-Elawier

Science Publishers B.V.

of Tokyo

518 is indispensable that the accuracy of numerical simulation be examined by comparing the numerical results with those from wind tunnel tests or field experiments. The numerical method cannot be developed and refined without the support of accurate experiments. On the other hand, the growing precisibn of prediction by numerical simulation also increases the incentive for new research into experimental methods. Therefore the two different methods of research should proceed in concert and in cooperation with each other. BRIEF REVIEW OF NUMERICAL METHOD Various types of Turbulence Models It is generally impossible to analyze directly the Navier-Stokes equations in their original form by the numerical method for a flowfield with large Reynolds number because that type of simulation, called direct simulation, calls for an unreasonably large number of grid points which can hardly be accomodated by existing super computers. Therefore, several types of averaged Navier-Stokes equations are usually used for numerical analysis. The first averaging type is the Reynolds averaging and the second is the filtering type. The various turbulence models are based on one of these types of averaging. Those turbulence models shown below are compared and evaluated from the viewpoint of wind engineering. (1) O-equation model This model is based on Reynolds averaged equations. As the modelling method is rather crude, it may be used very conveniently to give a rough prediction at an early stage of research. It is highly efficient in that it needs very little CPU time. (2) k - e two equation model This model, also based on the Reynolds averaged equations, is most widely applied to engineering problems. This model has a good reputation for reliability and seems to be most promising for present application to many problems in the field of wind engineering (Murakami et al.[ 1983], Yeung et al. [ 1985], Paterson et al.[ 1986], Mathews et al.[ 1987a, 1987b], Baetke et al.[ 1987], Murakami et al.[ 1988b, 1988e, 1988f]). The model equations are shown in Table 1. (3) Differential s t r e s s model, algebraic s t r e s s model These models are based on the Reynolds averaged transport equations for the Reynolds stress u,uj. Because these models are more sophisticated than those shown above, they require more sophisticated techniques for computation. As these models are expected to be able to analyze the turbulent flowfield very precisely, it is hoped to apply

519

these models to wind engineering problems in the future (Murakami et al.[ 1988], Kondo et al.[ 1988]. (4) Large Eddy Simulation Large Eddy Simulation(hereafter abbreviated as LES ) is based on the filtered Navier-Stokes equations. It may be applied to basic flows in the field of wind engineering (Deardorff [ 1970], Schumann [ 1975], Horiuti [ 1982], Moin et al.[ 1982], Murakami et al.[ 1984], Murakami et al.[ 1987]). This model has the advantage of generating a time-dependent flowfield. Its disadvantage is that it requires great deal of CPU time. The model equations are shown in Table 2. (5) Direct simulation This method, which uses t h e original Navier-Stokes equations, is applied to primitive, simple types of flowfields with a relatively low Reynolds number(below 104 at present), for example channel flows, etc (Kim et at.[ 1987]). Table 1

Model e q u a t i o n s f o r

k-e

model

(based on Reynolds averaged equations)

a-~'=o

~?

'

(l)

ax, ='

+

, p

I ,~(-~-~7.~)I

(2)

a h + a k~l= a__a__(2_ a k . 6b -~) * ~ --~- -c at ~x, a x , a,

(3)

a c a ~-~J a . . . . ( * . . ~ ) + C , ± ~ gb _C, ~' at F a x ~ ax, a, ~ k ~ k

(4)

k, ~ =C,,---

(5) au,

au~

a, =i.0. ¢,=l.& C,=0.09. C,=1.44. C,=1.92

Table 2

Model equations for LES (based on filtered Navier-Stokes equations)

B/~: =0 azl ~ ; + a f f , ~ _ ~ _ . ~ . _ _Pa 2 ..

at

(6)

azj--

=

!

ax,(p'3 ~ + a s ,

here v,cs=(Cz/P./-~) eu

(7)

am az~ 4 au~ ~z~

, C=O.1

(8)

520 In a strict sense, it is not easy to apply this method to problems of wind engineering where the Reynolds number is always very large. Sometimes quasi-direct simulation is conducted using the convection term of the third-order upwind scheme (Tamura et al.[ 1987]). Method of Discretization for Space Derivative Various methods of space discretization are used in computational fluid dynamics. (1)FDM (finite difference method) (2)CVM (control volume method) (3)FEM (finite element method) (4)spectral method In the case of FDM, it is easy to assess the accuracy of discretization and therefore it is also easy to develop a scheme with a higher-order accuracy. This method is most popular. CVM may be regarded as an integrated version of FDM. In CVM, the basic concept of the conservation of physical quantity is very clear. Because of this advantage this method has become the leading method in the field of computational fluid dynamics. When the flowfield is composed of the complex geometry, the generalized curvilinear coordinate is often used in FDM or CVM. FEM is good in matching the discretization to a complex-shaped boundary. In the case of FEM, large CPU memory storage is usually required in comparison with FDM. The spectral method is good in saving CPU time and also in accuracy. But it is difficult to apply this method to problems in wind engineering because it is suitable only for flowfields with very simple boundary conditions.

11 11

& ~ <-- <- T i

I [

,

.I. , ~ ,

~ ..~...'~

' ~ ~1 ~ " .~,

1 1

t 1

I!

i

"Vx~

0 1) F i r s t - o r d e r upwind, scheme with

coarse grid

grid

poi nt s=8X8x8. ( i nco,-rect: rosul

Fig. 1

difference scheme 12) C e n t r a d with iner grid 8rid

t )

,ol nts= correct

10×16x 18. result)

Indoor t u r b u l e n t flowfie]d w i t h d i f f e r e n t schemes for convection term (based on 3D k - ~ model)

521 Finite Difference Scheme for Convection Term and Numerical Viscosity In the space discretization of each term of the basic equation, the discretization of the nonlinear convection term is most important because the truncation terms of the discretized expression for the convection term work often as the numerical viscosity. They have great influence on the results of numerical simulation. The well-known forms of the finite difference scheme for convection term in a transport equation for scalar variable ¢ , are as follows(here velocity U is assumed to be positive and constant). (1)central difference scheme

-( °a-~-~-)' = - u I~ '"'-~2h ' - ' I "-U (¢9 ¢, 1+l,/,,,,h

=

I

.....

0 x " 5~ *]~-6¢, h (2)upwind difference scheme (first-order) -(08-~ ) , = - U

I# ' - # h 0 ¢

~-''1 ]

--U I(-~---~-)-~#" h + . . . . .

~-+.

....

)

) (9)

-(10)

(3)QUICK scheme (second-order upwind scheme)

-(~},=-U

IE~--(3~ , . , + 3 ¢ , - 7 ~ , - , + ¢ ,-=)l 0 ~ . ] .... =+1¢,., h =+. ) (11) - - U I(~-~-xJ+~-~ J~ . ~ .... The derivative of eventh-order, for example ¢ "', ¢ .... , works well the numerical viscosity because the e v e n t h - o r d e r derivatives possess the same function as the diffusion term which is composed of second derivative. The central difference scheme has no numerical viscosity of eventh-order. The f i r s t - o r d e r upwind scheme has the large numerical viscosity of ~ "U h / 2. As the QUICK scheme is very stable and has the relatively small numerical viscosity of -q~ .... U h 3/16, it is applied very widely to engineering problems with successful results. Fig. 1 compares numerical analysis of an indoor turbulent flowfield using different types of convection schemes based on a k - ¢ model. Fig. 1(2) (Murakami et al.[ 1986]) is exactly similar to the experimental results. The completely different flowfield reproduced in Fig. l(1)(Ishizu [ 1986]) shows the effect of numerical viscosity caused by the f i r s t - o r d e r upwind scheme with coarse mesh.

Macroscopic Wall Boundary Condition The value of the Reynolds number of the flowfield in wind engineering is generally high owing to the large length scale imposed, which makes it very difficult to set the mesh interval near the solid walls fine enough to resolve the viscous sublayer. For example, the first grid point from the wall is usually located at the area of 50-200 by value of wall unit y÷ (=u'y/v). Therefore the treatment of the tangential velocity component near the wall boundary is very

522

i m p o r t a n t . In t h i s t y p e of flow simulation, it has been proved t h a t t h e n o - s l i p b o u n d a r y condition does not p r o v i d e good r e s u l t s . Some t y p e of a r t i f i c i a l macroscopic b o u n d a r y conditions, such as log law t y p e o r power law t y p e etc, should t h e r e f o r e be i n t r o d u c e d to c o m p e n s a t e for t h e effect of t h e viscous sublayer. One of t h e most p o p u l a r b o u n d a r y conditions is t h e g e n e r a l i z e d log law, shown below.

U, (r

Ip

)(Co w

1/2 k ~)1~2=1g__ ~c

n (

E ~ h ,(C2"2k v

,)1~2}

(12)

This equation d e s c r i b e s t h e r e l a t i o n between U 1, k 1, ( r / p )w ( t a n g e n t i a l v e l o c i t y near wall, k i n e t i c e n e r g y near wall and wall s h e a r s t r e s s , r e s p e c t i v e l y ). DIAGNOSTIC SYSTEM FOR NUMERICAL SIMULATION The numerical simulation of t u r b u l e n t flow c o n s t i t u t e s a large s y s t e m composed of a g r e a t many s u b s y s t e m s which range from computational mathematics to t h e fluid dynamics of turbulence s t a t i s t i c s . It is v e r y difficult to j u d g e t h e d e g r e e of accuracy of t h e simulation r e s u l t because t h e many s u b s y s t e m s of t h e simulation a r e usually left as black boxes. Here, it is t h u s most i m p o r t a n t to develop a method for examining t h e accuracy of t h e simulation result. The examination method is h e r e called t h e " d i a g n o s t i c s y s t e m ". In t h i s paper, we g i v e two examples of such diagnoses. F i r s t is an e s t i m a t i o n of e r r o r s caused by f i n i t e difference. The second is t h e e f f e c t of t h e fineness of mesh r e s o l u t i o n on t h e simulation r e s u l t s . Numerical simulations a r e h e r e c o n d u c t e d by a 3D k - E model (Murakami et al.[ 1988e]). A s t a g g e r e d g r i d is used (MAC method, Harlow et al. [ 1965]). A simultaneous i t e r a t i o n method for v e l o c i t y and p r e s s u r e is a d o p t e d (ABMAC method, Viecelli [ 1971]). The QUICK scheme is used for c o n v e c t i v e t e r m of t h e scalar t r a n s p o r t equations. A power law t y p e of macroscopic boundary c o n d i t i o n is used for t h e wall boundary. All f i n i t e d i f f e r e n c e s for space and time a r e t r e a t e d w i t h s e c o n d - o r d e r a c c u r a c y c o n c e r n i n g t h e e r r o r estimation. E s t i m a t i o n of T r u n c a t i o n E r r o r ~ h and Solution E r r o r ~" The d i s t r i b u t i o n of two t y p e s of e r r o r s , namely, s o l u t i o n e r r o r s and t r u n c a t i o n e r r o r s , are e s t i m a t e d by Richardson e x t r a p o l a t i o n . Solution e r r o r is t h e residual between t h e exact s o l u t i o n and t h e s o l u t i o n g i v e n by t h e f i n i t e d i f f e r e n c e method. In t h e method used here, we assume t h a t t h e s o l u t i o n g i v e n by t h e f i n i t e d i f f e r e n c e method can be e x p r e s s e d as a Taylor series. If t h e f i n i t e d i f f e r e n c e scheme has a second o r d e r accuracy, t h e s o l u t i o n e r r o r e , of mesh s i z e h can be easily e s t i m a t e d w i t h a small amount of algebra as follows (Caruso et al.[ 1986], Murakami et al.[ 1988a]), ~ = ( U i - U 2J )/3, (13) w h e r e U , is t h e s o l u t i o n of mesh s i z e h , and U 2~ is t h e s o l u t i o n of mesh s i z e 2 k and ~", is t h e e s t i m a t e d s o l u t i o n e r r o r . T r u n c a t i o n e r r o r s a r i s e in f i n i t e d i f f e r e n c e approximation. The

523

t r u n c a t i o n e r r o r r , is defined as the residual obtained by s u b s t i t u t i n g t h e exact s o l u t i o n of t h e d i f f e r e n t i a l e q u a t i o n U e× into t h e f i n i t e difference equation, i.e., r , = L , [ U E x ] - L [ U =x] (14) where U e x is t h e exact solution, L is the spatial d i f f e r e n t i a l operator, and L ~ is t h e spatial difference o p e r a t o r of mesh size h . Here t h e exact solution U e× is estimated by U

U=x

~x=U

(15)

~+ e ~ .

is used in place of exact solution U Ex in eq.(14) to estimate

A ........................... ! o ~ ........................... • . ~ upper bound a ry . . . . . . . . . . . . . . ,,:,,~,.~ .. ~ ! !.! ! ! !.! ! ! ,., , . . . .H ~ _ , , z v , hl~'~[d.f~ ................ H + H + H + V , ; .... downstream

m e ." . . . . . . . . . . . .+. .~., ., ~. . . . ..................... : ............ . . . . . . . . . . . ~.~.~) I........-~ • ' [,):~:Z.~, ', I ', ', I ', ', ',

' ' ' ', ', ', ', ', . . . . .

~,,~,,

: ', ', ', ', ', ', ', ', ', ', ', ', I ',

, , , I ', ', ', ', ', ', ',

,~!iii!!!!!iiiiii!!~,

~!!i!!!!!!iii! 7

"Eround

cubic / model

Fig. 2

Mesh i (vertical plane)

z~

0

(l ) t ~ n c ~ t i o n

error

~-

:.:.:.:.:.:.." ".::::~,

::::iiiiii;. . . . . !i!i!iiiiiiiii!i:i:i:i:!!i!!i!~iiiii i i!iii ! !!! iiK":.:-.":~o. o [ ~'"'" "'"""""'; ..... ~ xN~

(2)

solutioa e ~ c o r

Non-dimensional herefore, O7~. o f t h e

Fig.3

~,~,~ '::;~

:

'")-,~ ~ ..

~~ • u,c

' '":':*:.:.:.:.:.:.:v.-. .-~ ...... .....:

T i

velocity st the maximum v a l u e velocity value.

a

helsht T i

of

of does

H

I not

is l.O.X~ e×ceed J

D i s t r i b u t i o n of e s t i m a t e d e r r o r s f o r absolute v e l o c i t y (J U=+V=+W =)

/

524

t h e t r u n c a t i o n error: r ,=L,[U~xl-t[Ut~x] (16) If a very a c c u r a t e solution is obtained with mesh size h , the s o l u t i o n e r r o r e% is estimated by eq.(13) to be v e r y small at all mesh points. T r u n c a t i o n e r r o r s are regarded as the source of s o l u t i o n e r r o r s and are supposed to be c o n c e n t r a t e d in t h e region where large g r a d i e n t s of variables exist. If we can decrease t r u n c a t i o n e r r o r in t h i s region, s o l u t i o n e r r o r will become smaller in the whole computational domain. Fig. 3 i l l u s t r a t e s d i s t r i b u t i o n s of e r r o r s with Mesh 1 (h ~=Hb/6) p r e s e n t e d in Fig. 2(Murakami et al.[ 1988c]). Fig. 3(1) shows t h e d i s t r i b u t i o n of t r u n c a t i o n e r r o r r ,. r ~ c e n t e r s in the region in f r o n t of the model. The d i s t r i b u t i o n of solution e r r o r b-, is given in Fig. 3(2). High values of t r u n c a t i o n e r r o r r i s i n g a r o u n d t h e windward c o r n e r are convected and diffused; as a result, r a t h e r high values for s o l u t i o n e r r o r are d i s t r i b u t e d in t h e whole flowfield around t h e model. However, the maximum value of ~ ' , does not exceed 10% of t h e velocity value. It is clarified from t h e s e r e s u l t s t h a t t h e mesh

~! :~, H]I

!

.H%tt,flHBtttlHltH,I< ~t .~[

II[I;

r l .... '~,',~" "lllllZ~_~b<m"d'" ~(z)~Z .....

I i I . . .I. --l>oundlry

....~

N! ;! II IIIIIIIIIIIIIIIIIIIlUlIIIIIIII I I I I IIIIII IIII I I I I 11- I....... ,,,,,IttI'HH{ ........ [I I*:':~'I I I I I I I IIIIIIIIIIIIIInall I I I I I I

.:::~:. :~ :! :~:~:~: :~~!!!~!-!-?.:.!!':!'.'ii i i i ~[ ~ ::::::::::::::::::: : : : : : :

!

:::::

:

.":'.:L: ~ ~ ! ! ! ! ~ H ~ : ' " " ' " ~ ! ~ ! !

i i i i

/

~bic

.....

!

!

I

: :: :: ::

~'Svou.nd

. . . .

model

(1)

Mesh 2 Is 5(x)x37(y)x2 I(z)=34,965 h ~ = H b /6, h ~ is the m i n i m u m

ize adjacent to the model.

~ ,,,,,,,IHt~,.~..~u~,,,, ~..~| II ':.':"2~:"fllllllllllllll[I ~,h?~,?~.-41111111111111]l ~.~ I IIIIIIIIIgillllllllllllllll ~".~ I IIIIIIIIlillllllllllllllll fX~'I'dLr4~l I I I Illllllllllll I I I I .....

III

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

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

(2)

Mesh

I I

i~

\

mesh~ / ¢

I I I I I I I I I I J I I I I I I II II II II II II II I

~ ~ , , ~ ~ , [ I

:.i~!:---=--= - ; - ~ = m - = ~ - -

. . . . . .

3

50(x)x49(y)x28(z)=68,600, h ~=24/Hb The mesh distribution is concentrated J near the windward corner. / Fig.4

t-~/ ~--I~ J_LI] LIII I-I~} I I .

Mesh dividings (vertical plane)

525 resolution around the windward corner of the model is one of the most important factors in predicting the flowfield with small error. Effect of Sensitivity of Mesh Resolution on Simulation Results Numerical simulation is conducted using two different types of mesh dividings and compared with wind tunnel experiments (Murakami et al. [ 1988e]). The mesh dividings used here are shown in Fig. 4. The mesh interval h i adjacent to the model is Hb/6 with Mesh 2. With Mesh 3, the mesh interval h i around the windward corner is Hb/24. The velocity vectors around the model are compared in Fig. 5(I),(2). The entire flow pattern in both vertical section and horizontal plane is well reproduced in the result using Mesh 2. However, the reverse flow on the roof and near t h e side wall which exists in the experimental data does not appear in the results of Mesh 2. The pressure coefficients Cp are compared in Fig. 6. In the case of Mesh 2 ( - O - ) , the positive pressures on the windward wall are overestimated and the negative pressures on the top and on the side walls are underestimated in comparison with the experimental data. The results using Mesh 3 are also presented in Figs. 5 and 6. With Mesh 3, the mesh distribution is concentrated near the windward corner. The reverse flows on the roof and near the side walls are clearly reproduced. Surface pressure distribution in Fig. 6( - ~ 7 - ) also corresponds very well to the experiments. But some differences in the flowfield of the wake remain to be improved. PREDICTION OF TIME-DEPENDENT FLOWHELD BY LES AND VISUAL ANIMATION In this section, the time-dependent flowfield predicted by 3D LES is presented first. Some results of visual animation are illustrated in the latter part. Building Model and Outline of Computational Method The arrangement of the building model is illustrated in Fig. 7. The cubic shaped blocks are set regularly. Some types of "street canyons" occur between the blocks (Murakami et al.[ 1987b]). The staggered grid system is used here. All spatial derivatives are approximated by central difference (second-order). The Adams-Bashforth scheme (second-order) is used for time-marching. Numerical integrations are conducted following the ABMAC method, a simultaneous iteration method for pressure and velocities. The mesh dividing is shown in Fig. 8. The number of grid points amounts to 39(x)x36(y)x21(z)=43,524. The mesh interval h I adjacent to the wall surface is Hb/20. One unit of block area ( one building plus half of the width of the surrounding streets) is selected for the simulations, and cyclic boundary conditions are adopted here.

526

j - - - - - - ~

-

(~) e x p e r i m e n t (MurakamJ

~

.

~

.

.

. et

experiment

al.. Re=7xlO')

-E

> ~" > > , >>> > ¢ -->----

-.~\

~-I~U~I'

4 .

.

.

.

.

' .......

:I1[

. . . . . . .

Mesh 2

(2) Mesh 2

~

(3) Mesh 3 (1) vertical

(3) Mesfl 3 (2) h o r i z o n t a l

plane at c e n t e r of the model Fig.5

plane

a t Z=H b / 2

Comparison of velocity vectors

-L.O 7*'

~ ~

CL

oxperi mont. Mesh 2

W.////i!//////~ ~ ~Mesh

3

4- LO +0.5

$

0.0

--0.5 - I.O

- 0 . ! ~

+1.0 ~0.5

(1)

0.0~ ~ / ~ .

vertical

plane Fig.6

t

--(].';-L.O

--L,O

(2)

horizontal a t l =Hb / 2

Comparison of surface pressure coefficients (based on k - e )

plane

527

bui Iding mode!

l~reei:

rind

Y

L_,

Fig.7

Arrangement of building models (shape of building is cubic)

4Hb wind

7Hb

1.5 4Hb

J /

,x (1)

L__. x ~ Hb "~ Hb :"

vertical Fig.8

plane

2Hb

(2) h o r i z o n t a l

plane

Finite difference mesh system (grid points=43,524)

The b o u n d a r y c o n d i t i o n s a r e as follows: (1) Downstream and s p a n w i s e d i r e c t i o n s : cyclic b o u n d a r y c o n d i t i o n s . Gross d o w n s t r e a m p r e s s u r e g r a d i e n t is imposed. Here, AP*=0.011. AP* means t h e p r e s s u r e d i f f e r e n c e n o n - d i m e n s i o n a l i z e d by pUb 2. Ub means time a v e r a g e d v e l o c i t y at inflow b o u n d a r y of c o m p u t a t i o n a l domain at h e i g h t of Hb. (2) Upper s u r f a c e of c o m p u t a t i o n a l domain: f r e e - s l i p b o u n d a r y condition. (3) Wall boundary: t h e p r o f i l e s of t h e t a n g e n t i a l v e l o c i t y c o m p o n e n t s a r e assumed to obey a log law e x p r e s s e d as eq.(12) near t h e wall. Normal v e l o c i t y component at t h e b o u n d a r y is assumed to be zero.

528

/

i. O0

displayed

area

V1 O.

O0

- I . O0 O. O0

25. O0

"~0.O0

~., 75. O0

experiment (Murakami e t a l . ,

(1)

:i:i:i:::~J ~!~i~i~!~ ~::::'2"

_

~

.... . . . . . .

~

--"

R e - 2 . S x l O ~)

l. O0 i!'

•

"

.

• fl

m

(i) experiment ~N!iii/~

-l. O0 ~

'

O. O0

50. O0

numerical

Fig.9

,

. . . . . .

e~ a [ . )

(Murakami

.__>

....

. . . . . . . .

I

25. O0

(2)

:

. . . . . . . .

c;" 75. O0

simula~:ion

Velocity fluctuation at point 1 (cf.Fig.8)

l.O

n~(n)

(u',)

i" " , II If!" II t _. c..+ 12) simulation

F i g . 11

"" .. :

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

Comparison

of

velocity

vectors

0.|

~-.,~ . . . . . . .

mean

~

O. Ol

I

O. 00! 0.01

•

.,++,,*l

Fig. I0

,

O. I

. . . . . . .

,~))~s

1.0

fl

+ " • ,,++~i,

IO.O ..~/Ub

Velocity spectrum of u component at point 1

r ='

Fig.12

II

r

Mean velocity vectors with finer mesh dividing (14XlO4 grid points)

Time-Dependent Velocity Fluctuation and Velocity Spectrum Fig. 9 illustrates the time-history of the velocity fluctuations at point 1 presented in Fig. 8. In this figure, t* denotes the non-dimensional time scale tUb/Hb. Fig. 10 shows the velocity spectrum of u component at the same point. The predicted spectrum corresponds very well with that of the experiment. Mean Velocity Field Time-averaged velocity vectors are illustrated in Fig. II. In general, the entire flow pattern from the numerical simulation agrees well with that from the experimental results. But slight differences do exist in the wake behind the model. The wake region is a little larger than that in the experimental result.

529

.~:. .". . . . .~ ~. :

.~.

d ~:" 0, ~ ~

i:....... ":~i~i:=i~

" -.~

'~ ~@:~i~ ~ ~i' ~::..:..~ '~N~"

~.....

' i:i:i:~:!:::::::::::::::::::::::i:i:i:i!:~:

i) experiment

i:!:i::::::.. •

(2)

......

s~mulat[on l~,XlO"grid points)

Fig. 14

Pressure coefficients horizontal plane) at z=Hb/2

-2.0~ ~ I H ~

....:::::::: '."!..:::. ~ . ~ Ii!ii i.ii~i iiiii~:i~iii.il li~ii!:i:i!iii:iiiliiiiil:i:ii!i!ii:i~~:l!:i:i:~:~:~:i:i:i:i:i:i:i:i:i:i:i:i!ii

-,o- 3Hb

(3) simulation (l~,xlO ~ grid point;s) Fig.13

Turbulent kinetic energy ]=

1.0

-1.I -2.0 Fig.15

~

~ ~

.~

-120

Preasure coefficients with various buildingdistances (simulation~

l}lb 3tlb etc means the distance between) the building models, cf. Fig.7. horizontal plane at z=llb/2

At the windward corner on the roof, a r e v e r s e flow should appear but is not reproduced in Fig. II(2). Fig. 12 shows mean velocity v e c t o r s with finer meshes (60(x)x51(y)x46(z)=140,760), in which the mesh interval adjacent to the model surface is Hb/40. The r e v e r s e flow on the roof surface is here well reproduced. Distribution of Turbulent Kinetic Energy The d i s t r i b u t i o n s of turbulent kinetic energy /z are compared in Fig. 13. Fig. 13(2) shows the result with about 4xlO 4 grid points and Fig. 13(3) shows the result with 14xlO 4 grid points. The result with 4xlO 4 grid points shows some difference when compared to the experimental result, especially for the k d i s t r i b u t i o n in the v i c i n i t y of the windward corner. The correspondence between the result obtained with finer meshes (14xlO 4 grid points) and t h a t from the experiment is confirmed to be v e r y good.

530

Surface Pressure Fig. 14 shows distributions of pressure coefficient Cp. In this study, the static pressure at point z=5Hb at the inflow boundary of the computational domain is used as the reference pressure Po for both the experiment and the numerical simulation. Pressure distribution is rather well reproduced by the numerical simulation. In Fig. 15 the simulation results of the surface pressure distributions for various building arrangements are compared, where the distance between the building models has been changed sequentially. The preasure distribution for the case of Hb is greatly different from the others. Visual Animation Numerical simulation produces a great amount of valuable information. However, no matter how valuable are the results of numerical simulation, they themselves are only a series of figures. The techniques of computer graphics are indispensable in assimilating massive amounts of information into pictures of the flowfield and in visualizing the characteristics of a turbulent flowfield. (1)Instantaneous flowfields of velocity, pressure and vorticity Figs. 16, 17, and 18 illustrate the instantaneous flowfields of velocity, pressure and vorticity respectively (Murakami [ 1988d]). Various scales of turbulent structures distributed within the separation zone and wake region are observed clearly. The distribution of the vorticity field shows structures similar to those of velocity and pressure but seems to be composed of finer structure. (2)time lines Fig. 19 illustrates the flowfields of time lines. This visualization corresponds to a wind tunnel experiment using the smoke wire method. In Fig. 19, 400 particles are arranged in each line. As the number of the particles is too large, the picture of their movements seems to be same as that of smoke diffusion. The fine structures of turbulence in the region of separation or wake are clearly expressed.

531

(1) vertical plane

(2) horizontal plane Fig. 16 Instantaneous flowfield of scalar velocity by computer graphics (14×10 4 grid points}

(i) vertical plane

(2) horizontal plane Fig. 17 Instantaneous pressure field

(i) vertical plane (I) vertical plane

(2) horizontal plane Fig.18 Instantaneous vorticity field

(2) horizontal plane Fig.19 Time lines

532

APPLICATION OF k - ¢ MODEL TO THE FLOWFIELD AROUND A BUILDING COMPLEX UNDER CONSTRUCTION Flow Around an Actual Building Complex The results of the numerical simulation of the airflow around four buildings to be located on an urban renewal site in a city located near Tokyo are presented. Figure 20 illustrates the site where four buildings are to be situated: A is a 4 storey building, B is 19 storeys, C is 20 storeys and D is 7 storeys. The buildings are to be used for a shopping center, apartments, offices and a parking garage, respectively. In this simulation, the standard numerical method based on a l - t model is used, the details of which are described in the section of "diagnostic system for numerical simulation".

(I) Preliminary study using simply-shaped building complex Before the simulation using the actual complicated building complex, a preliminary study is conducted as a first stage using a simply-shaped building complex. Fig. 21 shows a comparison of the simulation and experiment with this simply-shaped model. The correspondence between the two is very good. (2) Velocity vector field around an actual building complex Mesh dividings for the area shown in Fig. 20 are given in Fig. 22. The total number of grid points is 52(x-direction)×55( y direction)x42(z-direction) =120,120. The computational domain covers 350m in the x - direction, 370m in the y - d i r e c t i o n and 215m in the z - d i r e c t i o n . In the simulation, both the buildings on the site and the main outer surroundings have been reproduced, as can be seen in Fig. 22.

2.:

l

:.:

.

" : : ,

Ig

::::.::.:

(I) wind tunnel experiment

Fig.20 Building complex under construction and its surroundings (Kawasaki, Japan)

(2) simulation Fig.21 Velocity vector field around simply-shaped building model

533

The predicted velocity field at near ground level (z-l.5m), with the wind direction from 'the SSW, is illustrated in Fig. 23. Prevailing winds in the Kawasaki district are from the SSW. The numerical results reproduced the complicated flowfield at ground level. It is shown in Fig. 23(2) that a relatively strong wind blows into the open space between Buildings B and C. Velocity vectors at the x - x ' section and the y - y ' section are illustrated in Figs. 24 and 25. It is well recognized in these figures that a strong wind hits the windward surface of Building C and diverges'as the gust attacks the open space surrounded by the four buildings. The velocity vectors in the horizontal plane at z =41m (13th floor) are shown in Fig. 26. Strong gusts occur at the windward corners of two high-rise buildings and the space between the buildings. A reverse flow at the side wall of Building C can be seen in this figure.

•

.

I1~

•

.

.~II

I I I I i I

. . . .

i'i

/

s

L.

i

(1)vertical plane

/

* .~sr~wr.t~#

I I

I

t

i

i

tf~l~tttt~p,

I t

.t f

t

I

I I

i

s sw (1) flowfield around building complex and surroundings

II I J!T

t

tr t

Bl_.C

ah:l . . . . . I~H

liI!!,h

'~

h

SSw ,lilt

L.

~

/

(2)horizontal p!ane ~ISW

Fig.22 Mesh dividings (Region I)

1

,,

I.

4

I[

#ill%%S*ll

t

% % N

~

(2) flowfield near the building complex

Fig.23 Velocity vector field at ground level (z=l.Sm)

534

~ ~ ' ~ _ ~

X"

. ,

I

¥

V

Y' (section) Bldg.A X

i

I

"

X'

~---point ~x wind direction

wind direction

point of view Fig.24

of view

~ Bldg. ] g

Velocity v e c t o r field (vertical plane at x - x ' )

~T

Fig.25

Y Velocity v e c t o r field (vertical plane at y - y ' )

,

71 rl ~, . . . . . . . . . . . . .

[.~

pedestrian deck

Bldg.C _.

~_~

uare-shaped area at the corner of \ ~ln[dg:B indicats the Region 2 to be used ] the next stage of slmulation J

/ ssw Region

2

Fig.26

Flowfield at 13th floor of buildings B and C

535

-•

~ outdoors ~

[ ///

indoors

/// (1) vertical plane

handra i 1

indoors

/

1.3m

I

andrail

isw

Fig.27

Detailed view of balcony

(2) horizontal plane

l'i tl~,, . ~.

.~

Fig.28

I i ~~.,~.,,.~.~

! I ~t

'::

"k

Mesh dividings for Region 2

\\\\\X\\\\":~.%'~,'<,.~,,.-.~ ~ ,\ \ \ \ \\XN\Nl\N%%%%'~' ~

.

.

....

I l l flowfield st balcony corner w i t h o u t windbreak of solid fence

-- ~

: ~ . . ~ . ~

:

=

:::.:~.~[~:::

~

. . . . . . . . .

~

--

:

:

:

--: : : : :

(I) without windbreak

~llt= l|lll 'illf

)f

mdrail

.......

:..:.:,&.&,.,:,,,&,,:,, (2) with windbreak

(2)

flowfield at balcony c o r n e r w i t h w i n d b r e a k of s o l i d fence

Fig.29

Effect of windbreak on strong wind at corner of balcony (horizontal plane, z=5Ocml

Fig.30

Same as Fig.29, but vertical plane at y--y'

536

Strong Winds at the Balconies of Building B (Apartment Building) Several numerical simulations were conducted in order to investigate the characteristics of strong winds at the balconies of Building B and the shelter effects provided by the installation of a solid fence. Fig. 27 shows details of a balcony. The flowfield near the windward corner of Building B at the 13th floor (Region 2, shaded area in Fig. 26) is predicted using the mesh dividings shown in Fig. 28, which are much finer than those used in Region I. In these simulations, the technique of "adaptive grid" is utilized. The boundary conditions of the computational domain of Region 2 were interpolated with the results using coarse mesh dividing of Region I, which has a much larger computational domain. Velocity vectors around the corner are illustrated in Figs. 29 and 30, comparing the two cases with and without the windbreak. It is clear that strong winds occur at the windward corner of the balcony and that the solid fence, installed at the corner, is very useful as a windbreak. CONCLUSIONS

The current status of computational wind engineering is reviewed. (I) After examining the state of the art of numerical simulation, the following are recommended as preferable methods for application in the field of wind engineering: CVM, staggered grid, M A C method, k - ¢ model, QUICK scheme, etc. Besides the factors mentioned above, there are many others which must be studied in order to conduct an efficient and accurate simulation, for example, wall boundary conditions, method of time-advancement etc. (2) Numerical simulation cannot be entirely free from various numerical errors, and its reliability must be confirmed by means of experiments. The diagonostic system for assessing the results of numerical simulation should be developed from the viewpoint of engineering application. (3) Time-dependent flowfields given by LES and the techniques of visual animation based on them are very useful tools in turbulent flow analysis concerned with wind engineering and provide information hardly given by experimental techniques. (4) Numerical simulation can predict the 3D flowfield around a complicated building complex. By successively decreasing the computational domain, and with the aid of the adaptive grid technique, the micro-scale airflow near a building wall--for example the airflow within a balcony-- can be reproduced very' well. ACKNOWLEDGEMENTS The numerical study utilizing the k - e model was carried out in cooperation with I~,r.A. Mochida (Tokyo University) and Mr. Y. Hayashi (Tokyo University). The research concerned with LES was conducted in cooperation with Dr. A. Mochida and Mr. K. Hibi (Shimizu construction Co.,Ltd.). This report was made with the assistance of Dr. A. Mochida, Mr. Y. Hayashi and Miss H. Yasuda. The author would like to express his gratitude for their valuable contributions to this work.

537

NOMENCLATURE

f

: time-averaged or f i l t e r e d value of /

f':deviatlon from f H b : height o f building model h ~ : minimum mesh interval adjacent to solid wall ~ :mesh interval in ~ d i r e c t i o n A : g r i d scale (A=(AI'Az-AJ) ~ / ' ) : length scale of turbulence u~ : i = l ( s t r e a m w l s e ) , i = 2 ( s p a n w i s e ) , i = 3 ( v e r t i c a l ) , three components of the v e l o c i t y vector in the z~ d i r e c t i o n (.Tb : v e l o c i t y at inflow of computational domain at height of H b P :pressure P o : reference s t a t i c pressure Cp : pressure c o e f f i c i e n t , Cp=(P-Po)/(),~vU6 ~) k

: turbulent k i n e t i c energy,

k-~.~

k° k~ P~ ¢ e~

:$GS (subgrid s c a l e ) : t u r b u l e n t kinetic energy : k value at the f i r s t grid point near solid wall :Production of k : turbulence d i s s i p a t i o n rate : e value at the f i r s t grid point near solid wall : kinetic viscosity : eddy v i s c o s i t y • s , : $GS eddy v i s c o s i t y

u ° : f r i c t i o n velocity (u,~ = r~ =C i/,k, ) p ( r ) w : s h e a r s t r e s s at the wall Re

: Reynolds number ( R e = - ~ - - ~ )

~ ' h : s o l u t i o n e r r o r estimated by c a l c u l a t i o n of mesh size h ~ h : t r u n c a t i o n e r r o r of mesh size h All p r o p e r t i e s presented here are made non-dimensional by Hb and Ub.

REFERENCE Baetke. F., Werner, H. and Wenglem. H. 1987, "Computation of Turbulent Flow around a Cube on s Vector Computer', Prc.6th Syrup. on Turbulent Shear Flows. Caruso, S.C.. Ferziger, J,H. snd Uh~ger, J. 1986, "Adaptive ~3rld Techniques for Elliptic Fluid-flow problems', AIAA 24th Aerospace Science Meeting. Desrdrorff, J.W. 1970, "A Numerical Study ol Three-dimensional Turbulent Channel Flow at large Reynolds Numbers', ~.Fluid Mech., Vol.41, pp.453-480. Hariow, F.H., and Welch, J.E. 1965, "Numericsl Calculation of Time-dependent Viscous Incompresible Flow of Fluid with Free Surface', Phys.Fluids., Vol.8, pp.2182-2189. Horiutl. K. 1982, "Study of Incompressible Turbulent Channel Flow by Large Eddy Simulation', Theoretical and Applied Mechxnics~ VoL31, pp.402-407. Isbizu, Y. 1986, "Evaluation of Ventilation System through Three-dimensional Numerical Simulation', Trans.ASHRAE., Vol.30, pp.l-7 Kim, ]., Moin, P. and Moser, R. i987, "Turbulence Statistics in Fully Developed Channel Flow at Low Reynolds Number', ],Fluid.Mech.~ Voi.177 Launder. B.F- and Spalding. D.B, 1972, Mathematical Models of Turbulence~ Academic Press. Leonard, B.P. 1979, "A Stable and Accurate Convective Modelling Procedure based on Quadratic Upstream interpolation', Computer Methods in Applied Mechanics and Engineering, VoI.19, pp.59-99

538 Mathews, E.H. 1987a, "Prediction of the Wind-generated Pressure Distribution around Buildings", J.Wind Eng.lnd.Aerodyn.. Voi.25, pp.219-228 Mathew, E.ll. a n d ' q ~ l ~ y e ~ - - ~ . l ~ - - I ~ m p u t a t i o n of Wind Loads on a Semicircular Greenhouse around Buildings", Preprints.7tb Int.Conf.on Wind Eng., Vol. 14, pp81-89 Moin. P. and Kim, J. 1982, "Numerical investigation of Turbulent Channel Flow', J,Fluid.Mech., Vol.118, pp.341-377 " " Murakaml,~hida, A. and Hibi, K. 1983, "Numerical Simulation of Air Flow around Building IPartl-Part2)", Summaries of Technical Papers of Annual Meeting Architectural Institute of Japan, pp.487-490 Murakami, 5., Mochida, A. and Ilibi, K. 1984, "Numerical Simulation of Air Flow around Building by means of Large Eddy Simulation (Partl~Part2)", Summaries of Technical Papers of Annual Meeting Architectural Institute of ~apen, pp.257-260 Murakami, S., Hibi, K. and Mochida~. A. 1985, "Visualization o1 Computer-generated Turbulent FIowfield around Cubic Model , Proclnt.Symp. on Fluid Control and Measurement Murakami, S., aod Kato, S. 1986, "Discussion", Trans.ASHRAE., Vol.32, pp. ll5-117 Murakami, S., Mochida, A. and Hibi, K. 1987a, "Three-dimensional Numerical Simulation of Air Flow around a Cubic Model by means of Large Eddy Simulation", J.Wind Eng,lnd.Aerodyn., ¥ol. 25, pp.291-305 Murakami, S., Mochlda, A. and Hibi, K. 1987b, "Numerical Prediction of Velocity and Pressure Field around Building Model", Proc.7th lnt.Conf.on Wind Eng., Vol.2, pp.31-40 Murakami, S., Kato, S. and Nagano, S. IggRa, "Estimation of Error Caused by Coarseness of Finite-differencing", Journal of Architecture, Plamdug and Environmental Engineering (Trans.AIJ.I, No.385, pp.9-17 Mm'ak~ochida, A.. Oowada, J. and Hayashi, Y. 1988b, "Three Dimensional Numerical Simulation of Turbulent Flowfield around Building by means of k-= 2-equation Model (Part I )", Journal of Architecture, Planning and Environmental Engineering (Trans.AIJ.), No.392,pp. I I-21. Murak~ochida A. and Murakami, L. 1988c, "Numerical Simulation of Turbulent FIowfield around a Building with Adaptive Grid Technique--Study on passive method applied to 3D k-¢ two equation model", Journal of Architecture, Planning and Environmental Engineering (Trans.A1J.), No.~93, pp.l-9. Murakami. S. 1988d, "Visualization of Turbulent Flowfield Generated by Numerical Simulation", Proc.3rd lnt.Symp, on Refined Row Modelling and Turbulence Measurement, pp.I13-124. Murakami. S. and Mochida, A. 1988e, "3D Numerical Simulation of Airflow around a Cubic Model by means of k-= Model", J.Wind Eng.lnd.Aerodyn., Vol.31. No.2. Murakami, S. and Mochida, A. 1988f, "30 Numerical Yrediction of Turbulent Flow around Buildings by means of k-= Model", Building and Environment, VoL24, No.l. Murakami, S. and Kato, S. and Kondo, Y. 1988g, "Simulation ot Airflow by means of Reynolds Stress Model (Partl~Part2)", Summaries of Technical Papers of Annual Meeting Arcbitectual institute of Japan Murakami, S. 1988h, "Numerical Simulation of Turbulent Flowfield around Cubic Model", , No.37, pp.239-252.(Proc.Int.Colloquim on Bluff Body Aerodynamics and Its Kyoto, 1988) ' Nomura, T.. I~mtsuo, Y. and Kato, S. 1980, "Study on Computational Mesh of the MAC Method", Trans.AIJ., No.292, pp.61-72 Paters~and Apelt, C.J. 1986, "Computation of Wind Flows over Three-dimensional Buildings', J.Wind Eng.lnd.Aerodyn., Vo1.24, pp. 192-213 Smagorinsky, J.S. 1963, "General Circulation Experiments with the Primitive Equations : Part I, Basic Experiments', Monthly Weather Review, Vol.91, pp.99-164 Schumann, U. 1975, "Subgrid Scale "Model for Finite Difference Simulation of Turbulent Flows in Plane Channels and Annuli", ~ , Voi. 18, pp.376-404 Tamura, I., Kawahara, K. and S h i r a y a ~ "Numerical Study of Unsteady Flow Patterns and Pressure Distributions on a Rectangular Cylinder', Proc.Tth lnt.Conf.on Wind Eng., Vol.2, pp.41-50 Viecelli, J.A. 1971, "A Computing Method for Incompressible Flows Bounded by Moving Walls", J.Comput.Phys., Vol.8, pp. 119-143 Yeung, P.K. and Kot, S.C. 1985, "Computation of Turbulent Flows past Arbitary Twodimensional Surface-mounted Obstructions', J.Wind Eng.lnd.Aerodyn., Vol. 18, pp. ll7-190

~