Large Eddy Simulation of Unstably Stratified Turbulent Flow over Urban-Like Building Arrays

Abstract

Thermal instability induced by solar radiation is the most common condition of urban atmosphere in daytime. Compared to researches under neutral conditions, only a few numerical works studied the unstable urban boundary layer and the effect of buoyancy force is unclear. In this paper, unstably stratified turbulent boundary layer flow over three-dimensional urban-like building arrays with ground heating is simulated. Large eddy simulation is applied to capture main turbulence structures and the effect of buoyancy force on turbulence can be investigated. Lagrangian dynamic subgrid scale model is used for complex flow together with a wall function, taking into account the large pressure gradient near buildings. The numerical model and method are verified with the results measured in wind tunnel experiment. The simulated results satisfy well with the experiment in mean velocity and temperature, as well as turbulent intensities. Mean flow structure inside canopy layer varies with thermal instability, while no large secondary vortex is observed. Turbulent intensities are enhanced, as buoyancy force contributes to the production of turbulent kinetic energy.

1. Introduction

Global urbanization over the past decades has significantly changed the atmospheric environment of urban area. Urban environment problems caused by human activities arise and related researches become increasingly popular. Urban boundary layer (UBL) involves multiscale, multiphysics processes and thermal dynamics is one of the critical issues. Heterogeneous thermal distribution caused by solar radiation, anthropogenic heat, and building materials influences pedestrian comfort. Urban heat island [1] is induced by different heat capacities between urban and rural areas, and effect of buoyancy force on pollutant dispersion inside canopy layer is considerable. Therefore, it is necessary to study effect of thermal stratification, especially the unstable condition which takes up most time periods in daytime and a considerably large time periods during night [2].

Three-dimensional building arrays have been widely employed as models in studies of district-scale air flow and pollutant dispersion process. Compared with numerous experimental and numerical works on neutral flow over building arrays [3–7], only a few studies considered unstable thermal stratification [8] and the effect of buoyancy force is far from fully understood. Uehara et al. [9] conducted wind tunnel experiment of thermally stratified flow over squared array with whole ground heated uniformly. In Richards et al.’s wind tunnel study [10], one cube was mounted on the bottom wall and its leeward was heated uniformly. Kanda and Moriizumi [11] measured heat flux from building arrays placed on open ground, while it is difficult to get a full view of flow field from out-door measurement.

Numerical simulation provides a flexible and low-cost method for studies of urban atmospheric environment. So far, most of the numerical studies of thermal effect were conducted using Reynolds averaged Navier Stokes (RANS) simulation [12–16]. Kim and Baik’s RANS [17] simulation of two-dimensional street canyon with bottom heating shows counter-rotating vortex inside canyon, while no such phenomenon is observed in wind tunnel experiment [18] and field measurement [19]. With the development of computers, large eddy simulation (LES) becomes popular as it is able to capture main turbulent structures with acceptable computing time. Li et al. [20] and Cheng and Liu [21] studied the effect of thermal stratification on pollutant dispersion in two-dimensional street canyon with ground heating using LES. Boppana et al. [22] simulated Richards et al.’s experiment [10] and the simulated temperature field show fine agreement with the measured Cai [23] compared canyon flow under windward heating and leeward heating in his LES study of Kovar-Panskus et al.’s experiment [18]. Most of the previous numerical studies, both RANS and LES, simulated two-dimensional canyon and three dimensional configuration of building array was seldom considered. Meanwhile, domain heights in most of these studies are too low to study effect of buoyancy force on flow above canopy.

In this paper, a numerical model able to simulate thick boundary layer flow over three dimensional building arrays is established, which can be used in the study of the effect of buoyancy force on flow both within and above canopy layer. Large eddy simulation with dynamic subgrid scale model is employed. The mean flow and turbulence statistics presented in Uehara et al.’s experiment [9] are used to verify our numerical model and numerical method. Effects of thermal instability on mean velocity and turbulent intensities are further studied by comparing with neutral case in spatial distribution.

