#### 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 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.

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 m. Distances between buildings are equal to in streamwise direction and 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 and free stream temperature , which are the same as the ones of the unstable cases in the experiment. Reynolds number defined with free stream velocity and building height () equals about , while that for real urban area can be up to 10^{6}~10^{7}. 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**

##### 2.2. Governing Equations

For the filtered velocity and filtered pressure , continuum equation and momentum equations of incompressible fluid with Boussinesq hypothesis read and transport equation of filter temperature read Here, is subgrid scale stress and is subgrid scale thermal flux. is thermal expansion coefficient, and is the reference temperature. are streamwise, spanwise, and vertical directions respectively, and , and are velocity components in the three directions. Air density equals , kinetic viscosity coefficient equals and molecular Prandtl number is .

##### 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: and 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 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: The time scale is set to equal as proposed by Meveneau et al. [27], equals , and equals . 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 is determined with subgrid turbulence Prantel number , 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 where 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: where is the distance to the wall in wall units and and 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 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 is simulated after the flow approaching statistically steady state, where and friction velocity 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**

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**

##### 3.2. Turbulent Intensities

Vertical profiles of root mean square of velocity fluctuation 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**

Large gradient of 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 in neutral case and increases to about in unstable case. Above the canopy layer, decreases with height monotonously for both thermal conditions, while within the canopy, shows slightly uniform distribution vertically. The value is about for neutral case and for unstable one.

Figure 6 illustrates contours of 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. is lower inside canyon especially near the wall and changes little vertically throughout the canyon. Comparing Figures 6(a) and 6(b), increase of under unstable condition is obvious. Besides, large region at roof level extends into canyon due to stronger recirculation observed in Figure 4.

**(a) Neutral**

**(b) Unstable**

Similarly, rms of vertical velocity fluctuation is compared in Figure 7. The simulated results satisfy well with the experiment, except inside canopy layer where 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 is observed at roof level, while like , is enhanced by buoyancy force throughout the whole domain. Figure 8 illustrates contours of on the same vertical planes as in Figure 6. Due to recirculation inside canyon, near windward is quite different from that near leeward. Peak value of 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 . is smaller near roof and ground due to restriction of impermeable wall boundary condition.

**(a) Neutral**

**(b) Unstable**

**(a) Neutral**

**(b) Unstable**

Turbulent kinetic energy 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 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 illustrated in Figure 10(b), our LES performs generally well, while it is overestimated within canopy layer in comparison with the experiment. Peak of is observed near ground which is consistent with the location of large temperature gradient.

**(a) Neutral**

**(b) Unstable**

**(a) Mean temperature**

**(b) Temperature fluctuation**

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**

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.