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

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.


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][4][5][6][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][13][14][15][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 twodimensional 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 twodimensional 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.

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.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.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   equals 0.9 m.To limit the influence of outlet boundary, there is a distance of 10ℎ between last row and outlet, thus streamwise domain size   equals 3 m.Uehara et al. measured at a position with boundary depth equals 7ℎ, 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 Θ  = 332 K and free stream temperature Θ  = 292 K, 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 =   ℎ/V) equals about 10 4 , 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.

Governing Equations.
For the filtered velocity û and filtered pressure p, 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. = 1/Θ ref is thermal expansion coefficient, and Θ ref is the reference temperature., , and  are streamwise, spanwise, and vertical directions respectively, and , V, and  are velocity components in the three directions.Air density  equals 1.208 kgm −3 , kinetic viscosity coefficient V equals 1.5 × 10 −5 m 2 s −1 and molecular Prandtl number Pr is 0.72.

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 Ŝ = (û  /  + û  /  )/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 Δ 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   and   , which can be determined by the following equations: The time scale  is set to equal 1.5Δ(    ) −1/8 as proposed by Meveneau et al. [27],   equals ũ ũ − ũ û , and   equals |Ŝ | Ŝ }.The superscripts "⋅" and "⋅" represent filtering with Δ 1 and Δ 2 , respectively.Equation ( 5) can be easily discretized and coupled with governing equations, thus model coefficient  2  =   /  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   (= V  /Pr  ) is determined with subgrid turbulence Prantel number Pr  , which is set to equal 0.72.

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 V  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  = 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 where  is the distance of the first near-wall grid point to the wall.In the calculation, V  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  * = ( 0 /) 1/2 is determined by maximum total shear stress  0 .The following results are time-averaged.

Results and Discussion
Cheng and Liu [21] 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).
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.

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.
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.
Similarly, rms of vertical velocity fluctuation  rms = (    ) 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.
Turbulent kinetic energy TKE = ( 2 rms + V 2 rms +  2 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 twodimensional 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.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.
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.
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).

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.

Figure 2 :
Figure 2: Grid distribution in horizontal and vertical planes.

Figure 3 :
Figure 3: Vertical profiles of mean streamwise velocity normalized with free stream velocity   for both neutral and unstable cases.

Figure 4 :
Figure 4: Normalized velocity vector (/  , /  ) on vertical planes through the middle of the canyon.

Figure 5 :
Figure 5: Vertical profiles of streamwise turbulence intensity  rms normalized with free stream velocity   for both neutral and unstable case.

Figure 6 :
Figure 6: Contours of streamwise turbulence intensity  rms normalized with free stream velocity   on vertical plane through the middle of the canyon.

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

Figure 7 :
Figure 7: Vertical profiles of vertical turbulence intensity  rms normalized with free stream velocity   for both neutral and unstable case.

Figure 8 :
Figure 8: Contours of vertical turbulence intensity  rms normalized with free stream velocity   on vertical plane through the middle of the canyon.

Figure 9 :
Figure 9: Contours of turbulence kinetic energy normalized with free stream velocity   on horizontal plane through the middle of the canyon.

Figure 10 :
Figure 10: Vertical profiles of mean temperature and turbulent intensity for temperature  rms normalized with total temperature difference (Θ  − Θ  ).
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.
7 times of the building height.In Figure3, 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