Urban-Scale Computational Fluid Dynamics Simulations with Boundary Conditions from Similarity Theory and a Mesoscale Model

Abstract

Mesoscale numerical weather prediction models usually provide information regarding environmental parameters near urban areas at a spatial resolution of the order of thousands or hundreds of meters, at best. If detailed information is required at the building scale, an urban-scale model is necessary. Proper definition of the boundary conditions for the urban-scale simulation is very demanding in terms of its compatibility with environmental conditions and numerical modeling. Here, steady-state computational fluid dynamics (CFD) microscale simulations of the wind and thermal environment are performed over an urban area of Kozani, Greece, using both the k-ε and k-ω SST turbulence models. For the boundary conditions, instead of interpolating vertical profiles from the mesoscale solution, which is obtained with the atmospheric pollution model (TAPM), a novel approach is proposed, relying on previously developed analytic expressions, based on the Monin Obuhkov similarity theory, and one-way coupling with minimal information from mesoscale indices (V_y = 10 m, T_y = 100 m, L_*). The extra computational cost is negligible compared to direct interpolation from mesoscale data, and the methodology provides design phase flexibility, allowing for the representation of discrete urban-scale atmospheric conditions, as defined by the mesoscale indices. The results compared favorably with the common interpolation practice and with the following measurements obtained for the current study: SODAR for vertical profiles of wind speed and a meteorological temperature profiler for temperature. The significance of including the effects of diverse atmospheric conditions is manifested in the microscale simulations, through significant variations (~30%) in the critical building-related design parameters, such as the surface pressure distributions and local wind patterns.

1. Introduction

The urban-scale microclimate is dependent on urban environmental parameters, local building design and construction, and on the parameters related to mesoscale weather. This complex interaction necessitates advanced analysis tools, almost always at different scales, creating the issue of passing information from one scale to another in an efficient, accurate, and compatible way. Research on the subject is active, and there are several approaches that have been proposed in the literature. Coupling a mesoscale model with a computational fluid dynamics (CFD) model, for urban-scale flow simulations, has been a prominent choice for some time, e.g., with the MM5 model [1] or the WRF model [2], and others. The direct application to several engineering fields, such as wind energy assessment, pollutant dispersion, wind comfort studies, and urban flows, has been performed. One of the main advantages of this approach is that mesoscale physics are captured and passed to a CFD simulation, which is capable of including geometric and environmental details at the microscale. Among the many advanced turbulence modeling approaches that are available for CFD simulations, Reynolds-averaged Navier–Stokes (RANS) models allow for an acceptable balance between accuracy and calculation time, and represent a valid alternative to other methods, such as LES, that have higher resource demands [3]. However, regardless of the turbulence model being used, one of the main challenges in urban-scale CFD simulations remains the definition of appropriate values of atmospheric variables at the boundaries of the CFD computational domain [2,4].

A common choice of boundary conditions for the simulation of the atmospheric boundary layer (ABL) is the use of logarithmic or power law [5] wind profiles, often taking into account parameters from the Monin–Obukhov similarity theory [6]. However, the application of these profiles has not yet been site-specific [6], and they generally assume steady-state and horizontally uniform conditions. Furthermore, several issues arise during the application of boundary conditions for CFD simulations of atmospheric flows, not the least of which is the undesirable streamwise gradients in the vertical profiles of mean wind speed and turbulence. One of the first studies addressing the subject of proper boundary conditions for CFD studies of wind engineering problems at the urban scale was that of [7], which was specific to the widespread k-ε turbulence model. The study has been revisited [8] and a number of other turbulence modeling approaches have also been addressed. Although the requirements for the boundary conditions have been worked out, their application, especially within commercial CFD software, still remains a challenge in some cases [9,10], necessitating modification of the turbulence model constants [6]. This creates an issue of whether or not to use the same constants for the land surface and the building surfaces. On the other hand, there are reported studies [11] that indicate that the internally, building-generated turbulence may override deviations of up to 50% in the inflow values.

The issue of processing meteorological data from the mesoscale, for use in local-scale studies, has been an issue of research for many years [12]. Since the variability in atmospheric conditions is suitably considered at the mesoscale, a realistic approach to defining a boundary condition for detailed microscale CFD modeling is to use the mesoscale data [13]. Coupling to the mesoscale data has been applied in combination with advanced turbulence modeling, such as LES for the microscale [14,15], but it is used with RANS modeling as well (most frequently with variants of the k-ε model, e.g., RNG k-ε [1], k-ε [2], comprehensive k-ε [4], Sogachev k-ε [16]). Further effort is being made to bridge the gap between these two turbulence modeling approaches, specifically for meso- to micro-scale coupling [17]. An issue arising when creating boundary conditions for CFD is how to interpolate the data from the mesoscale spatial resolution, of the order of kilometers or hundreds of meters, to the urban scale. Several studies apply a simple linear interpolation [1,2], or a power law for the region close to the ground [16], while others apply more sophisticated methods, such as nudging [18], in order to increase the accuracy. For the issue of the temporal variation in the boundary values, an obvious approach is to perform unsteady microscale simulations with LES [14] or RANS [16] turbulence modelling, but this puts heavy demands on computational resources and will only deal with the time period that corresponds to the mesoscale simulations [18]. In some applications, e.g., for wind energy assessment, time averaging for long time periods, in the order of years, is needed, and so steady-state simulations are usually performed. In these cases, one approach is to use analytic expressions [19] for the interpolation process from the meso- to micro-scale. Validation of this approach [20] showed that the coupling procedure significantly improves the results, especially in terms of horizontal variation. An alternative approach has also been applied [21], whereby characteristic atmospheric conditions are identified from the mesoscale simulations and a sequence of representative steady-state solutions is calculated, thus creating a catalogue of precomputed conditions, corresponding to discrete, but locally relevant, atmospheric conditions.