2. Numerical Methods

2.1. Numerical Model and Boundary Condition

The aim of the present study is to simulate unstably stratified turbulent boundary layer flow over three-dimensional building array. Some researchers studying neutral turbulent flow over building arrays employ numerical models with shear-free and impermeable condition on the upper boundary and periodic condition in streamwise and spanwise directions [5, 24, 25]. Coceal et al.’s research [5] indicated that impermeable condition only influences flow close to upper boundary. However, vertical convection is enhanced and influence of upper boundary condition will be stronger when buoyancy force is exerted on flow. A numerical model with convective condition 𝜕𝜙/𝜕𝑡+𝑈_𝑛𝜕𝜙/𝜕𝑛=0 on the upper boundary is employed and the flow develops in streamwise direction as shown in Figure 1. Convective boundary condition is also used at the outlet. In spanwise direction, periodic condition is used as buildings are regularly distributed and the flow is equilibrium in this direction. Turbulent boundary layer develops slowly and an extremely long domain is required to obtain a fully developed turbulent flow from laminar flow. To save computing time, simulation is conducted previously in a domain with periodic boundary condition in both streamwise and spanwise directions and building layout in this periodic domain is the same as in Figure 1. Fully developed turbulence obtained in the periodic domain is used as an inlet flow. The ground is heated with a fixed temperature Θ_𝑔 higher than air temperature thus, heat flux from ground induces thermal instability. The walls of buildings are set to be adiabatic, which is proper as model material is styrofoam in the experiment.

Figure 1 

Three-dimensional schematic of numerical model.

Configuration of the array considered in the present paper is the same as Uehara et al.’s scaled model [9]. The array is squared and composed of buildings with uniform height ℎ=0.1 m. Distances between buildings are equal to ℎ in streamwise direction and 0.5ℎ in spanwise; thus, the flow belongs to skimming flow [26]. 10 rows and 6 columns of buildings are mounted on the ground as in Figure 1. The spanwise domain size 𝐿_𝑦 equals 0.9 m. To limit the influence of outlet boundary, there is a distance of 10h between last row and outlet, thus streamwise domain size 𝐿_𝑥 equals 3 m. Uehara et al. measured at a position with boundary depth equals 7 h, thus vertical domain size 𝐿_𝑧 is set to be 0.7 m.

Both neutral and thermally unstable cases are simulated to study the effect of buoyancy force. The free stream velocity 𝑈_𝑓 equals 1.5 m/s. In the unstable case, ground temperature Θ_𝑔=332K and free stream temperature Θ_𝑓=292K, which are the same as the ones of the unstable cases in the experiment. Reynolds number defined with free stream velocity and building height (Re=𝑈_𝑓ℎ/𝑣) equals about 10⁴, while that for real urban area can be up to 10⁶~10⁷. It is difficult for laboratory experiment to reach such a high Reynolds number. In present paper Uehara et al.’s experiment is simulated, and investigation of real urban scale arrays is left for further studies.

Grid distribution is illustrated in Figure 2 on horizontal and vertical planes. Number of grid points on each building is 16 in all three directions. Xie and Castro [24] compared the results with different resolutions and indicated that 16 grids are sufficient to simulate detailed flow field and turbulent statistics. Finer grids are used near the walls and ground, where high shear layers or flow separation appear. Above the building cluster, the grids are stretched with height to save computing time. The total number of grid points are 384 in 𝑥-direction, 144 in 𝑦-direction, and 80 in 𝑧-direction.

(a) Horizontal plane

(b) Vertical plane

Figure 2 

Grid distribution in horizontal and vertical planes.

2.2. Governing Equations

