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

Authors: Bobin Wang; Guixiang Cui
Journal: Journal of Applied Mathematics
 (2013)

---

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


---

## Body

## 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 ∂ϕ/∂t+Un∂ϕ/∂n=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 Θg 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 h=0.1 m. Distances between buildings are equal to h in streamwise direction and 0.5h 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 Ly 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 Lx equals 3 m. Uehara et al. measured at a position with boundary depth equals 7 h, thus vertical domain size Lz 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 velocityUf equals 1.5 m/s. In the unstable case, ground temperature Θg=332K and free stream temperature Θf=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=Ufh/v) equals about 104, while that for real urban area can be up to 106~107. 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 Figure2 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 x-direction, 144 in y-direction, and 80 in z-direction.Grid distribution in horizontal and vertical planes.
(a)
Horizontal plane
(b)
Vertical plane


### 2.2. Governing Equations

For the filtered velocityu^i and filtered pressure p^, continuum equation and momentum equations of incompressible fluid with Boussinesq hypothesis read
(1)∂u^i∂xi=0,(2)∂u^i∂t+∂u^iu^j∂xj=-1ρ∂p^∂xi+∂∂xj(τij+v∂u^i∂xj)+αgθ^δi3,
and transport equation of filter temperature θ^ read
(3)∂θ^∂t+∂θ^u^i∂xi=∂∂xi(πi+vPr∂θ^∂xi).
Here, τij=u^iu^j-uiuj^ is subgrid scale stress and πi=θ^u^i-θui^ is subgrid scale thermal flux. α=1/Θref is thermal expansion coefficient, and Θref is the reference temperature. x,y,andz are streamwise, spanwise, and vertical directions respectively, and u,v, and  w are velocity components in the three directions. Air density ρ equals 1.208kgm-3, kinetic viscosity coefficient v equals 1.5×10-5m2s-1 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:(4)τij=2vtS^ij+τkkδij3,πi=kt∂θ^∂xi,
and S^ij=(∂u^i/∂xj+∂u^j/∂xi)/2 is the filtered strain tensor. The flow is heterogeneous in horizontal plane and eddy viscosity coefficient vt=Cs2Δ2|S^| is determined dynamically with procedure proposed by Meneveau et al. [27]. The dynamic procedure is based on the scale similarity between two filter scale Δ1 and Δ2. Δ1 equals Δ and Δ2(>Δ1) 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 ILM and IMM, which can be determined by the following equations:
(5)∂ILM∂t+u^j∂ILM∂xj=1T(LijMij-ILM),∂IMM∂t+u^j∂IMM∂xj=1T(MijMij-IMM).
The time scale T is set to equal 1.5Δ(ILMIMM)-1/8 as proposed by Meveneau et al. [27], Lij equals u^~iu^~j-u^iu^j~, and Mij equals 2{Δ22|S~|S~ij-Δ12|S^|S^ij~}. The superscripts “·^” and “·~” represent filtering with Δ1 and Δ2, respectively. Equation (5) can be easily discretized and coupled with governing equations, thus model coefficient Cs2=ILM/IMM 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 kt(=vt/Prt) is determined with subgrid turbulence Prantel number Prt, 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)∂∂xn(v+vt)∂u^i∂xn=F^i,
where
(7)F^i=1ρ∂p^∂xi+∂u^i∂t+∂∂xju^iu^j,n denotes wall-normal direction and i,j belong to the other two directions. The eddy viscosity vt is computed using simple mixing-length eddy viscosity model with near-wall damping in Wang and Moin [30] as follows:
(8)vtv=κzw+(1-e-zw+/A)2,
where zw+ is the distance to the wall in wall units and κ=0.4 and A=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 τwi is obtained.Consider(9)τwi=μ∂u^i∂xn|xn=0=ρ∫0δ(dxn/v+vt){u^i,δ-F^i∫0δxndxnv+vt},
where δ is the distance of the first near-wall grid point to the wall. In the calculation, vt is obtained with (8) and τwi 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 1000T* is simulated after the flow approaching statistically steady state, where T*=h/u* and friction velocity u*=(τ0/ρ)1/2 is determined by maximum total shear stress τ0. The following results are time-averaged.