It is often the case that a parametric study needs to be conducted for a specific area, and this should include a range of atmospheric conditions. The possibilities are infinite and it is extremely demanding to perform mesoscale simulations in order to include all the possible conditions. There is, therefore, a need to define a range of atmospheric conditions as inputs to microscale simulations, ideally detached from the mesoscale simulation, but still representative of realistic local conditions. In the present study, we investigate the use of analytic expressions [12], based on the Monin–Obukhov similarity theory, for prescribing boundary conditions of the main variables on the lateral boundaries of a CFD computational domain for an urban area. The expressions rely on sparse information from the mesoscale; in fact, only three parameter values are used, instead of detailed vertical profiles, thus simplifying the coupling between the two models. Furthermore, the mesoscale information defines locally relevant meteorological states, to which the analytic expressions will correspond and can potentially be used to perform quasi-steady microscale simulations that are defined by, but computationally detached from, these mesoscale states.

The work being presented is part of a coordinated effort involving numerical weather predictions at the mesoscale, CFD modeling at the building scale, and experimental measurements of wind speed and temperature profiles at a number of locations above the SW area of the city of Kozani in Northern Greece (Figure 1). The main objective of the research effort is to develop a numerical method for microscale simulations in urban areas that can be used to provide representative, discrete microclimate conditions at the building scale, from which information regarding the thermal and wind environment around the building may be drawn. A significant challenge arising through this objective is to perform the microscale simulations asynchronously from the mesoscale ones, performed here using the TAPM model [22], while retaining the locally relevant meteorological information in the simplest possible form, in order to facilitate the coupling procedure. Here, we evaluate a novel methodology, implementing analytic expressions, instead of interpolating from mesoscale data, to derive boundary conditions for representative atmospheric conditions. The atmospheric conditions are defined by minimal, pointwise, mesoscale indices. We also evaluate the implementation of the standard form of two RANS turbulence models. We do not modify the turbulence models, in order to evaluate the effects of the explicit inclusion of the buildings’ geometry [11]. Through comparison of the vertical profiles with those from the more common interpolation practices, and with measurement data, our combination of analytic expressions from the Monin–Obukhov theory, with empirical relations for meteorological parameters, proves to be a reliable and versatile approach for defining boundary conditions for microscale simulations. Mesoscale information may be defined via simple indices, and significant effects on the urban-scale flow and temperature fields were observed.

2. Materials and Methods

2.1. Measurements

Measurements were conducted at two locations in the city of Kozani, Greece (Figure 1) from January 2014 to January 2015. Temperature measurements were performed at location (1) (40.298809°, 21.799331°) while wind speed measurements were performed at location (2) (40.301718°, 21.800835°) (Figure 1). For temperature measurements, a meteorological temperature profiler (MTP-5, R.P.O. ATTEX) was used to measure 5 min average atmospheric air temperature at increments of 50 m up to a height of 600 m. For wind speed and direction, an acoustic sounder (BMETEK, Phased Array Doppler SODAR PCS.2000-16/MF) was used, providing 10 min average measurements at increments of 20 m from 40 m up to 300 m above ground. A detailed description of the experimental measurements of wind speed and air temperature may be found in [23,24]. Results of these measurements will be compared with the following CFD and mesoscale simulations.

2.2. Mesoscale Modeling

For the mesoscale numerical weather prediction, the meteorological module of the air pollution model (TAPM) [22] was used. The meteorological component of TAPM is a three-dimensional, nestable, prognostic model solving incompressible, non-hydrostatic, primitive equations with a terrain-following vertical coordinate. The Exner pressure function is split into hydrostatic and non-hydrostatic components, and the non-hydrostatic component necessitates solution of a Poisson equation. At the surface, a vegetative canopy, soil scheme as well as an urban scheme are used. Radiative fluxes are included at the surface and at upper levels. Based on surface similarity and urban parameterizations of surface parameters, such as roughness length and displacement height, the model estimates surface turbulent fluxes. The standard E-ε turbulence closure is used for calculation of eddy diffusivity in the ABL. Non-local buoyancy flux is explicitly considered in TAPM under convective conditions. More details on TAPM equations and parameterization, including the numerical methods used to solve the model equations, can be found in a series of publications related to modeling [25], testing for different seasons [26], or even year-long [27,28] and as an analysis tool [29]. Specifically for the present study, TAPM was configured over a 30 × 30 km region with five (nested) grids of 21 × 21 horizontal points with 30, 10, 3, 1 and 0.3 km grid spacing, respectively. Forty-five vertical levels were defined, starting at 10 m, with 5 levels between 0 and 100 m, and gradually increasing to 8 km. TAPM was used to model the meteorology of 9 August 2014 for the region of Kozani. The topography of the area was imported in high resolution (90 m, STRM3) and in order to improve simulations of urban environments, TAPM was modified by incorporating four urban land surface types (low, medium, and high density and CBD, i.e., a city’s commercial center, usually characterized by tall office and residence buildings) replacing the existing single urban surface [29]. Areas with over 15 dwellings per 10,000 m² were considered high density, medium-density areas corresponded to 10–15 dwellings per 10,000 m², and low-density areas to below 10 dwellings per 10,000 m². The values of the parameters σ_U (urban heat capacity), a_U (urban albedo), A_U (urban anthropogenic heat flux, W m⁻²), k_U (conductivity), and zo_U (urban roughness length) are different for each land use category. Details and validation of the TAPM configuration and runs have been presented in [23] where comparisons were found agreeable between meteorological forecasts and observations for both surface and ABL meteorological parameters, as follows: wind speed and potential temperature vertical profiles, friction velocity (u_*), surface turbulent heat flux (Hs) and mixing height (Z_i) of the ABL. Furthermore, both the index of agreement (IOA) and the root mean square error (RMSE) showed good agreement.

2.3. Urban Microscale Modeling