For the filtered velocity 𝑢_𝑖 and filtered pressure 𝑝, continuum equation and momentum equations of incompressible fluid with Boussinesq hypothesis read 𝜕𝑢_𝑖𝜕𝑥_𝑖𝜕𝑢=0,_𝑖+𝜕𝑢𝜕𝑡^𝑖𝑢_𝑗𝜕𝑥_𝑗1=−𝜌𝜕𝑝𝜕𝑥_𝑖+𝜕𝜕𝑥_𝑗𝜏_𝑖𝑗𝜕𝑢+𝑣^𝑖𝜕𝑥_𝑗+𝛼𝑔𝜃𝛿_𝑖3, (1)(2) and transport equation of filter temperature 𝜃 read 𝜕𝜃+𝜕𝜃𝑢𝜕𝑡𝜕𝑥^𝑖_𝑖=𝜕𝜕𝑥_𝑖𝜋_𝑖+𝑣𝜕𝜃Pr𝜕𝑥_𝑖. (3) Here, 𝜏_𝑖𝑗=𝑢_𝑖𝑢_𝑗−𝑢_𝑖𝑢_𝑗 is subgrid scale stress and 𝜋_𝑖=𝜃𝑢_𝑖−𝜃𝑢_𝑖 is subgrid scale thermal flux. 𝛼=1/Θ_ref is thermal expansion coefficient, and Θ_ref is the reference temperature. 𝑥,𝑦,and𝑧 are streamwise, spanwise, and vertical directions respectively, and 𝑢,𝑣, and  𝑤 are velocity components in the three directions. Air density 𝜌 equals 1.208kgm⁻³, kinetic viscosity coefficient 𝑣 equals 1.5×10⁻⁵m²s⁻¹ and molecular Prandtl number Pr is 0.72.

2.3. Subgrid Scale Model

The subgrid scale stress and subgrid scale thermal flux are assumed to be in eddy viscosity and eddy diffusion form, which can be written as: 𝜏_𝑖𝑗=2𝑣_𝑡𝑆_𝑖𝑗+𝜏_𝑘𝑘𝛿_𝑖𝑗3,𝜋_𝑖=𝑘_𝑡𝜕𝜃𝜕𝑥_𝑖, (4) and 𝑆_𝑖𝑗=(𝜕𝑢_𝑖/𝜕𝑥_𝑗+𝜕𝑢_𝑗/𝜕𝑥_𝑖)/2 is the filtered strain tensor. The flow is heterogeneous in horizontal plane and eddy viscosity coefficient 𝑣_𝑡=𝐶²_𝑠Δ²|𝑆| is determined dynamically with procedure proposed by Meneveau et al. [27]. The dynamic procedure is based on the scale similarity between two filter scale Δ₁ and Δ₂. Δ₁ equals Δ and Δ₂(>Δ₁) is set to equal 2Δ in the present paper. Unlike Germano et al.’s procedure [28], which averages on homogeneous planes, Meneveau et al.’s method averages along trajectory of fluid particles and results in two time-integral values 𝐼_𝐿𝑀 and 𝐼_𝑀𝑀, which can be determined by the following equations: 𝜕𝐼_𝐿𝑀+𝑢𝜕𝑡^𝑗𝜕𝐼_𝐿𝑀𝜕𝑥_𝑗=1𝑇𝐿_𝑖𝑗𝑀_𝑖𝑗−𝐼_𝐿𝑀,𝜕𝐼_𝑀𝑀+𝑢𝜕𝑡^𝑗𝜕𝐼_𝑀𝑀𝜕𝑥_𝑗=1𝑇𝑀_𝑖𝑗𝑀_𝑖𝑗−𝐼_𝑀𝑀. (5) The time scale 𝑇 is set to equal 1.5Δ(𝐼_𝐿𝑀𝐼_𝑀𝑀)^−1/8 as proposed by Meveneau et al. [27], 𝐿_𝑖𝑗 equals 𝑢_𝑖𝑢_𝑗−𝑢_𝑖𝑢_𝑗, and 𝑀_𝑖𝑗 equals 2{Δ²₂|𝑆|𝑆_𝑖𝑗−Δ²₁|𝑆|𝑆_𝑖𝑗}. The superscripts “⋅̂” and “⋅̃” represent filtering with Δ₁ and Δ₂, respectively. Equation (5) can be easily discretized and coupled with governing equations, thus model coefficient 𝐶²_𝑠=𝐼_𝐿𝑀/𝐼_𝑀𝑀 can be obtained. Meneveau et al.’s procedure is suitable for heterogeneous flow and previous investigations [29] indicate that dynamic model performs better than standard Smagorinsky model. See more about the details of the model in [27]. As for subgrid scale thermal flux, eddy diffusivity coefficient 𝑘_𝑡(=𝑣_𝑡/Pr_𝑡) is determined with subgrid turbulence Prantel number Pr_𝑡, which is set to equal 0.72.