## 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 ∂ϕ/∂t+Un∂ϕ/∂n=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 Θg 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 h=0.1 m. Distances between buildings are equal to h in streamwise direction and 0.5h 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 Ly 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 Lx equals 3 m. Uehara et al. measured at a position with boundary depth equals 7 h, thus vertical domain size Lz 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 velocityUf equals 1.5 m/s. In the unstable case, ground temperature Θg=332K and free stream temperature Θf=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=Ufh/v) equals about 104, while that for real urban area can be up to 106~107. 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 Figure2 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 x-direction, 144 in y-direction, and 80 in z-direction.Grid distribution in horizontal and vertical planes.
(a)
Horizontal plane
(b)
Vertical plane


## 2.2. Governing Equations

For the filtered velocityu^i and filtered pressure p^, continuum equation and momentum equations of incompressible fluid with Boussinesq hypothesis read
(1)∂u^i∂xi=0,(2)∂u^i∂t+∂u^iu^j∂xj=-1ρ∂p^∂xi+∂∂xj(τij+v∂u^i∂xj)+αgθ^δi3,
and transport equation of filter temperature θ^ read
(3)∂θ^∂t+∂θ^u^i∂xi=∂∂xi(πi+vPr∂θ^∂xi).
Here, τij=u^iu^j-uiuj^ is subgrid scale stress and πi=θ^u^i-θui^ is subgrid scale thermal flux. α=1/Θref is thermal expansion coefficient, and Θref is the reference temperature. x,y,andz are streamwise, spanwise, and vertical directions respectively, and u,v, and  w are velocity components in the three directions. Air density ρ equals 1.208kgm-3, kinetic viscosity coefficient v equals 1.5×10-5m2s-1 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:(4)τij=2vtS^ij+τkkδij3,πi=kt∂θ^∂xi,
and S^ij=(∂u^i/∂xj+∂u^j/∂xi)/2 is the filtered strain tensor. The flow is heterogeneous in horizontal plane and eddy viscosity coefficient vt=Cs2Δ2|S^| is determined dynamically with procedure proposed by Meneveau et al. [27]. The dynamic procedure is based on the scale similarity between two filter scale Δ1 and Δ2. Δ1 equals Δ and Δ2(>Δ1) 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 ILM and IMM, which can be determined by the following equations:
(5)∂ILM∂t+u^j∂ILM∂xj=1T(LijMij-ILM),∂IMM∂t+u^j∂IMM∂xj=1T(MijMij-IMM).
The time scale T is set to equal 1.5Δ(ILMIMM)-1/8 as proposed by Meveneau et al. [27], Lij equals u^~iu^~j-u^iu^j~, and Mij equals 2{Δ22|S~|S~ij-Δ12|S^|S^ij~}. The superscripts “·^” and “·~” represent filtering with Δ1 and Δ2, respectively. Equation (5) can be easily discretized and coupled with governing equations, thus model coefficient Cs2=ILM/IMM 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 kt(=vt/Prt) is determined with subgrid turbulence Prantel number Prt, 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)∂∂xn(v+vt)∂u^i∂xn=F^i,
where
(7)F^i=1ρ∂p^∂xi+∂u^i∂t+∂∂xju^iu^j,n denotes wall-normal direction and i,j belong to the other two directions. The eddy viscosity vt is computed using simple mixing-length eddy viscosity model with near-wall damping in Wang and Moin [30] as follows:
(8)vtv=κzw+(1-e-zw+/A)2,
where zw+ is the distance to the wall in wall units and κ=0.4 and A=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 τwi is obtained.Consider(9)τwi=μ∂u^i∂xn|xn=0=ρ∫0δ(dxn/v+vt){u^i,δ-F^i∫0δxndxnv+vt},
where δ is the distance of the first near-wall grid point to the wall. In the calculation, vt is obtained with (8) and τwi 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 1000T* is simulated after the flow approaching statistically steady state, where T*=h/u* and friction velocity u*=(τ0/ρ)1/2 is determined by maximum total shear stress τ0. 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 Figure3, vertical profiles of time averaged streamwise velocity U 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 Uf. 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).Vertical profiles of mean streamwise velocity normalized with free stream velocityUf for both neutral and unstable cases.
(a)
Neutral
(b)
UnstableFor 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 Figure3(b). Consequently, reverse flow in the low canopy is stronger in unstable case.Spatial distribution of mean flow structure can be better viewed in Figure4, which illustrates normalized velocity vector (U/Uf,W/Uf) 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.Normalized velocity vector(U/Uf,W/Uf) on vertical planes through the middle of the canyon.
(a)
Neutral
(b)
Unstable