Typically, urban topography and buildings are not included in mesoscale modeling. Instead, their macroscopic effect is included through a canopy model and/or local surface parameters. In the present case, one-way coupling is assumed between the mesoscale and microscale simulations, i.e., the mesoscale results are taken into account in the CFD simulations, but not vice versa. Land use data for TAPM was available at a spatial resolution of 3 × 3 km while at the innermost grid resolution of 300 × 300 m, the model was modified by incorporating four urban land surface types, distributed according to the available information on population distribution within the city [24]. For extracting information at higher spatial resolution, the CFD model must include details of the urban geometry. In order to achieve this, a plan view of the city was geolocated in a 3D design software and an area of 700 × 750 m was chosen around the sites where experimental measurements were available. It should be noted that this corresponds to a region slightly larger than 2 × 2 of the finest (0.3 km) grid cells of the mesoscale domain, i.e., ~4 locations of mesoscale information throughout the microscale computational domain. Within the chosen computational domain, variation in terrain elevation was below 2% and was neglected, assuming a flat terrain. Buildings’ plan views were traced and extruded to the required height (Figure 1) with a mean building height of 10 m, in agreement with the TAPM land use data [30] that gave land use categories of medium to high urban density for the chosen area [24]. The final computational domain was 750 × 700 × 100 m including a 100 m fetch around the buildings in order to keep the computational domain boundaries away from building effects (Figure 2).

For the CFD calculation, the 3D geometry was processed using an in-house algorithm to define the solid boundaries within a Cartesian grid of 142 × 77 × 152 ≈ 1.6 M cells (Figure 2). Grid resolution was uniform at 5 m in the two horizontal directions. In the vertical direction, resolution was 1 m up to a height of 50 m and then gradually coarser up to the top of the domain. There was no effort made to align building geometries with the grid and so the Cartesian spatial discretization leads to a step-like representation of building faces that are not aligned with the grid. This is not expected to alter the larger scale characteristics being discussed here, but would be important for more sensitive analyses related, e.g., to building ventilation or infiltration [31].

The steady-state volume averaged Navier–Stokes equations are solved in three dimensions on a Cartesian collocated grid implementing Rhie–Chow corrections [32] to avoid checkerboard pressure effect when using the SIMPLE pressure correction algorithm [33]. Second-order upwind discretization of the convection terms was used and discretization of diffusion and other terms was also second order. For turbulence closure, both the k-ε [34] and the k-ω SST model [35] have been used with enhanced wall functions [36]. The standard k-ε model is a robust general-purpose turbulence model, but is well known to introduce increased turbulence production, particularly near stagnation points. The SST k-ω model is a blending between the k-ε and k-ω models and is frequently proposed as a reliable alternative. For buoyancy terms, Boussinesq approximation has been used and turbulence production, due to buoyancy, is implemented using the procedure of [37]. Previous application of this CFD model to external flows past buildings and urban geometries has proven it to be reliable [38,39].

Special focus has been placed on the boundary conditions for the CFD computational domain. At outlet boundaries, zero-gradient conditions are applied, along with a mass flow rate correction, to aid convergence. Throughout the top boundary, in order to preserve mesoscale conditions, the uppermost values of the inlet velocity and turbulence intensity are imposed, as suggested by [12,13]. On the ground, wall functions were applied [36], but no special treatment for shear stress or modification of turbulence model parameters was applied. Although this is an active subject of research, it applies to regions that do not include detailed geometric features of the buildings [13]. In fact, there are studies that indicate that the effects of the buildings in generating turbulence overshadow even large variations in upstream turbulence levels [11]. In the present case, it was left to the momentum equations and the turbulence model to generate the atmospheric boundary layer shear stress within the urban setting. Although the CFD code is capable of dealing with shortwave [38] and longwave radiation [40], at this stage, radiation effects have been neglected in the microscale simulation and the temperature on the lower boundary (ground and buildings) has been considered constant, its difference with the air temperature depending only on atmospheric conditions.

The application of inlet profiles for the CFD studies is one of the major subjects of the present work. For the CFD simulations, velocity, temperature, turbulence kinetic energy and turbulence kinetic energy dissipation rate must be prescribed at the inlet plane of the computational domain. The inlet plane is chosen among the lateral boundaries of the computational domain according to the desired wind direction, which in turn determines the values of the horizontal components of the wind speed. So, for a wind direction from the northwest both the north and west boundary would be considered inlet planes and the velocity components set to give the correct angle of 45° to these planes. The mean value for the vertical component can safely be assumed as zero for the present case of flat terrain. All parameters are assumed to be uniformly distributed across the inlet plane in the horizontal direction and vary only in the vertical direction, i.e., distance from the ground. The following two different methods were examined for prescribing the vertical variation (profiles) in values along the inlet boundaries: (a) results of the TAPM mesoscale calculations were interpolated to the CFD grid and (b) the parameterizations proposed by COST 710 WG3 [12], based on Monin–Obukhov similarity and supplemented here for the calculation of mixing height [41] for unstable [42] or stable [43] and nocturnal [44] atmosphere.

For the case of directly applying mesoscale results as inlet boundary conditions to the CFD simulations, a linear interpolation of the TAPM results for velocity, temperature and turbulence kinetic energy from the mesoscale to the CFD grid was performed. This is a common approach [1,2], but it should be kept in mind that the first computational grid point for the mesoscale calculations is at a height of 10 m from the ground and the CFD grid starts at 1 m. Alternatives of interpolating using a logarithmic or power law [5] were not applied so as to allow for a clear comparison with the COST 710 parameterization, which will briefly be described following this.