2.4. Wall Function

To take into account large pressure gradient induced by separation flow behind bluff bodies, wall function proposed by Wang and Moin [30] is used, which solves boundary layer equations on embedded mesh between wall and first near-wall grid.

Consider 𝜕𝜕𝑥_𝑛𝑣+𝑣_𝑡𝜕𝑢_𝑖𝜕𝑥_𝑛=𝐹_𝑖, (6) where 𝐹_𝑖=1𝜌𝜕𝑝𝜕𝑥_𝑖+𝜕𝑢_𝑖+𝜕𝜕𝑡𝜕𝑥_𝑗𝑢_𝑖𝑢_𝑗, (7)𝑛 denotes wall-normal direction and 𝑖,𝑗 belong to the other two directions. The eddy viscosity 𝑣_𝑡 is computed using simple mixing-length eddy viscosity model with near-wall damping in Wang and Moin [30] as follows: 𝑣_𝑡𝑣=𝜅𝑧⁺_𝑤1−𝑒^{−𝑧+𝑤/𝐴2}, (8) where 𝑧⁺_𝑤 is the distance to the wall in wall units and 𝜅=0.4 and 𝐴=19 are the model coefficients. If unsteady and convective terms are ignored in (7) and only pressure gradient, which is determined by outer flow and equals pressure gradient at the first near wall grid, is considered, an explicit expression for local shear stress 𝜏_𝑤𝑖 is obtained.

Consider 𝜏_𝑤𝑖𝜕𝑢=𝜇^𝑖𝜕𝑥_{𝑛𝑥𝑛=0}=𝜌∫^𝛿₀𝑑𝑥_𝑛/𝑣+𝑣_𝑡𝑢_𝑖,𝛿−𝐹_𝑖^𝛿₀𝑥_𝑛𝑑𝑥_𝑛𝑣+𝑣_𝑡, (9) where 𝛿 is the distance of the first near-wall grid point to the wall. In the calculation, 𝑣_𝑡 is obtained with (8) and 𝜏_𝑤𝑖 can be computed directly by one-dimensional integration along normal direction in (9). Such wall model has also been adopted by Boppana et al. [22] in their large eddy simulation of turbulence over single cube with heat flux from surface, and near wall temperature distribution is improved compared with standard wall function.

2.5. Numerical Scheme

Finite volume method (FVM) [31] is employed to discretize governing equations (1)–(3). The spatial numerical scheme is quick and third-order Runge-Kutta scheme is used for time advancement. The discretized equations are solved using simple method on nonstagger grids with momentum interpolation. A total of 1000𝑇_∗ is simulated after the flow approaching statistically steady state, where 𝑇_∗=ℎ/𝑢_∗ and friction velocity 𝑢_∗=(𝜏₀/𝜌)^1/2 is determined by maximum total shear stress 𝜏₀. The following results are time-averaged.

3. Results and Discussion