### 3.2. Turbulent Intensities

Vertical profiles of root mean square of velocity fluctuationurms=(u′′u′′¯)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.Vertical profiles of streamwise turbulence intensityurms normalized with free stream velocity Uf for both neutral and unstable case.
(a)
Neutral
(b)
UnstableLarge gradient ofurms 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.18Uf in neutral case and increases to about 0.22Uf in unstable case. Above the canopy layer, urms decreases with height monotonously for both thermal conditions, while within the canopy, urms shows slightly uniform distribution vertically. The value is about 0.08Uf for neutral case and 0.1Uf for unstable one.Figure6 illustrates contours of urms 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. urms is lower inside canyon especially near the wall and changes little vertically throughout the canyon. Comparing Figures 6(a) and 6(b), increase of urms under unstable condition is obvious. Besides, large urms region at roof level extends into canyon due to stronger recirculation observed in Figure 4.Contours of streamwise turbulence intensityurms normalized with free stream velocity Uf on vertical plane through the middle of the canyon.
(a)
Neutral
(b)
UnstableSimilarly, rms of vertical velocity fluctuationwrms=(w′′w′′¯)1/2 is compared in Figure 7. The simulated results satisfy well with the experiment, except inside canopy layer where wrms 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 wrms is observed at roof level, while like urms, wrms is enhanced by buoyancy force throughout the whole domain. Figure 8 illustrates contours of wrms on the same vertical planes as in Figure 6. Due to recirculation inside canyon, wrms near windward is quite different from that near leeward. Peak value of wrms 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 wrms. wrms is smaller near roof and ground due to restriction of impermeable wall boundary condition.Vertical profiles of vertical turbulence intensitywrms normalized with free stream velocity Uf for both neutral and unstable case.
(a)
Neutral
(b)
UnstableContours of vertical turbulence intensitywrms normalized with free stream velocity Uf on vertical plane through the middle of the canyon.
(a)
Neutral
(b)
UnstableTurbulent kinetic energyTKE=(urms2+vrms2+wrms2)/2 normalized with free stream velocity Uf 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 wrms 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 Figure9. Vertical profile of time averaged temperature is illustrated in Figure 10(a) in the form of (Θ-Θg)/(Θf-Θg), which is the temperature difference between ground and certain height normalized with bulk temperature difference Θf-Θg. Note that the ground is hottest and larger value of (Θ-Θg)/(Θf-Θg) 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.Contours of turbulence kinetic energy normalized with free stream velocityUf on horizontal plane through the middle of the canyon.
(a)
Neutral
(b)
UnstableVertical profiles of mean temperature and turbulent intensity for temperatureθrms normalized with total temperature difference (Θf-Θg).
(a)
Mean temperature
(b)
Temperature fluctuationContours of normalized mean temperature(Θ-Θf)/(Θg-Θf) and normalized vertical turbulence heat flux w"θ"¯ 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.Contours of normalized mean temperature(Θ-Θf)/(Θg-Θf) and normalized turbulent heat flux w′′θ′′¯ on vertical plane through the middle of canyon.
(a)
Temperature
(b)
Heat fluxTurbulent heat flux, which is in the form of second-order moment of velocity and temperature fluctuationui′′θ′′¯, 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 w′′θ′′¯ is illustrated in Figure 11(b). Although temperature is constant on the ground, w′′θ′′¯ is quite nonuniform in the canyon. Positive w′′θ′′¯ means that heat is transferred upward by turbulence, while negative one means that heat is transferred in opposite direction. There is positive w′′θ′′¯ 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, w′′θ′′¯ 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).