The parameterized approach that relies only on sparse point values and meteorological indices from the mesoscale results is particularly interesting in situations when parametric studies are to be performed at a design or simulation stage of any urban computational wind application, e.g., with regard to varying wind directions, stability conditions, etc. These inlet boundary conditions may be used to pre-define case-specific meteorological conditions that can subsequently be used independently or with minimal mesoscale input [21]. Here, the COST 710 parameterizations [12] were used to calculate analytic expressions. For the calculation procedure, point data from the TAPM results are needed, but are limited to a value V_{y = 10} of the wind speed at y = 10 m, a value T_y=100 of the air temperature at y = 100 m and the values of the Obukhov length (L_*) and bulk Richardson number (Ri), as stability indices [12], as follows:

$${\mathrm{L}}_{*}=\frac{{\mathrm{u}}_{{}_{*}}^3{\mathrm{c}}_{\mathrm{p}}\mathsf{\rho }\mathsf{\theta }}{\mathsf{\kappa }\mathrm{g}\mathrm{H}}, \mathrm{R}\mathrm{i}=\frac{\mathrm{g}}{\mathsf{\theta }}\frac{\mathsf{\Delta }\mathsf{\theta }}{\mathsf{\Delta }\mathrm{y}}{\left(\frac{\mathrm{y}}{{\mathrm{u}}_{\mathrm{y}}}\right)}^2$$

(1)

where θ is potential temperature, H is the thermal exchange with the ground, c_p is the specific heat capacity of air under constant pressure and ρ is air density. If the thermal exchange with the ground (H) is non-zero, Monin–Obukhov theory predicts that the Obukhov length scale is a measure of the influence of buoyancy as it expresses the ratio of mechanically generated-to-turbulence generated by buoyancy. On the other hand, the Richardson number expresses the ratio of buoyancy-to-flow shear. Both of these are indicators of atmospheric stability. In the equations following, only the L_* value is necessary with the Ri number used as a confirmative index.

With these three values (V_y=10, T_y=100, L_*) from the mesoscale data, the calculation procedure for the inlet vertical profiles of velocity, turbulence kinetic energy and temperature can proceed. A schematic of the flow of data is presented in Figure 3.

The methodology is summarized below based on [12], who recommend the Monin–Obukhov stability theory (MOST) as essentially the only theory tested widely and applied to describe vertical distributions of mean horizontal wind speed u(y) potential temperature θ(y) and turbulence σ_x, σ_y in atmospheric boundary layer processes (y < 200 m), as follows:

$$\begin{array}{c}\mathrm{u}(\mathrm{y})=\frac{{\mathrm{u}}_{*}}{\mathsf{\kappa }}\left[\ln \frac{\mathrm{y}}{{\mathrm{y}}_{\mathrm{o}}}-{\mathsf{\psi }}_{\mathrm{m}}\left(\frac{\mathrm{y}}{{\mathrm{L}}_{*}}\right)\right]\\ \mathsf{\theta }(\mathrm{y})-\mathsf{\theta }({\mathrm{y}}_{\mathrm{o}})=\frac{{\mathsf{\theta }}_{*}}{\mathsf{\kappa }}\left[\ln \frac{\mathrm{y}}{{\mathrm{y}}_{\mathrm{o}}}-{\mathsf{\psi }}_{\mathrm{h}}\left(\frac{\mathrm{y}}{{\mathrm{L}}_{*}}\right)\right]\end{array}$$

(2)

In (2), ψ_m = 0 and ψ_h = 0 will lead to the logarithmic law of the wall, which is valid for neutral conditions. The MOST theory results in expressions for these functions and essentially modifies turbulence generation and destruction, depending on atmospheric stability. So, for an unstable atmosphere, the following equations are used:

$$\begin{array}{c}{\mathsf{\psi }}_{\mathrm{h}}(\frac{\mathrm{y}}{{\mathrm{L}}_{*}})=2\ln \left[0.5\left(1+\sqrt{1-16\raisebox{1ex}{$\mathrm{y}$}\!\left/ \!\raisebox{-1ex}{${\mathrm{L}}_{*}$}\right.}\right)\right]\\ {\mathsf{\psi }}_{\mathrm{m}}(\frac{\mathrm{y}}{{\mathrm{L}}_{*}})=\ln \left[\frac{1+{{\mathrm{x}}^{\prime }}^2}{2}{\left(\frac{1+{\mathrm{x}}^{\prime }}{2}\right)}^2\right]-2\mathrm{ }\arctan ({\mathrm{x}}^{\prime })+\frac{\mathrm{\pi }}{2}\\ {\mathrm{x}}^{\prime }={\left(1-16\frac{\mathrm{y}}{{\mathrm{L}}_{*}}\right)}^{1/4}\end{array}$$

(3)

and for stable conditions, the following equation is used:

$${\mathsf{\psi }}_{\mathrm{m}}(\frac{\mathrm{y}}{{\mathrm{L}}_{*}})=-5\frac{\mathrm{y}}{{\mathrm{L}}_{*}}, {\mathsf{\psi }}_{\mathrm{h}}(\frac{\mathrm{y}}{{\mathrm{L}}_{*}})=-5\frac{\mathrm{y}}{{\mathrm{L}}_{*}}$$

(4)

where the following applies:

$$\begin{array}{l}{\mathsf{\theta }}_{*}=-\frac{\mathrm{H}}{{\mathrm{c}}_{\mathrm{p}}\mathsf{\rho }{\mathrm{v}}_{*}}\\ \begin{array}{ccc}\mathrm{H}={\mathrm{c}}_{\mathrm{p}}\mathsf{\rho }\overline{{\mathrm{v}}^{\prime }{\mathrm{T}}^{\prime }},& {\mathrm{v}}_{*}={\left(\frac{\mathrm{g}\mathrm{H}{\mathrm{y}}_{\mathrm{i}}}{{\mathrm{c}}_{\mathrm{p}}\mathsf{\rho }{\mathrm{T}}_{\mathrm{o}}}\right)}^{1/3},& {\mathrm{v}}_{*}={\mathrm{u}}_{*}{\left(-\frac{\mathrm{h}}{\mathsf{\kappa }{\mathrm{L}}_{*}}\right)}^{1/3}\end{array}\end{array}$$