Cheng and Liu [21] simulated similar array with LES. Height and streamwise distance of Cheng and Liu’s array is the same as those in Uehara et al.’s experiment while spanwise width is infinite, thus the array is two dimensional. Some other authors [12, 13, 20] have also simulated such two dimensional arrays with RANS or LES, while their domain heights are insufficient to obtain a fully developed boundary layer above canopy. Present LES results are compared with Uehara et al.’s experiment, Cheng and Liu’s LES (CLES) in the following.

3.1. Mean Velocity

Uehara et al. measured mean flow at the center of the canyon where depth of boundary layer is about 7 times of the building height. In Figure 3, vertical profiles of time averaged streamwise velocity 𝑈 obtained at the center of canyon before the 8th row of buildings are compared with Uehara et al.’s measurement and CLES for both neutral and unstable cases. The mean velocity is normalized with free stream velocity 𝑈_𝑓. In Figure 3(b), profiles obtained at the centers of canyons before 6th and 7th rows are also illustrated. Difference among the three streamwise positions is negligible and results presented in the following are all obtained at the 8th row. Profiles within and immediately above the canopy layer are also compared between neutral and unstable conditions directly in Figure 3(b).

(a) Neutral

(b) Unstable

Figure 3 

Vertical profiles of mean streamwise velocity normalized with free stream velocity 𝑈𝑓 for both neutral and unstable cases.

For both neutral and unstable cases, our LES results are in fine agreement with the experiment and slightly better than CLES both within and above canopy layer. Flow speed is reduced within the canopy and strong shear layer is observed at roof level for both thermal conditions. Difference between two conditions can be observed. Buoyancy force enhances vertical convection and high momentum fluid in the out flow is transferred to roof level, which results in stronger roof level shear layer in unstable case as seen in Figure 3(b). Consequently, reverse flow in the low canopy is stronger in unstable case.

Spatial distribution of mean flow structure can be better viewed in Figure 4, which illustrates normalized velocity vector (𝑈/𝑈_𝑓,𝑊/𝑈_𝑓) on vertical planes through the middle of the canyon. The flow is from left to right and recirculation is observed inside canopy. Under unstable condition, the recirculation is slightly enhanced and the center moves to the upstream. The center is 0.3 h from windward and 0.7 h above the ground under neutral condition and moves left by 0.08 h and up by 0.06 h under unstable condition. However, no large counter-rotating structure, as presented in RANS results [12, 14], is observed near windward in Figure 4(b). Multiple vortex structures were not observed in Uehara et al.’s experiment (see Figure  12(c) in reference [9]), field measurements [15, 19] and Li et al.’s LES [20] of street canyon with ground heating. Louka et al. [15] simulated their experiment with RANS and concluded that their RANS model overestimated thermal effect on mean flow structure inside canyon.

(a) Neutral

(b) Unstable

Figure 4 

Normalized velocity vector (𝑈/𝑈𝑓,𝑊/𝑈𝑓) on vertical planes through the middle of the canyon.

3.2. Turbulent Intensities

Vertical profiles of root mean square of velocity fluctuation 𝑢_rms=(𝑢^′′𝑢^′′)^1/2 normalized with free stream velocity are compared for both neutral and unstable cases in Figure 5. Comparison with Uehara et al.’s experiment indicates that our numerical simulation is able to capture streamwise turbulent fluctuation both within and above the canopy layer, for both neutral and unstable conditions. Results of CLES are considerably smaller under both thermal conditions. This discrepancy may results from different morphologies between experiment and CLES or different numerical models between CLES and our LES.

(a) Neutral

(b) Unstable

Figure 5 

Vertical profiles of streamwise turbulence intensity 𝑢rms normalized with free stream velocity 𝑈𝑓 for both neutral and unstable case.

