In this paper, numerical studies on the flow and diffusion fields around a building and cooling tower, and on the unstably-stratified flow over heated urban areas ...
NUMERICAL SIMULATION OF TURBULENT DIFFUSION IN CITIES
S. MURAKAMI, A. MOCHIDA ANDY. TOMINAGA
Institute of Industrial Science University of Tokyo 7-22-1, Roppongi, Minato-ku, Tokyo 106, Japan
ABSTRACT. Examples of numerical simulations on turbulent diffusion fields in cities are presented including comparison with wind tunnel experiments. The first example is the computation of the flow and diffusion fields around a cooling tower by k-E turbulence model. A composite grid technique, which connects two grids supported by a fortified solution algorithm, is adopted in this computation. The second example is the computations of flow over a heated urban area by Large Eddy Simulation (LES) and k-E model. In the k-E model computations, two different modellings for buoyancy effect are adopted; first is the model proposed by Rodi (1984), and second is the one proposed by Viollet (1987). The last example is the LES computations of turbulent diffusion of buoyant (lighter-than-air) and heavy gases discharged in the wake region behind a building model. Results given from two types of Smagorinsky models, i.e., the standard one and its modified version proposed by Mason et al. (1989, 1990), are compared here. In the concluding remarks, the relative abilities of various turbulence models are discussed.
1. Introduction The micro-climate around buildings is most important from the viewpoint of air quality. Daily life produces all kinds of gaseous contaminants which are discharged into the outside air around buildings ; the air quality in this space is thus highly dependent on the diffusion mechanism by which this contaminated outside air is convected and diluted into the atmosphere. Hence, it is necessary to clarify the structures of the flow and diffusion fields around buildings in order to predict the air quality in out-door living spaces. The flow and diffusion fields around buildings are very complicated and three-dimensional. The analysis of such fields was usually performed by means of wind tunnel tests, but in recent years, methods of numerical simulation have developed rapidly as a new analytic tool (Murakami (1990c, 199lb, 1992c)). In this paper, numerical studies on the flow and diffusion fields around a building and cooling tower, and on the unstably-stratified flow over heated urban areas are presented.
2. Brief Review of Turbulence Models 2.1. KOLMOGOLOV MICROSCALE AND GRID DISCRETIZATION
The flowfields treated in civil engineering are usually expressed by the incompressible Navier-Stokes equations (hereafter N-S equations). This set of nonlinear partial differential equations is infamous for the difficulty with which they are solved. Furthermore, the flows treated in this field are always highly turbulent, and we must thus rely on the theory of turbulence statistics when solving these N-S equations by means of numerical methods. 681 J. E. Cerm11k et al. (eds.), Wind Climate in Cities, 681-701.
© 1995 Kluwer Academic Publishers.
682
Turbulent flowfields are characterized by various scales of turbulent eddies. Hence the most important decision in numerical simulation is how finely we should resolve these eddies in numerical analysis. The amount of CPU time is completely dependent on this judgment. Resolving small eddies requires much CPU time. The smallest eddy in turbulent flow, the socalled "Kolmogorov microscale 11" is very small, around 0.1-l mm for natural wind. It will be evident that discretization of the flowfield corresponding to the Kolmogorov microscale 11 is very difficult in practice. Numerical analysis of a flowfield, whether the numerical method be finite difference or finite element, begins with the discretization of the flowfield; computation is conducted using these discretized expressions in the N-S equations etc. While the amount of CPU time increases with greater numbers of these discretizations of space and time, the accuracy of the prediction is governed by the fineness of discretization. As was stated above, turbulent flow can be reproduced with complete accuracy only if the flowfield is discretized corresponding to the Kolmogorov microscale. However, such a fine discretization is almost impossible to obtain in practice. The size of the microscale(11/l, l : integral scale) decreases as the Reynolds number (=ul /v, hereafter Re, u:characteristic velocity scale) increases. Consequently the number of grid points must be increased in proportion to the value of Re9/4, according to the theory of turbulence statistics. The numerical method based on this concept of microscale discretization is called "Direct Numerical Simulation" (hereafter DNS, Kim et al.(l987)). While this simple and clear method uses no turbulence modelling, it requires a huge number of grid points if the Reynolds number of the flowfield is large. The Reynolds number of the flowfields treated in civil engineering is about 1QL 107, so DNS of such a flowfield requires around 10 1L 10 16 grid points. On the other hand, the maximum number of grid points ever used is around JQ6-JQ7 (Le et al.(l993)). Therefore it can easily be understood why it is practically impossible to solve a fully turbulent flow by means of DNS even if computer performance is advanced greatly in the future. 2.2. TURBULENCE MODELS- CURRENT STATUS AND FUTURE POTENTIAL
The numerical analysis of turbulent flow faces a serious obstacle of the Kolmogorov microscale, the concept of which comes from the theory of turbulence statistics. In order to overcome this obstacle, various types of turbulence models have been invented. Turbulence modelling is designed to simulate the averaged flowfield (coarse graining, cf. Fig.l) rather than the original flowfield. Such small scale eddies, which are difficult to resolve, are thus neglected in the process of coarse graining. There are two types of coarse graining ; ensemble (or time) averaging and space averaging. These types of averaging and turbulence models are illustrated in Fig.2, including DNS, which uses no turbulence model. By averaging a flowfield with steep fluctuations, the averaged fluctuations of the flow variables become less steep and the relatively coarse grid discretization can resolve the flowfield with a certain accuracy. The small eddies neglected by coarse graining can be incorporated in the simulation through the turbulence modelling. In this section, the concept of turbulence modelling is described briefly. 2.2.1. Turbulence Models Based on Reynolds-Averaged Equations. The ensemble (or time) averaged forms of the N-S equations are called the Reynolds equations. These equations express only the movements of large scale eddies. Since the flowfield contains only larger eddies, numerical analysis can be conducted with a rather coarse grid discretization. Most CFD simulations in the field of engineering have been based on this type of modelling partly because of inadequate computer performance. The most famous and most widely-used is "the k-£ two equation turbulence model" (hereafter k-£ model, Launder and Spalding( 1972)). However, it has become clear that the k-£ model possesses some shortcomings. The high-precision version of the Reynolds-averaging approach is the Differential Second-moment Closure ( hereafter
683
DSM). Applications of DSM and its simplified version, ASM (Algebraic Second-moment Closure Model, Launder, Reece Rodi (1975), Rodi (1976), Murakami et al. (1990a, 199la, 1993)) are increasing year by year. 2.2.2. Turbulence Models Based on Space-Filtering. Another method for coarse graining is "space-filtering", the concept of which is illustrated in Fig.!. The degree of averaging can be adjusted easily by changing the width of the filter. Since eddies smaller than the width of the filter(eddies of subgrid scale) are not analyzed directly and eddies larger than the filter width(eddies of resolvable or grid scale) are simulated, this type of modelling is called "Large Eddy Simulation" (hereafter LES, Smagorinsky(1963), Deardorff{l970), Schumann{l975), Moin et al.(1982), Murakami et al.(1987), Germano(1992)). In LES, the effect of small eddies neglected by filtering is incorporated in the simulation through the turbulence model of subgrid scale (SGS) Reynolds stress. For modelling the SGS Reynolds stress, the SGS eddy viscosity model proposed by Smagorinsky has been most commonly used (Smagorinsky(l963)). Since Reynolds-averaging is not applied to the basic equations used in LES, the scale of resolved eddies in LES is much smaller than those of the k-e model or DSM. Thus LES predicts a timedependent flowfield with small fluctuations. Consequently the effect of modelling in subgrid scale on the turbulence structure of the flowfield is relatively small, compared to that of k-e. In fact, LES has succeeded in reproducing the properties of a highly anisotropic flowfield (Murakami et al. (1990b, 199Ia, 1992a, 1992c)).
!.]
space
Fig. 1 Filtering of flowfield
I
Navier-Stokes equations
l
I
space av~aging (filtenng)
JI
···-·········· ································· ...........
l
time averaging ~ J (ensemble averaging)
I I
direct numerical( Integration
I
Direct numerical simulation
.............. ...............1 Large eddy simulation !························
I _
k- • moo a Differential secondmoment closure model
I
length
energy dissipation scale ( n : kotmogorov scale)
(larger than n , smaller than l
)
1·....................................... characte~stlc length scale tor turbulence ( t : Integral scala) characteristic length scale for mean ftow (L: channel width, building height etc.)
Fig. 2 Methods of numerical simulation and turbulence models
scale
l
684
In recent years, more refined SGS models for LES have been developed in the field of meteorology and mechanical engineering, for example, extended versions of the standard Smagorinsky model (Mason (1989), Mason and Derbyshire (1991), Yoshizawa (1991 )). new versions of the SGS eddy viscosity models (Gemano (1991), Horiuti (1993)), the SGS oneequation model (Deadorff(1980), Moeng(\984), Yoshizawa and Horiuti (1985), Nieuwstadt et al. (1991)) or the SGS second-moment closure model (Dcadorff (1973, 1974), Schmidt and Schumann ( 1989)). To sum up these various surveys, LES is most promising for the simulation of the flowfields treated in wind engineering. The weak point of LES is the large amount of CPU time required for computation, which is much larger than that fork-£. While the k-£ model gives an ensemble or time-averaged flowfield, LES must reproduce and store a long sequence of time-dependent flowfields. However this weak point will be overcome to a certain degree by future growth in computer hardware technology. In the following chapters, examples of numerical simulations by LES and k-£ model are presented.
3. Computation of Flow and Diffusion Fields around Cylindrical Cooling Tower with Composite Grid Technique The first topic to be discussed here is the analysis of the flow and diffusion fields around a cylindrical cooling tower with a composite grid technique. The same model was computed by Majumdar and Rodi using the ordinary polar cylindrical grid , and also by Demuren and Rodi using the ordinary rectangular grid. Here the standard k-£ model is used for turbulence modeling as it was in the previous work (Majumdar and Rodi (1989), Demuren and Rodi( 1987) ). 3.1. OUTLINE OF COMPUTATION
A polar cylindrical grid is very suitable for the accurate definition of the boundary conditions around cylinder walls, and in this regard it is much better than the rectangular grid. However, the disadvantage of the polar grid was also pointed out by Majumdar and Rodi (1989) as follows: the spacing of the polar grid becomes rather coarse and cannot be aligned with the streamline for regions far from the model. In the present computation, a composite grid technique, which connects two or more grids supported by a fortified solution algorithm, is adopted. Details of the composite grid technique adopted here are given by Ishida et al. (1993) and Mochida et al.(1993). Fig. 3 illustrates the grid layout. A rectangular grid (grid B) is used to cover the whole computational domain and a polar cylindrical grid (grid A) is applied to the region near the cooling tower. Thus grid A overlaps grid B. It is well known that computation in the vicinity of the pole is very difficult with a polar cylindrical grid since the radius approaches zero there. Thus, polar grid A does not
Fig. 3 Grid laymg ( Grid A is overlapped onto grid B )
685
inch.:de the region near the pole.The region near the pole is only covered by rectangular grid B, as shown in Figs. 3. A blending technique is adopted over the overlapping regions in order to connect smoothly the dependent variables on rectangular grid B and cylindrical grid A (Delsem et al. (1986), Fujii et al.(l988), Ishida et al. (1993), Mochida et al. ( 1993)). The boundary conditions used in this case are illustrated in Table l. The concentration field is treated as passive scalar and the value 0.5 is adopted for the turbulent Prandtl number (Majumdar and Rodi (1989)). The QUICK scheme is applied to the convection terms of all governing equations. For the time advancement, the first order fully implicit scheme is chosen.
_. ---4>
6
5 4
3 2 '1
-3 -2 -1
0
2
3
6
(1) Measurements (Viollet)
3.2. RESULTS AND DISCUSSION
The distributions of the mean velocity vectors and mean concentration are illustrated in Fig. 4. The results of the present method are compared with those of the measurements by Viollet and the computation of Majumdar and Rodi (1989). It is shown in Fig. 4 (3) that a smooth connection is given by the present method for the distribution with the aid of the blending technique. The distribution of mean velocity vectors predicted by the present computation agrees rather well with the measurements (Fig. 4 (1)). The distributions of mean velocity and given by the present method are observed to be reasonable. Table l Boundary conditions
inflow
outflow
=xf"', ==O k(J:,)=0.0025 i(J:,J=, not only on the mean momentum equations but also on the equations of k and E, as is summarized in Table 3. The treatment of Gk ( buoyancy production ) in e equation is different in the two models. The shear layer treated here is horizontal and unstable. The coefficient C3 for Gk in e equation is greatly different for both models : C3=1.44 for Viollet model (eq. (14)) and 0.29 for Rodi model (eq. (15)). In the three-dimensional LES computation, the computational domain covered 8.0Lo (x 1-direction), 2.0Lo (x2-direction) and 2.0Lo(x3-direction). This domain was discretized into 53(xi)x21(x2)x17(x3) meshes. In the computations of the k-E model, only the XJ-X3 plane was treated. The size and mesh spacing of the computational domain (XJ-X3 plane) in the k-e model computations were the same as those in the LES computation. Boundary conditions are described in Table 4. A second-order centered difference scheme was used for all spatial derivatives in the case of LES. The QUICK scheme was used for the convection terms in the cases of the k-E model. For time advancement, the Adams-Bashforth scheme was used for the convection terms, the Crank-Nicolson scheme for the diffusion terms for LES, and the fully implicit scheme for the k-E model. inflow (lineA)
lineB
lineC
lineD
lineE Roughness
x, '-....Hot panel
1.25L 0
0.875L 0 0.875L 0
1.3L 0
~---------8.0L 0,--------------+I
Fig. 6 Model flow field of heated urban area
688 Table 2 Model equations for LES for non-isothermal flow ( standard Smagorinsky model)
au; au;u; a (- 2 ) a au; au 1 -8t +-8- = -8 - p+0Xj +3 ksas +8Xi (v+vsas)(8X, )-Arileo13 _ Xi X;
•••
(1)
••• (2)
OU;=O OX;
aeat+ aeu; a( 1 ae) axj = OXj_}RePr ~sd OXj
••• (3)
••• (5)
llsas=(Csh)2{8u;(8u;+8u;)}"2 ••• (4) Ar=gf3Lo OXt 8x1 OX; ' 2 llsas asas=p-, Prsas=O.S • • • (6) Cs = 0.1 Tsas
Table 3 Model equations for k-e model for non-isothermal flow ( Viollet model and Rodi model)
a ( 2 ) a a a ••• (7) aX; P+-3 k +8Xj-v,-Aro,. X, 1
... (8) ... (9)
... (10)
Type 1 (Viollet model) :
ae 1>-=ae a [lie ae] e e2 -+O (unstable) (present computation) : C3=C 1=1.44
• d.;:;;o
(stable) : Cs=O
Type 2 (Rodi model) :
aeat +ae a-.[ll, de] e ' e - +C,-k (P.+G.)(1+C3'R!J-C,-k 0x =8x, a, 8x, 2
1
... (15)
Cs'=O.B ·Vertical shear layer:
RJ=O (Cs=C,=1.44)
• Horizontal shear layer (present computation)
·. RJ=-_0__ P.+G•
(C3=0.29)
689 4.3. RESULTS AND DISCUSSION
4.3.1. Mean Velocity Distribution (Fig. 7). The mean velocity profiles given from the experiment, LES and the two versions of the k-e model are almost identical. 4.3.2. Mean Temperature Distribution (Fig. 8). The result of LES based on the standard type of Smagorinsky model shows very good agreement with the experiment for the flowfield analyzed here, although the buoyancy production of SGS component of turbulent kinetic energy ksas is not considered in this model (cf. the Appendix ). On the other hand, distribution as is shown in Fig. 8. The mechanism which causes the difference in the results from the Violtet and Rodi models is discussed in section 4.3.5. Table 4 Boundary conditions LES
k-•
-rrot. (measured profile in wind tunnel) = 0, k(IJ : (measured value in wind tunnel),
=xj'2, =0, =0
=O v,(I,) = k(x,)'"f(IJ, t(rJ = C,k(I,)"'!e(:t,) err,)= -k(I,)'"'(o/ox,)'' (-/()
E
Fig. 11 Profiles of k ( k-E) 4.3.3. Distributions of Velocity Fluctuation 112 in the cases of the isothermal and non-isothermal flowfields measured in the wind tunnel. In the case of the non-isothermal flowfield, the positions of peak value of ( -2Ar). The results of LES reproduce this tendency rather well (Fig. 9 (3)). The result given by the Rodi type k-E model also reproduces the increase of 1/2 cannot be reproduced well by the Viollet type k-E model (Fig. 9 (2)).
691
4.3.4. Distributions of (Fig. 10). Fig. 10 indicates that the distribution of obtained by LES is generally less steep in comparison with those given by the two types of k-E model. Judging from the better agreement of with the experiment in the LES result (Fig. 8), the distribution of in LES seems to be more reasonable than those of the k-E model computations. The most important point to be noted in the result of LES is that the positions of the peak value of move to the upper region as the flow goes downstream. This tendency is reproduced rather well by the Rodi model. On the other hand, the result of the Viollet model shows a different distribution, with the peak positions observed only in the area X3/Lo/2 (PsiPa=l.O)
( 1) experiment
(2) Case 2
Fig. 14 / (pJpa=l.7)
(2) Case 2
Fig. 17 / 2 (p 5 /pa=0.3)
The results of two types of Smagorinsky model, the standard ( type 1) and its extended version ( type2), are compared with the experiment for the case of p siP a= 0.3 ( Fig. 15). The result of the experiment shows that decreases due to the buoyancy effect in the area near the ground in the case of PsiPa = 0.3 (Fig. 15 (1)) in comparison with the case of PsiPa = 1.0 (Fig. 13 (1)). It is shown in Fig. 15 (2) that this tendency is not sufficiently reproduced in the result of Case3 (type I) and is relatively large in this area as well as in the case of PsiPa=l.O. As is illustrated in Figs. 15 (3) and 19, this discrepancy is improved to some extent in the result from the extenlkd version of Smagorinsky model (type2), although value is still overestimated at the center section. The difference observed in the results of Cases 3 and 4 indicates the effectiveness of the extended SGS models as is di~cussecl further in sections 5.3.3 and 5.3.4.
693
Table 5 Model equations of LES with buoyancy effects 00 '=0 ••• (22)
ax,
(au; au;)] a(- 2 ) a [ 1 au;+au;u;_ ax, --ax,P+]kscs+ax, (Re+vSll,) ax/ax, -Frdco., at
l ...
(23)
(24) ac +aU; c = j__[(-1-+a ) ac at ax, ax, ReSc ''"ax, ••• (25) a,,=vSlJJScSlls (ScSlls=0.5) Frd =-Lip H, g ••• (26) p, ' (a) Type 1 (the standard type of the Smagorinsky model)
(au; au )']'/'
• • • (28) liscs= (cS f , h) ,S • • • (27) S= [ -1 -+~ 2 ax, ax, ••• (29) !.= l-e:rp(-y"/25) Cs=O.l2, (b) Type 2 (the extended version of the Smagorinsky model) • • • (30) '.Jscs=(Cs f,h)' 11.81) ••• (38) (u.: friction velocity)
u.
Table 7 Cases analyzed in Chapter 5 p./p,
easel Case2 (heavy gas) Case3 (buoyant gas) Case4 (buoyant gas)
Frd
u,.,
c,
1.0
0.0
eq.(27) (Typel)
0.12
1.7
+8.6
eq.(27) (Typel)
0.12
0.3
-8.6
eq. (27) (Typel)
0.12
eq.(30) (Type2)
0.12 (neutral, stable) 0.32 (unstable)
0.3
-8.6
Table 7 lists the four cases analyzed here. Three different values of the densimetric Froude Number (Frd=-g~pH/pa2) were specified : -8.6 for cases 3, 4; 0 for case 1 ; 8.6 for case 2. The Standard Smagorinsky model (Type!) was applied in Cases 1-3 and the value of 0.12 was selected for the Smagorinsky constant Cs. The extended version of the Smagorinsky model (Type2) was applied to the buoyant gas case ( Case4, PsiPa = 0.3, Frd=-8.6). The results of Case4 (type2) were compared with those of Case3 (type!). It should be noted that the stability of the flowfield changes locally in the case of buoyant (lighter-than-air) gas emission ( Cases 3 and 4). The value of Cs was switched locally in accordance with the stability of the flowfield in Case 4 (type2) : Cs=0.12 in the neutral (Rf=O) and stable (Rf>O) regions and Cs=0.32 in the unstable (Rf (Figs. 16-18 and 20). Figs. 16- 18 illustrate the distributions of are almost exactly similar to these of for each case of PsiPa· The results given from type! show fairly good agreement with the experiment for all cases except for the discrepancy observed in the values of given from the experiment decreases rapidly as the measuring point approaches ground level (Fig. 20), while observed in the experiment (Fig. 20) 5.3.3. Distributions ofTime-Averaged Total Kinetic Energy of Velocity Fluctuations ( Periodic +Stochastic) ( Fi!f. 2 I). The result of the experiment shows a maximum peak of at x3::0.3H just behind the building model, as is shown in Fig. 2 1. The height of this peak is located in the high concentration region where /;:: 10, as is shown in Fig. 19.
696
---,~
Fig. 23
Fig. 22 distributions (case 4)
«I» distributions (case 4)
1.2
1.0.
x/H
0.8 .. hne
0.6 0.4
0 /
Fig. 24 Lateral profiles of /
-0.4
Fig. 25 Vertical profiles of /
The distribution of in Case3 (type!) is quite different from that of the experiment (Fig. 21 (1)), while the result in Case4 (type2) corresponds fairly well with the experiment (Fig. 21 (2)). As is shown in Fig. 22, the time-averaged value of the flux Richardson number shows negative values in the area where the maximum peak is observed. This means that the flowfield is unstable and the increase of is mainly caused by the buoyancy production of stochastic turbulent kinetic energy k, as denoted below. Fig. 23 displays the distribution of in Case4(type2), which was calculated by eqs. (31)-(33) in Table 5 as a function of the flux Richardson number Rf. The distribution of «I» shows a maximum peak of about 1.2 in the unstable flow region just above the gas exit. Furthermore, the value of Cs is set at a large value, 0.32, in the unstable flow region in type2 (cf. Table 7). Therefore, the value of VsGs is estimated to be large in type2 in comparison with that in type 1. The large value of v SGS in type2 apparently reduces the grid scale (GS) component of stochastic turbulent kinetic energy k, but increase the subgrid component of k, ksGS (=VsGs 2 /(Ckh) 2 ,Ck=0.094). In Case4 (type2), the influence of the increase of ksGs is larger than the influence of the decrease of the GS component of k. Hence, the time-averaged total kinetic energy of velocity fluctuations (periodic+stochastic) also increases due to the increase of the stochastic turbulent kinetic energy k in Case4 (type2, Fig. 21 (2)), and the result given by type2 corresponds much better with the experiment than does type 1. Fig. 24 compares lateral profiles of behind the building at XJ/H:=0.25. The value of given by type2 becomes smaller than that in typel near the centerline. This is because the highly contaminated air in the area 0.2::>x3/H::>0.4 at the center section was well diffused in the lateral direction in the result of type2 by the diffusion effect caused by the larger value of in this area. By this decrease of at the center region, the total amount of the advection and diffusion of highly contaminated air from the area where shows a maximum peak (x3/H:=0.2) (cf. Fig. 19) toward the area near ground is also decreased for type2, and thus the distribution of from type2 corresponds better with the experiment than does the result from type!, as is shown in Fig. 19. The difference in the between types 1 and 2 is closely related to the difference in the distribution, as is discussed futher in the next section.
697 5.3.4. Distributions of (Fig. 25). Vertical profiles of behind the model are compared in Fig. 25. The distribution of given by type2 shows a smaller value in comparison with that of type! near the ground. In particular, given by type2 becomes negative in the vicinity of the ground surface (Fig. 25), where d/dX3>0 ( PsiPa) is larger than that of the negative value of the underlined terms in eq. (39) in the result of type l. On the other hand, as is shown in Fig. 25, given by type2 has a negative value near the ground, where the positive value of the buoyancy production term (-Frd), the absolute value of which is propotional to increases as the measuring point approaches the ground in type!, while ) also increases in type! and decreases in type2 in the area near the ground because Frd is negative ( -8.6) in this case. The modelling for the buoyancy effects on SGS turbulent kinetic energy, ksas. does not have direct influence on . However, the production term of is also greatly affected by the modelling for the buoyancy effects, as is shown in Fig. 20, through the influence of the modelling on the values of , , etc.
698
6. Concluding Remarks The technique of numerical simulation has developed greatly in these past ten years. However, it is not frequently applied to the problems of turbulent diffusion in cities. There exist many difficulties if we want to simulate the contaminant diffusion around buildings. These are mainly caused by the following two factors: 1) the flowfield is very complicated ; 2) the size of the space is large. The flowfield in cities is characterized by the following factors : 1) an elliptic flow composed of separation, reattachment, circulation etc., since it is usually a flowfield around bluff bodies ; 2) buoyant gas flow or heavy gas flow ; those are often thermally stratified. When the space is large, we are forced to conduct a simulation with coarse dividing mesh. The coarse mesh always causes various types of numerical errors. In order to analyze such a complicated flowfield, we should select the proper turbulence model. The relative abilities of various turbulence models are compared in Table 8. An improvement in prediction accuracy usually requires an increase of CPU time and calculation instability. Thus a strategy for selecting a turbulence model may be proposed according to the purpose of the analysis and the accuracy required. The flow pattern around a bluff body can be predicted by k-E with a certain accuracy. However, for reproducing surface pressure and turbulence statistics such as , analysis by DSM (or ASM) or LES is effective. LES provides greater estimation accuracy than the other turbulence models but requires much more CPU time. In this paper, we have presented some examples of numerical simulation of turbulent diffusion, including comparisons with wind tunnel experiments. The correspondence between the simulation results and the experiments is not perfect, but we do think this technique very promising from the view-point of practical application in the future, when the technique of numerical simulation will become a powerful tool for the analysis of wind climate in cities. Table 8 Relative comparison of various turbulence models for practical modelling of flowfields related turbulent diffusion TurbuleJ'lce model
Standard
Low-Re. No.
k-•
k-•
DSM (ASM)
Low-Re. No. DSM (ASM)
LES
0
0
0
0
0
0
0 0
0 0
0 0
0 0
Standard
------ ------- ------ ------- -------------------------------wall function non-slip wall function non-slip wall function wall boundary condition 1. Simple flows (channel flow, PIP9 flow, etc) 2.
(local equilibrium is valid) low With streamhne curvature 1) weak curvature 2) strong curvature
(flow around bluff body)
0 X. c,
0 X'
c,
LES
------non- slip
3. Impinging flow
X,t>
)(. c,
L>.O
6.0
0
0
4. Flowfield with low Reynolds No.
X ,6
L>.O
X • c,
L>.O
0
0
0
0 0
0 0
0 0
0 0
0
X.
0
X • c,
0
X 0
X 0
0 0
0 0
0 0 0 X X 0 . ' l:l : tnsuffiCtently functional
0 0
0 0
0 0
5. Non - tsothermal flow
1)
0
weak stratification
.c X • c, X
2) strong stratification 6. Convective heat transfer at wall
x.c
c,
7. Unsteady flow, unsteady diffusion
a.
1)
highly unsteady
2)
vortex shedding
1)
normal
2)
swirl
X X
X
6
Jet
0 : functtons well
X . funct•ons poorly
699
APPENDIX
DERIVATION OF THE EXTENDED VERSION OF SMAGORINSKY MODEL INCORPORATING SUBGRID SCALE BUOYANCY EFFECT (Mason(1989), Mason and Derbyshire( 1990)) The transport equation of SGS component of turbulent kinetic energy for the flowfield with the buoyancy effects is expressed as follows. (40) Dksos/Dt = Dksos + Pksos + ~sos -Esos where Pksos :shear production of ksos (=VsosS 2), Gksos : buoyant production of ksos , Dksos : diffusion term of ksos , Esos : SGS energy dissipation rate. Here, the local equilibrium for is assumed. (41) Esas = Pksos + Gksas =VsasS 2(1-Rf) where Rf = - Gksos I Pksos . It is assumed here that the major portion of the total turbulence kinetic energy k (subgrid+resolvable) dissipates in the subgrid scale. Esos =E (42) Dimensional analysis gives Vsos= C1 1 4/3 Esasl/3 (43) here l means the SGS length scale. The extended Smagorinsky model is given from eqs. (41) and (43)as follows. Vsos= (Cs l ) 2(1-Rf)l/2-s , (44) where Cs=c 13/4. The standard Smagorinsky model can be obtained in the same manner as that for eq. (44) by neglecting Gksos in eq. (41). vsos= ccs t )2S (45) For the analysis of the flowfield presented in Chapter 5, the SGS length scale was set equal to the mesh scale h ,which was multiplied by the Van Driest type wall damping factor f 11 in order to account for the near wall effect. (46) For the unstable condition (Rf