(5)

As these functions are dependent on L_*, and L_* is dependent on the friction velocity u_*, a non-linear equation arises. Although empirical relations exist to relate L_* with y_o, the equation may easily be solved numerically as well. For turbulence kinetic energy k = (2σ_x² + σ_y²)/2, again there are different expressions depending on atmospheric stability, as follows:

$$\begin{array}{l}\mathrm{when} -1000<{\mathrm{L}}_{*}< 0 \left(\mathrm{unstable}\right)\\ \hspace{1em}\begin{array}{l}{\mathsf{\sigma }}_{\mathrm{x}}(\mathrm{y})={\mathrm{u}}_{*}{\left[0.35{\left(-\frac{\mathrm{h}}{\mathsf{\kappa }{\mathrm{L}}_{*}}\right)}^{2/3}+\left(2-\frac{\mathrm{y}}{\mathrm{h}}\right)\right]}^{1/2}\\ {\mathsf{\sigma }}_{\mathrm{y}}{}^3(\mathrm{y})={\left[1.6{\mathrm{u}}_{*}^2\left(1-\frac{\mathrm{y}}{\mathrm{h}}\right)\right]}^{3/2}+1.2{\mathrm{w}}_{*}^3\left(\frac{\mathrm{y}}{\mathrm{h}}\right){\left(1-0.9\frac{\mathrm{y}}{\mathrm{h}}\right)}^{3/2}\end{array}\}\\ \mathrm{when} \left|{\mathrm{L}}_{*}\right| > 1000 \left(\mathrm{neutral}\right)\\ \hspace{1em}{\mathsf{\sigma }}_{\mathrm{x}}(\mathrm{y})={\mathsf{\sigma }}_{\mathrm{y}}(\mathrm{y})=1.3{\mathrm{u}}_{*}\exp \left(-2\frac{\mathrm{f}\mathrm{y}}{{\mathrm{u}}_{*}}\right)\\ \mathrm{when} 0 < {\mathrm{L}}_{*}< 1000 \left(\mathrm{stable}\right)\\ \hspace{1em}{\mathsf{\sigma }}_{\mathrm{y}}(\mathrm{y})={\mathsf{\sigma }}_{\mathrm{x}}(\mathrm{y})=1.3{\mathrm{u}}_{*}\left(1-\frac{\mathrm{y}}{\mathrm{h}}\right)\end{array}$$

(6)

where f is the Coriolis parameter (f = 9.33·10⁻⁵ rad/s for Kozani at φ = 40°) and h is the mixing height. The mixing height is essentially the top of the atmospheric boundary layer, i.e., the boundary layer height (y_i) for neutral conditions [41], as follows:

$${\mathrm{y}}_{\mathrm{i}}\approx \mathrm{h}=0.2\frac{{\mathrm{u}}_{*}}{\mathrm{f}}$$

(7)

For unstable conditions, there is usually an inversion capping the mixed region close to the ground and the mixing height is limited to the height of this inversion. The height of the inversion varies throughout the day, beginning close to the ground in early morning, with buoyancy-driven turbulence generation due to solar heating of the ground, and then increasing. An equation for its daily variation may be derived from a thermal energy balance [45], as follows:

$$\mathrm{h}(\mathrm{t})=\sqrt{\raisebox{1ex}{$2{\displaystyle \underset{{\mathrm{t}}_{\mathrm{o}}}{\overset{\mathrm{t}}{\int }}\mathrm{H}\mathrm{d}\mathrm{t}}$}\!\left/ \!\raisebox{-1ex}{${\mathrm{c}}_{\mathrm{p}}\mathsf{\rho }\left({\mathrm{\gamma }}_{\mathrm{d}}-\mathrm{\gamma }\right)$}\right.},$$

(8)

where (γ_d = 9.86 × 10⁻³ Κ/m) is the dry adiabatic lapse rate and (γ) is the lapse rate in the atmosphere (γ = −ϑΤ/ϑz) at sunrise and H is the thermal flux to the ground. Here, when H is unknown the mixing height is approximated by (7). Finally, for stable conditions mixing is solely due to friction with the ground while the stable layer tends to dampen it out and thus the mixing height is dependent on the Obukhov length [43], as follows:

$${\mathrm{h}}_{\mathrm{e}\mathrm{q}}\approx 0.4\sqrt{\frac{{\mathrm{u}}_{*}}{\mathrm{f}}{\mathrm{L}}_{*}}.$$

(9)

The effect of surface roughness is taken into account in the above methodology through the roughness length, where y_o = 0.8 m has been used, based on the land categories from TAPM (medium to high urban density, Ref. [30]).

For both the analytic expressions and the case where profiles are derived directly from TAPM results, inlet values for turbulence kinetic energy dissipation rate (ε)—for the k-ε model—or specific dissipation rate (ω)—for the k-ω SST model—must still be calculated from the available values of turbulence kinetic energy. For the dissipation rates, profiles are defined through [13].

$$\begin{array}{ccc}\mathsf{\epsilon }=\raisebox{1ex}{${\mathrm{u}}_{*}^3$}\!\left/ \!\raisebox{-1ex}{$\mathsf{\kappa }\mathrm{y}$}\right.,& \mathsf{\omega }=\raisebox{1ex}{$\mathsf{\epsilon }$}\!\left/ \!\raisebox{-1ex}{${\mathrm{C}}_{\mathsf{\mu }}\mathrm{k}$}\right.,& {\mathrm{u}}_{*}^{}=\end{array}{\mathrm{C}}_{\mathsf{\mu }}^{1/4}{\mathrm{k}}^{1/2},$$

(10)