Large gradient of 𝑢_rms is observed at roof level, which is consistent with strong shear layer of mean velocity in Figure 3. The peak value is obtained just above the canopy, which equals about 0.18𝑈_𝑓 in neutral case and increases to about 0.22𝑈_𝑓 in unstable case. Above the canopy layer, 𝑢_rms decreases with height monotonously for both thermal conditions, while within the canopy, 𝑢_rms shows slightly uniform distribution vertically. The value is about 0.08𝑈_𝑓 for neutral case and 0.1𝑈_𝑓 for unstable one.

Figure 6 illustrates contours of 𝑢_rms on vertical planes through the middle of the buildings. Intensive turbulent fluctuation is observed above canyon at roof-level, as strong shear layer between out flow and canyon flow contributes to the production of turbulent kinetic energy. 𝑢_rms is lower inside canyon especially near the wall and changes little vertically throughout the canyon. Comparing Figures 6(a) and 6(b), increase of 𝑢_rms under unstable condition is obvious. Besides, large 𝑢_rms region at roof level extends into canyon due to stronger recirculation observed in Figure 4.

(a) Neutral

(b) Unstable

Figure 6 

Contours of streamwise turbulence intensity 𝑢rms normalized with free stream velocity 𝑈𝑓 on vertical plane through the middle of the canyon.

Similarly, rms of vertical velocity fluctuation 𝑤_rms=(𝑤^′′𝑤^′′)^1/2 is compared in Figure 7. The simulated results satisfy well with the experiment, except inside canopy layer where 𝑤_rms is smaller than measured value under unstable condition, and vertical turbulent intensity is much less than measurements in CLES again. Unlike streamwise component, no large gradient of 𝑤_rms is observed at roof level, while like 𝑢_rms, 𝑤_rms is enhanced by buoyancy force throughout the whole domain. Figure 8 illustrates contours of 𝑤_rms on the same vertical planes as in Figure 6. Due to recirculation inside canyon, 𝑤_rms near windward is quite different from that near leeward. Peak value of 𝑤_rms is located in the vicinity of windward, where the downdraft flow penetrates into canyon. Vertical shear layer between downdraft flow and windward contributes to the production of 𝑤_rms. 𝑤_rms is smaller near roof and ground due to restriction of impermeable wall boundary condition.

(a) Neutral

(b) Unstable

Figure 7 

Vertical profiles of vertical turbulence intensity 𝑤rms normalized with free stream velocity 𝑈𝑓 for both neutral and unstable case.

(a) Neutral

(b) Unstable

Figure 8 

Contours of vertical turbulence intensity 𝑤rms normalized with free stream velocity 𝑈𝑓 on vertical plane through the middle of the canyon.

Turbulent kinetic energy TKE=(𝑢²_rms+𝑣²_rms+𝑤²_rms)/2 normalized with free stream velocity 𝑈_𝑓 is illustrated on horizontal planes, which are at a height of 0.5 h. Results are compared between neutral and unstable condition again and variation of turbulence in spanwise direction can be observed, which reveals difference between flow over two-dimensional array and three-dimensional array. High turbulent kinetic energy region is observed in the vicinity of windward where intense vertical turbulence intensity 𝑤_rms exists. Meanwhile, turbulence in the crossroad is strong as the flow in this region is not directly restricted by buildings. TKE is larger throughout the plane for unstable condition reasonably.

3.3. Temperature and Heat Flux

Distributions of temperature and heat flux are discussed in this subsection. Since CLES cannot provide mean temperature and temperature fluctuation, present LES results and measured results are compared in Figure 9. Vertical profile of time averaged temperature is illustrated in Figure 10(a) in the form of (Θ−Θ_𝑔)/(Θ_𝑓−Θ_𝑔), which is the temperature difference between ground and certain height normalized with bulk temperature difference Θ_𝑓−Θ_𝑔. Note that the ground is hottest and larger value of (Θ−Θ_𝑔)/(Θ_𝑓−Θ_𝑔) indicates lower temperature. The simulated mean temperature satisfies well with has the temperature measured result. About 80% of bulk temperature difference is inside canopy layer, and there is a thin layer with large temperature gradient near ground, which takes up most proportion of canopy layer temperature difference. As for the fluctuation of temperature 𝜃_rms illustrated in Figure 10(b), our LES performs generally well, while it is overestimated within canopy layer in comparison with the experiment. Peak of 𝜃_rms is observed near ground which is consistent with the location of large temperature gradient.