## 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 Figure3, vertical profiles of time averaged streamwise velocity U 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 Uf. 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).Vertical profiles of mean streamwise velocity normalized with free stream velocityUf for both neutral and unstable cases.
(a)
Neutral
(b)
UnstableFor 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 Figure3(b). Consequently, reverse flow in the low canopy is stronger in unstable case.Spatial distribution of mean flow structure can be better viewed in Figure4, which illustrates normalized velocity vector (U/Uf,W/Uf) 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.Normalized velocity vector(U/Uf,W/Uf) on vertical planes through the middle of the canyon.
(a)
Neutral
(b)
Unstable


## 3.2. Turbulent Intensities

Vertical profiles of root mean square of velocity fluctuationurms=(u′′u′′¯)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.Vertical profiles of streamwise turbulence intensityurms normalized with free stream velocity Uf for both neutral and unstable case.
(a)
Neutral
(b)
UnstableLarge gradient ofurms 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.18Uf in neutral case and increases to about 0.22Uf in unstable case. Above the canopy layer, urms decreases with height monotonously for both thermal conditions, while within the canopy, urms shows slightly uniform distribution vertically. The value is about 0.08Uf for neutral case and 0.1Uf for unstable one.Figure6 illustrates contours of urms 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. urms is lower inside canyon especially near the wall and changes little vertically throughout the canyon. Comparing Figures 6(a) and 6(b), increase of urms under unstable condition is obvious. Besides, large urms region at roof level extends into canyon due to stronger recirculation observed in Figure 4.Contours of streamwise turbulence intensityurms normalized with free stream velocity Uf on vertical plane through the middle of the canyon.
(a)
Neutral
(b)
UnstableSimilarly, rms of vertical velocity fluctuationwrms=(w′′w′′¯)1/2 is compared in Figure 7. The simulated results satisfy well with the experiment, except inside canopy layer where wrms 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 wrms is observed at roof level, while like urms, wrms is enhanced by buoyancy force throughout the whole domain. Figure 8 illustrates contours of wrms on the same vertical planes as in Figure 6. Due to recirculation inside canyon, wrms near windward is quite different from that near leeward. Peak value of wrms 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 wrms. wrms is smaller near roof and ground due to restriction of impermeable wall boundary condition.Vertical profiles of vertical turbulence intensitywrms normalized with free stream velocity Uf for both neutral and unstable case.
(a)
Neutral
(b)
UnstableContours of vertical turbulence intensitywrms normalized with free stream velocity Uf on vertical plane through the middle of the canyon.
(a)
Neutral
(b)
UnstableTurbulent kinetic energyTKE=(urms2+vrms2+wrms2)/2 normalized with free stream velocity Uf 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 wrms 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 Figure9. Vertical profile of time averaged temperature is illustrated in Figure 10(a) in the form of (Θ-Θg)/(Θf-Θg), which is the temperature difference between ground and certain height normalized with bulk temperature difference Θf-Θg. Note that the ground is hottest and larger value of (Θ-Θg)/(Θf-Θg) 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.Contours of turbulence kinetic energy normalized with free stream velocityUf on horizontal plane through the middle of the canyon.
(a)
Neutral
(b)
UnstableVertical profiles of mean temperature and turbulent intensity for temperatureθrms normalized with total temperature difference (Θf-Θg).
(a)
Mean temperature
(b)
Temperature fluctuationContours of normalized mean temperature(Θ-Θf)/(Θg-Θf) and normalized vertical turbulence heat flux w"θ"¯ 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.Contours of normalized mean temperature(Θ-Θf)/(Θg-Θf) and normalized turbulent heat flux w′′θ′′¯ on vertical plane through the middle of canyon.
(a)
Temperature
(b)
Heat fluxTurbulent heat flux, which is in the form of second-order moment of velocity and temperature fluctuationui′′θ′′¯, 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 w′′θ′′¯ is illustrated in Figure 11(b). Although temperature is constant on the ground, w′′θ′′¯ is quite nonuniform in the canyon. Positive w′′θ′′¯ means that heat is transferred upward by turbulence, while negative one means that heat is transferred in opposite direction. There is positive w′′θ′′¯ 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, w′′θ′′¯ 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.



---

*Source: 783187-2013-12-05.xml*