where u_* is the friction velocity and C_μ is a model constant of the standard k–ε model (=0.09). Typically, in order for the boundary conditions and the turbulence model to be compatible, the value of the von Karman constant (κ) should satisfy the relation [7].

$${\mathsf{\kappa }}^2=\left({\mathrm{C}}_{\mathsf{\epsilon }2}-{\mathrm{C}}_{\mathsf{\epsilon }1}\right){\mathsf{\sigma }}_{\mathsf{\epsilon }}\sqrt{{\mathrm{C}}_{\mathsf{\mu }}},$$

(11)

which leads to κ = 0.4326 when C_ε2 = 1.92, C_ε1 = 1.44, σ_ε = 1.3, C_μ = 0.09. This deviates from the commonly used value of κ = 0.42, also used here, but this is not expected to lead to any significant errors [8].

3. Results and Discussion

Simulations were performed for 9 August 2014, a day on which measurements of temperature, and wind speed and direction were available, as were TAPM mesoscale calculation results. The date was chosen at random, within the time period that the research project was running, and it included a sufficient range of atmospheric conditions from unstable to neutral to stable. It should be noted, however, that the methodology can be applied to any date or standard meteorological condition that is adequately described by the MOST theory. We chose four different atmospheric conditions, determined from the available data for the specific day, and chosen based on the calculated bulk Richardson number and Obukhov length, to correspond to conditions B, C, D, and E on the classic Pasquil-Gifford categorization scheme ([46] from A, extremely unstable, to G, very stable, with D being the neutral class), i.e., two unstable conditions (B, C), a neutral (D), and a stable (E) condition, respectively. The details of these four conditions are presented in

Table 1. Characterization of four stability conditions over Kozani, Greece at different times on 9 August 2014, based on TAPM mesoscale calculations.

Kozani 09/08/2014	Uy = 10 m (m/s)	L* (m)	Ri	Pasquil
06:00–07:00	3.5 (NW)	999	0.00164	D
07:00–08:00	3.1 (N)	−426.66	−0.010	C
13:00–14:00	2.6 (N)	−44.64	−0.009	B
23:00–00:00	3.5 (N)	270.28	0.00651	E

, and it should be noted that the wind speed at y = 10 m does not vary much (within 2.6–3.5 m/s) for the four cases.

The four conditions in Table 1 provided us with two sets of four inlet profiles for each of the following: velocity, temperature, turbulence kinetic energy, and turbulence kinetic energy dissipation, to use for the CFD calculations. The first set derives from the interpolation of the mesoscale data onto the CFD grid, and the second from the COST 710 methodology presented in the previous section, to create analytic expressions of vertical profiles. It should be noted that the computational cost of deriving the vertical profiles from the COST 710 methodology is of the same order as that of interpolating mesoscale data onto a microscale, CFD, grid, i.e., no iterative procedures are necessary and so there is no advantage in this respect. However, the COST 710 methodology relies only on the three pointwise values to create the profiles, whereas the interpolation presupposes a costly mesoscale simulation. Obviously, a mesoscale simulation provides a higher level of detail, in terms of location-specific atmospheric conditions, but if only representative conditions are desired, a hypothetical set of meteorological conditions can be used to create the profiles with the proposed methodology in minimal time and computational effort.

The results from the mesoscale calculations and CFD calculations at point 2 (x = 387 m, z = 289 m, in Figure 1), using inlet profiles based on the COST 710 methodology (taking into account the mesoscale values of V_y = 10, T_y = 100, L_*), are presented in Figure 4, for the wind speed at the four chosen atmospheric conditions (B, C, D, E in Table 1). One of the advantages of the use of the analytic expressions is implemented here as well, i.e., a parametric study and sensitivity analysis, at each atmospheric condition, using the same inlet profiles for the four main wind directions (north, south, east, and west (N, S, E, W)). This is similar to a screening procedure that is often applied using worst case scenarios in dispersion modeling [47], where, however, the CFD simulations are absent. For comparison purposes, the results are presented in a non-dimensional form relative to V₁₀, which is just above the building height. Before looking at the CFD calculations, it is important to note that, in spite of the favorable results of the validation [24], some discrepancies are evident between the SODAR measurements and the mesoscale numerical results, in terms of wind speed as well as wind direction. We chose to implement the mesoscale results, instead of the measurements, for creating the analytic expressions for the inlet profiles, as this would be more representative of a situation where the methodology is applied without available measurements. As shown in the graphs, the imposed inlet profiles (marked as inlet_B, C, D, or E) are in good agreement with the mesoscale results, indicating that they may be considered as an alternative to an interpolation procedure. A marked effect is observed when the wind direction is altered, keeping the inlet profile and the atmospheric condition the same. The development of the boundary layer profile within the CFD domain is highly dependent on the wind direction, obviously due to the distance of the measurement point 2 from each boundary, as well as the non-uniform distribution of the buildings constituting the urban geometry (Figure 1). It is the northern and eastern wind profiles that change the least, since they are closer (north), or approach from a less-dense urban area (east). The southern profile is the most affected, approaching from a denser part of the city and exhibiting the highest shear near the ground, thus reducing the ground velocity values and affecting the vertical distribution of the non-dimensional velocity. Information regarding the surface conditions is passed to the calculated inlet profiles through the roughness length (y_o) and the calculation of the friction velocity u_* from the velocity value at a specified height, given by measurements or mesoscale calculations. However, the results in Figure 4 are a clear indication of how sensitive the development of the boundary layer is to the local distribution of urban geometry, and how important its representation at a high spatial resolution becomes in microscale modeling efforts.

In Figure 5, Figure 6 and Figure 7, the effect of the turbulence model is examined, this time interpolating the results from mesoscale simulations to create inlet boundary conditions (marked as “TAPM-inlet” in the figures) for the CFD simulations. It is interesting to note that there is little effect on the profile from the inlet to the measurement point for both wind speed and temperature (Figure 5, Figure 6 and Figure 7). Apart from the application of the top boundary condition, in order to preserve the profiles [12,13], this is also due to the fact that the wind direction is mostly from the north in the mesoscale calculations (Table 1), which corresponds to the shortest distance from the measuring point, and therefore minimizes the effect of the urban fetch.