(a) Neutral

(b) Unstable

Figure 9 

Contours of turbulence kinetic energy normalized with free stream velocity 𝑈𝑓 on horizontal plane through the middle of the canyon.

(a) Mean temperature

(b) Temperature fluctuation

Figure 10 

Vertical profiles of mean temperature and turbulent intensity for temperature 𝜃rms normalized with total temperature difference (Θ𝑓−Θ𝑔).

Contours of normalized mean temperature (Θ−Θ_𝑓)/(Θ_𝑔−Θ_𝑓) and normalized vertical turbulence heat flux 𝑤^”𝜃^” on the same vertical planes are illustrated in Figure 11. Temperature is higher in lower canopy and decreases with height reasonably. Horizontal variation of temperature is obvious. The air is hotter near leeward and cooler near windward, as low temperature out flow is transported into canyon by recirculation near the windward and heat transfer is less efficient at the corner of leeward.

(a) Temperature

(b) Heat flux

Figure 11 

Contours of normalized mean temperature (Θ−Θ𝑓)/(Θ𝑔−Θ𝑓) and normalized turbulent heat flux 𝑤′′𝜃′′ on vertical plane through the middle of canyon.

Turbulent heat flux, which is in the form of second-order moment of velocity and temperature fluctuation 𝑢_𝑖^′′𝜃^′′, takes up most proportion of heat flux inside canopy layer. When the ground is heated, there is vertical heat flux from the canopy to the out flow. Contour of vertical turbulent heat flux 𝑤^′′𝜃^′′ is illustrated in Figure 11(b). Although temperature is constant on the ground, 𝑤^′′𝜃^′′ is quite nonuniform in the canyon. Positive 𝑤^′′𝜃^′′ means that heat is transferred upward by turbulence, while negative one means that heat is transferred in opposite direction. There is positive 𝑤^′′𝜃^′′ at roof-level above the canyon where heat is transferred out of the canyon. In the vicinity of windward, heat is effectively transferred upward by active vertical fluctuation. In the vicinity of leeward, 𝑤^′′𝜃^′′ is in negative value or small positive value. Heat exchange of leeward is much less effective than that of windward, which results in higher temperature near leeward in Figure 11(a).

4. Conclusion

In this paper, flow field both within and above three dimensional array is investigated with large eddy simulation. Thermal instability is induced by ground heating and buoyancy force is taken into account using Boussinesq hypothesis. Subgrid scale model and wall function employed are suitable for heterogeneous flow with separation behind buildings. Both neutral and unstable cases are simulated to explore the effect of buoyancy force. Overall, the simulated results shows fine agreement with the experiment in both mean flow and turbulent intensities, and perform better than former LES study. Our numerical model is suitable for further studies of more complex flow cases.

It is shown that buoyancy force changes pattern of mean flow inside canopy layer while no secondary reverse circulation is observed, which coincides with wind tunnel experiment and field measurement. The effect of buoyancy force is overestimated by RANS simulations and LES is a better choice for unstably stratified turbulent flow. Moreover, buoyancy force contributes to the production of turbulent kinetic energy, and turbulent intensities are enhanced, especially the vertical component.

Thermal dynamics of high Reynolds number flow over urban-scale buildings with heterogeneous temperature distribution induced by solar radiation, shading, or nonuniform building layout, will be investigated in further studies.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

The authors thank the National Natural Science Foundation of China (Grants 11132005, 10925210, and 11002081), MOST-2011BAK07B01-03, and the National Laboratory for Information Science and Technology for the financial support.