From the comparison of the results between the two turbulence models, it seems that there is negligible difference in either the mean wind speed profile or the temperature profile, but there is a significant difference in the results from the calculation of turbulence kinetic energy (Figure 6). Here, the inlet turbulence kinetic energy profiles are not conserved up to the measurement point for either of the two turbulence models. It is obvious that the Dirichlet-imposed boundary condition retains the inlet turbulence quantities along the upper boundary, but there is a shortcoming with the generation of turbulence at the bottom boundary. Although one would have expected the absence of imposed shear stress to lead to an even greater deficit of turbulence kinetic energy near the ground, it seems that the presence of the buildings generates a significant amount of turbulence and the overall level remains relatively high. This has been documented before in the literature [11] and seems to be verified here, even for the wind direction with the least effect of the buildings before reaching the measurement point. An exception to this is the highly unstable atmospheric condition (B), where the simulated turbulence kinetic energy production, due to mechanical shear and buoyancy, is not enough to retain the inlet value predicted by the mesoscale model. With regard to the two turbulence models, the differences are more pronounced here and although the k-ε model seems to manage to retain higher levels of turbulence, this should be observed with caution, as it might be attributed to the well-established shortcoming of this model to over predict turbulence kinetic energy production in stagnation regions on bluff bodies [13], such as those that arise here for the wind flow past buildings. Further effort needs to be put into this aspect of the modeling procedure, perhaps with some modification of the turbulence model parameters [6]. The SST model has been preferred for the rest of the calculations as, although no clear advantage could be discerned, it is the model with the stronger physical basis.

A comparison between the use of the analytic profiles, according to the COST 710 methodology [12], is summarized in Equations (2)–(11), and the interpolation of the mesoscale results as inlet boundary conditions is presented in Figure 8, Figure 9, Figure 10 and Figure 11 for the four chosen atmospheric conditions. The measured values are also shown on the graphs, where available. The mesoscale information (TAPM results) that was used for creating the COST 710 inlet profiles was the wind speed value at y = 10 m, the temperature value at y = 100 m, and the stability parameters (L_*, Ri) from Table 1.

As observed from Figure 8, Figure 9, Figure 10 and Figure 11, the analytic inlet profiles for wind speed and temperature agree well with the mesoscale results at these points (y = 10, 100 m, respectively), but in the case of the wind speed, the profiles diverge further away from the ground, while the temperature profile is hardly affected. Up to a height of y = 20–30 m, i.e., 3–4 mean building heights, the differences in the results for the wind speed profile are quite small, which is quite promising, since this is the region that is most affected by the urban roughness.

The effect of the inlet profiles on turbulence is interesting. In Figure 8—the neutral condition—turbulence is preserved near the top boundary for both types of inlet conditions. However, the constant turbulence levels prescribed by the COST analytic expressions are not preserved, as the discretely simulated urban roughness does not generate enough turbulence near the ground. On the other hand, although the high level of turbulence near the ground, applied through the mesoscale results, is also reduced, some of this energy is transferred higher above the ground and actually increases the turbulence levels there compared to the inlet.

In Figure 9 and Figure 10—the unstable conditions—the situation with regard to the mean wind speed and temperature is the same as in the neutral case, but, here, the inlet turbulence values are not sustained in either case, as it appears that a major turbulence production mechanism is missing from the CFD modeling procedure. This is more pronounced in Figure 9 and Figure 10, and is most probably related to insufficient modeling of turbulence production, due to buoyancy [37].

Under stable conditions (Figure 11), where it is the mechanically generated turbulence that dominates [41], there is an overall increase in turbulence kinetic energy compared to the inlet values of both the mesoscale results and the analytic expressions, most probably because of the high values of shear at the inlet, and possibly due to an inadequate modeling of turbulence suppression from stratification effects. The region very close to the ground (below 20–30 m) (Figure 11c) is a possible exception, where, again, turbulence production from the buildings cannot sustain the inlet values.

Overall, the discrepancies between the CFD calculated profiles over the urban region and the experimental measurements are of the same order as those of the mesoscale calculations. When focusing on the region close to the ground, for all the atmospheric conditions examined here, it seems that the temperature and wind speed profiles are well represented (Figure 8a,b, Figure 9a,b, Figure 10a,b and Figure 11a,b) by implementing the COST analytic expressions instead of profiles that are directly interpolated from the mesoscale results. For turbulence production, further effort is needed in order to correctly model buoyancy-related turbulence, but it seems that in the immediate vicinity of the buildings, it is the turbulence generated by the flow around them that dominates, and so some leverage may be permissible for the inlet values.

When coupling mesoscale and microscale simulations, the main goal is to take advantage of the microscale simulation for increased spatial resolution and detail, while relying on the mesoscale simulation to realistically provide for the atmospheric conditions. In Figure 12, Figure 13 and Figure 14, contours of pressure and temperature, along with wind speed vectors, are plotted from the microscale simulations, also in a close-up for the area where the experimental measurements took place (location 2 of Figure 1). It should be kept in mind that the whole area depicted in Figure 12, Figure 13 and Figure 14 corresponds to roughly four points of the mesoscale simulation in the horizontal plane and five in the vertical direction. The three figures correspond to three different atmospheric conditions (B, D, E of Table 1), and the effect of these is immediately evident in the pressure contours and the velocity vectors.

The pressure pattern shown in Figure 12, for the neutral atmospheric condition, clearly shows the stagnation regions appearing on the upstream surfaces of the buildings, not only in the first row of buildings exposed to the NW wind, but also for buildings in large open spaces (see the left side of the insert), where channeling of the flow from upstream buildings (see the right side of the insert) leads to an increase in flow momentum. Low-pressure regions on building roofs and corners (light-blue colored) are also evident on several buildings.

Comparing the results from different atmospheric conditions, the higher pressure regions of Figure 14 arise from the high shear and wind speeds appearing near the buildings at this stable atmospheric condition E. On the other hand, the neutral condition D (Figure 12), with the same value of wind speed at y = 10 m (Table 1), has markedly lower pressure values near the buildings, as noted especially in the close-ups of both Figure 12 and Figure 14. Even between the B and E conditions, where the wind direction is the same (Figure 13 and Figure 14, respectively), the difference in wind speed (2.6 m/s and 3.5 m/s for B and E, respectively) leads to pressure differences of up to ~30% on the building where the measurements took place.

The pressure and wind speed distribution shown in Figure 12, Figure 13 and Figure 14 cannot be reproduced through mesoscale simulations, and is of critical importance to several applications related to the microclimate of a specific neighborhood, such as cooling loads or the natural ventilation potential of a given building, urban canyon ventilation for air quality, small wind turbine siting, etc.

4. Conclusions

Microscale simulations have been performed using a computational fluid dynamics (CFD) approach, for the calculation of building-scale urban wind patterns and thermal environments. This is a step in the development process for the provision of detailed data and information required as input in building-scale applications that rely on the microclimate. The simulations were performed based on a one-way coupling to mesoscale simulations, with the information from the mesoscale applied to generate analytic expressions for the boundary conditions required by the CFD model.

The steady-state CFD simulation was performed at a 5 m horizontal and 1 m vertical spatial resolution, over a large urban area of the city of Kozani, Greece. For turbulence modeling, the standard k-ε and k-ω SST turbulence models were applied. Mesoscale numerical weather calculations were performed with the meteorological module of the atmospheric pollution model (TAPM), at a spatial resolution of 300 × 300 m (inner grid), and the first grid point above ground at 10 m. The measurements of the vertical profiles of wind speed and temperature have been provided by SODAR and a meteorological temperature profiler.

One of the main points of focus of the study was the method of passing mesoscale information to the microscale simulation procedure, i.e., the boundary conditions for the CFD. In a novel approach, previously developed analytic expressions, based on the Monin–Obukhov theory, were used to generate profiles, dependent on minimal mesoscale information, for the vertical distribution of wind speed, temperature, turbulence kinetic energy, and its dissipation rate. This procedure was compared to the conventional approach of directly interpolating data from the coarse mesoscale grid. The measurements of wind speed and temperature were performed for the purpose of the study, and were also used to evaluate the approach. Two turbulence models were applied in their standard implementation, allowing for the discrete representation of the building geometries, to account for turbulence processes at the microscale, rather than modifying the models to account for atmospheric conditions.

The advantage of the new approach is that locally relevant mesoscale information is passed to the microscale simulation based on a minimal set of index values (V_y=10, T_y=100, L_*), and so parametric studies and hypothetical situations may be studied, without relying on a cumbersome, direct, real-time coupling with the mesoscale calculations. The present approach cannot substitute such a direct coupling, especially in terms of realistically reproducing transient physics, as it is highly dependent on the Monin–Obukhov similarity theory. For example, inversions, nocturnal jets, and other outlier phenomena, such as hurricanes, etc., cannot be taken into account. Furthermore, in order to ensure applicability of the model, some consistent mesoscale information must be available in the form of the three index values (V_y=10, T_y=100, L_*). If these are not compatible amongst themselves, or if a transient or extreme meteorological event is present, the theory upon which the analytical expressions are based will fail. However, given the consistent set of these values, the majority of basic atmospheric conditions are easily accounted for, and so the method does provide flexibility as a design tool for microscale studies.

The application of the method proved promising. Mean wind speed and temperature profiles were found to be in close agreement with the direct interpolation of data from the mesoscale (Figure 4 and Figure 8, Figure 9, Figure 10 and Figure 11). Turbulence kinetic energy production proved to be more of a challenge, and the hypothesis that building-generated turbulence will dominate seems to depend highly on the choice of turbulence model. The results in Figure 6 show that the k-ε model produces more turbulence near the ground than the SST model, even though their inlet values are the same. From the results in Figure 8c, Figure 9c, Figure 10c and Figure 11c, it seems that turbulence far from the ground is consistently underestimated, except for in the case of stable atmospheric conditions (Figure 11c). This is most probably an indication of under-estimation of the damping occurring in the stable case, while in the unstable and neutral cases, it seems that ground level turbulence is underestimated and buoyancy does not generate or transfer turbulence to higher altitudes (Figure 9c and Figure 10c). On the other hand, mechanical turbulence, which dominates in the case of the neutral conditions (Figure 8c), does seem to be faring better at keeping up with the inlet values, although it still underestimates them. The representation of atmospheric turbulence further from the ground needs more modeling effort, possibly modifying the turbulence models in way that will not affect the local flow characteristics around buildings.

Further tuning and development of the turbulence modeling for non-neutral conditions is obviously one of the next steps of development. Furthermore, turbulence model modifications for the upstream fetch before the presence of the buildings is an active area of research. We have already (results yet to be published) successfully applied the methodology to rural areas, where the area of interest is a small cluster of buildings and the upstream fetch is less affected by microscale topography, e.g., clusters of buildings on islands. Nevertheless, even at the present stage of implementation, the new approach proves to be a useful tool for microscale studies in urban areas as well. It provides flexibility in defining different atmospheric states and does not rely on a mesoscale simulation, as long as preliminary parametric studies are being performed with regard to wind direction, atmospheric stability conditions, etc. The results show that the inclusion of mesoscale information in the form of (V_y=10, T_y=100, L_*) is important, as it may lead to significant variations in critical building-related design parameters, such as surface pressure distributions and local wind patterns.