Steady RANS model of the homogeneous atmospheric boundary layer

Journal of Wind Engineering & Industrial Aerodynamics 173 (2018) 289–301

Contents lists available at ScienceDirect

Journal of Wind Engineering & Industrial Aerodynamics journal homepage: www.elsevier.com/locate/jweia

Steady RANS model of the homogeneous atmospheric boundary layer Mihael Cindori, Franjo Jureti
c, Hrvoje Kozmar *, Ivo Dzijan Faculty of Mechanical Engineering and Naval Architecture, University of Zagreb, Ivana Luci
ca 5, 10000 Zagreb, Croatia

A R T I C L E I N F O

A B S T R A C T

Keywords: Homogeneous atmospheric boundary layer Neutral thermal stratification Flat terrain Cubic building model Surface pressure distribution Steady Reynolds-averaged Navier-Stokes equations Standard k-ε turbulence model Source term modeling Computational wind engineering

A computational model was proposed to simulate the homogeneous atmospheric boundary layer (ABL) above rural, suburban and urban flat terrains for the neutral thermal stratification of the atmosphere. This is achieved using the steady Reynolds-averaged Navier-Stokes equations and the standard k-ε turbulence model. The major feature of the proposed approach is the implementation of the wind-source term in the momentum equation, which is adjusted using the experimental wind-tunnel data. Profiles of the mean velocity, turbulent kinetic energy and Reynolds shear stress are homogeneous along the computational domain. The mean discrepancy in the longitudinal direction is 0.1% for the mean velocity profile, 0.6% for the Reynolds shear stress profile, 0.3% for the turbulent kinetic energy profile. The preliminary results show good agreement between the computed profiles and the ABL wind-tunnel simulations, as the mean discrepancy between the experimental and computational results is 4.5% for the mean velocity profile, 2.5% for the Reynolds shear stress profile, 6% for the turbulent kinetic energy profile. The surface pressure distribution at the cubic building model subjected to the suburban ABL flow indicates characteristic patterns in accordance with the respective flow phenomena, i.e. flow separation and reattachment, stagnation zone, leeward suction. The results show that the calculated wind-source term does not have a significant influence on the aerodynamic loads on the cubic building model, whereas it ensures homogeneity of the boundary layer flow properties throughout the computational domain in agreement with the wind-tunnel results.

1. Introduction Homogeneity of the atmospheric boundary layer (ABL) is a key issue in computational simulations of wind engineering problems. Previous studies investigated this issue from various perspectives, such as the modification of boundary conditions and inflow profiles, or by using a precursor domain. The current study applies a source-term modeling technique to ensure homogeneity throughout the computational domain, as well as to ensure boundary layer characteristics that closely agree with those obtained in an experiment. In their seminal work, Richards and Hoxey (1993) proposed inflow velocity and turbulence profiles for the turbulent kinetic energy and dissipation rate, that together with suggested near-wall treatment and constant shear stress value at the top of a domain ensure homogeneous flow and zero pressure gradient in the computational domain. However, the calculated turbulent kinetic energy profiles are constant along the height of the domain, which is not consistent with available wind-tunnel measurements. Hargreaves and Wright (2007) confirmed a necessity of setting the shear stress boundary condition at the top of the computational domain

to ensure homogeneous flow in their computational model of the neutrally stratified atmospheric boundary layer. Yang et al. (2009) proposed a set of inflow boundary conditions derived from the k-ε turbulence model transport equations that provide homogeneous ABL for the rural terrain, while O'Sullivan et al. (2011) suggested to use periodic boundary conditions in order to generate inflow profiles consistent with the rough wall function. They proposed the boundary conditions at the top of the domain for turbulent kinetic energy and dissipation rate. Richards and Norris (2012, 2013, 2015) concluded that without the defined shear stress at the top of the domain, the atmospheric flow is driven only by the pressure drop that can produce serious streamwise inhomogeneities. Consequently, they derived mean velocity and turbulent inflow profiles for the standard k-ε and SST k-ω turbulence models that ensure streamwise homogeneity for the equilibrium pressure-driven boundary layer, as well as the decay of turbulent kinetic energy along the height of the domain. Parente et al. (2011a, 2011b) introduced an approach based on the modification of the standard k-ε turbulence model equations to allow specification of an arbitrary inflow profile for the turbulent kinetic energy. Balogh and Parente (2015) further improved the methodology by

* Corresponding author. E-mail address: [email protected] (H. Kozmar). https://doi.org/10.1016/j.jweia.2017.12.006 Received 28 June 2017; Received in revised form 4 December 2017; Accepted 4 December 2017 Available online 8 January 2018 0167-6105/© 2017 Elsevier Ltd. All rights reserved.

M. Cindori et al.

Journal of Wind Engineering & Industrial Aerodynamics 173 (2018) 289–301

Nomenclature

z

Cp Pressure coefficient Cμ, C1, C2, σ k, σ ε Constants of the k-ε turbulence model eavg Mean percentage error ei Percentage error Iu Turbulence intensity k(z) Turbulent kinetic energy kC Turbulent kinetic energy in the near-wall cell L Length of the computational domain ntotal Total number of sampling points Pk(z) Turbulence production term PC Turbulence production in the near-wall cell peff Effective kinematic air pressure pinlet Pressure at the inlet of the computational domain poutlet Pressure at the outlet of the computational domain pref Reference pressure uðzÞ, vðzÞ, wðzÞ Mean velocity components in the x-, y-, z-direction uC Mean velocity value in the near-wall cell center u', v', w' Fluctuating velocity components in the x-, y-, z-direction -ρw0 u0 ðzÞ Turbulent Reynolds shear stress uref Mean reference velocity in the reference height Uref Reference velocity uτ Friction velocity WSx(z) Longitudinal component of the wind-source term WSinit x ðzÞ Initial longitudinal component of the wind-source term WSnew x ðzÞ Wind-source term in the current iteration WSold x ðzÞ Wind-source term in the previous iteration x Distance in the main flow direction y Spanwise distance from the test-section centerline

zc zref z0 z0f

α

δ δP ε(z)

εc

κ

μ ν νeff νT ðzÞ ρ τzx ðzÞ τCFD zx ðzÞ τexp zx ðzÞ τold zx ðzÞ τw

ϕi ϕmajor ϕminor

ψ

Vertical distance from the wind-tunnel floor and the bottom of the computational domain Distance between the near-wall cell center and the wall Reference height Aerodynamic surface roughness length in computational simulations Aerodynamic surface roughness length in full-scale Power-law exponent Atmospheric boundary layer thickness Pressure-gradient length scale Turbulent kinetic energy dissipation Turbulence dissipation in the near-wall cell Von K
arm
an constant Dynamic laminar viscosity Kinematic viscosity Effective kinematic viscosity Kinematic turbulent (eddy) viscosity Air density Shear stress Computed Reynolds shear stress Measured Reynolds shear stress Shear stress in the previous iteration Shear stress at the wall Value of the flow property Major value of the flow property Minor value of the flow property Ratio between the wind-source magnitude and the pressure gradient magnitude

Balogh et al. (2012) modified the rough wall function together with the turbulence model to provide homogeneous conditions to simulate the ABL flow over a hilly object. A theoretical study on the k-ε turbulence model with respect to the standard Reynolds shear stress distribution in the wind tunnel was studied by Jureti
c and Kozmar (2013). The results showed that the standard turbulence models can reproduce the ABL when subjected to appropriate loading, which is based on controlling the pressure drop in the domain to obtain the correct velocity at the measurement point. The present study focuses on developing a simple and robust methodology for Computational Fluid Dynamics (CFD) simulations of the homogeneous ABL, where the profiles of the relevant ABL parameters alter negligibly along the computational domain. This is achieved by using the standard k-ε turbulence model and steady RANS equations. The present study is based on the incorporation of the wind-source term in the momentum equation when performing computational simulations of the neutrally stratified ABL, as originally suggested by Jureti
c and Kozmar (2013), which approach was subsequently justified by Cai et al. (2014). The computational model in the present study is adjusted to enforce the correct Reynolds shear stress distribution in agreement with the wind-tunnel experiments. This approach yields the appropriate profiles of the mean velocity, turbulent kinetic energy and Reynolds shear stress in the created computational simulations of the neutrally stratified ABL. Preliminary 2D computational simulations were performed for the rural, suburban and urban type of terrain, by using the von Neumann boundary condition at the inlet and the outlet of the domain in order to ensure zero gradients of the characteristic flow parameters that are in agreement with the experimental results of the ABL wind-tunnel simulations. This approach is further tested by preliminarily studying surface pressure distribution at the cubic building model subjected to the suburban ABL simulation. The influence of the wind-source term on building

defining a realistic inflow profile for the turbulent kinetic energy. Yang et al. (2017) proposed a new set of inflow boundary conditions for the SST k-ω model. These conditions ensure homogeneity of the flow in the computational domain. The precursor domain technique is often used to recreate atmospheric conditions in an empty domain prior to computational study of obstacles placed in the ABL flow. For instance, Pontiggia et al. (2009) performed a simulation with the periodic conditions to generate appropriate atmospheric conditions for gas dispersion investigation, while O'Sullivan et al. (2011) used precursor simulations in an empty domain with applied periodic conditions to overcome streamwise gradients provoked by the inconsistencies between the boundary conditions, the inflow profiles and the wall function. Jureti
c and Kozmar (2014) used wind tunnel geometry with von Neumann boundary conditions applied at the inlet and the outlet boundaries to obtain homogeneous mean velocity, turbulent kinetic energy and Reynolds shear stress profiles with various RANS turbulence models. Blocken et al. (2007a) investigated the influence of the sand-grain roughness coefficients in the wall function on the streamwise homogeneity. They propose to always perform sensitivity tests in an empty domain to avoid homogeneity uncertainties. Moreover, Blocken et al. (2007b) concluded that the sand-grain roughness wall functions are not suitable for ABL modeling. They also showed that the streamwise gradients, which are mainly the result of flow inhomogeneity, can be influenced by the inappropriate wall-function formulation, near-wall grid resolution, roughness height and the applied boundary conditions. Parente et al. (2011a) proposed a modified version of the Richards and Hoxey (1993) rough wall function based on the aerodynamic roughness length. The new formulation allows more flexibility in the near-wall mesh generation. 290

M. Cindori et al.

Journal of Wind Engineering & Industrial Aerodynamics 173 (2018) 289–301

(1)–(4),

aerodynamics is investigated. The obtained pressure distribution is compared in simulations with and without the wind-source term. 2. Computational model The model equations are developed using the incompressible steady RANS equations, the standard k-ε turbulence model (Jones and Launder, 1972), and the Boussinesq (1877) hypothesis,

∂ui ¼ 0; ∂xi

(1) 




∂ui 1 ∂peff ∂ ∂ui ∂uj ðν þ νT Þ uj ¼ þ þ ∂xj ρ ∂xi ∂xj ∂xj ∂xi uj

∂k ∂ ¼ ∂xj ∂xj

∂ε ∂ uj ¼ ∂xj ∂xj




νþ



 ;






(2)

(4)

k2

ε

(5)



Pk

¼ νT

ρ

;

∂ui ∂uj þ ∂xj ∂xi



(13)

νT ðzÞ

duðzÞ τzx ðzÞ ¼ : dz ρ

(15)

(16)

(6)

Eqs. (15) and (16) indicate that the preliminary wind-source term calculation can be performed as a shear stress derivation as follows,

(7)

WSinit x ðzÞ ¼



∂ui : ∂xj

(14)

where,

The following homogeneity conditions are adopted for the model equations:

1 d τexp zx ðzÞ ; ρ dz

(17)

 The thermal stratification of the atmosphere is adopted to be neutral.  The convection terms for the stationary ABL flow are zero,

where the initial wind-source term value is calculated using the data reported in Kozmar (2011a) to reproduce the experimental Reynolds shear stress wind-tunnel data τexp zx ðzÞ. Thus Eq. (17) represents a definition of the wind-source term implemented in the momentum equation. If the wind-source term WSx(z) in Eq. (15) is set to zero, Eq. (15) is reduced to the form of the momentum equation derived by Richards and Hoxey (1993). The wind-source term in Eq. (11) allows the user to control the distribution of the flow properties by adjusting the source until the computed shear stress τzx ðzÞ matches the experimental values τexp zx ðzÞ to the required accuracy. The properly adjusted wind-source term should be able to match the experimental measurements, as well as provide no pressure gradient in the computational domain. Eqs. (3) and (4) reduce to Eqs. (12) and (13), and the turbulent viscosity is calculated as,

∂ui ¼ 0; ∂xj

νT ðzÞ ¼ Cμ

 Since the mean ABL flow is predominantly one-dimensional (e.g. Dyrbye and Hansen, 1997), the dominant velocity component is the horizontal uðzÞ velocity and the significant derivatives are with respect to the vertical direction only, x1 ¼ x; x2 ¼ y; x3 ¼ z; u1 ¼ uðzÞ; u2 ¼ v ¼ 0; u3 ¼ w ¼ 0; kðzÞ; εðzÞ; νT ðzÞ; Pk ðzÞ:

(8)

Accordingly, Eq. (1) is beforehand satisfied.

uj

(9)

 A new source term WSx(z) is introduced to achieve the distribution of the flow properties which is in close agreement with the available measurements.  The viscous stress is neglected, as the turbulent Reynolds shear stress largely prevails for the most part of the computational domain, 

2  duðzÞ Pk ðzÞ ¼ νT ðzÞ : dz



∂ui ∂uj þ ¼ 0: ∂xj ∂xi

kðzÞ2 : εðzÞ

(18)

Production of the turbulent kinetic energy Pk is modeled using Eq. (7). After the Boussinesq hypothesis is introduced and by assuming the homogeneous turbulence (turbulence production is equal to turbulence dissipation), Eq. (7) reduces to,

and all streamwise gradients are neglected.

ν

  d νT ðzÞ d εðzÞ Pk ðzÞεðzÞ εðzÞ2  C1 þ C2 ¼ 0:  dz σ ε dz ρkðzÞ kðzÞ

  d duðzÞ  ¼ WSx ðzÞ; νT ðzÞ dz dz

Eq. (4) is the transport equation for the turbulence dissipation rate ε, where the turbulent viscosity νT and the production term Pk are respectively modeled as,

νT ¼ C μ

(12)

The wind-source term is implemented in the momentum equation to maintain the inlet flow profiles throughout the computational domain, i.e. to allow for the appropriate distributions of the relevant flow properties, e.g. velocity, shear stress and turbulent kinetic energy. Eq. (11) is a starting point for the wind-source term calculation. By taking into account Eq. (14), Eq. (11) reduces to,

Eqs. (1) and (2) are continuity and momentum equations, respectively, and Eq. (3) is the transport equation for the turbulent kinetic energy k defined as, 1 k ¼ u0 i u0 i : 2

  d νT ðzÞ dkðzÞ Pk ðzÞ   þ εðzÞ ¼ 0; dz σ k dz ρ

dpeff ¼ 0: dx

(3)



ν ∂ε εPk ε2 νþ T  C2 : þ C1 σ ε ∂xj kρ k

(11)

Eqs. (11)–(13) ensure homogeneous ABL flow and the zero pressure gradient is assumed in order to determine the initial wind-source term value,



νT ∂k Pk þ  ε; σ k ∂xj ρ

  d duðzÞ 1 dpeff  þ ¼ WSx ðzÞ; νT ðzÞ dz dz ρ dx

(19)

In order to complete the standard k-ε turbulence model, the values of the model constants need to be determined. As given in Pope (2000), the Cμ constant represents the ratio between the Reynolds shear stress and the turbulent kinetic energy. For simple turbulent shear flows, under the

(10)

The following model equations are accordingly developed from Eqs.

291

M. Cindori et al.

Journal of Wind Engineering & Industrial Aerodynamics 173 (2018) 289–301

assumption of the homogeneous turbulence, Cμ can be calculated using the following expression,

 du PC ¼ τw  ; dz C

kðzÞ 1 ¼ pffiffiffiffiffiffi : w0 u0 ðzÞ Cμ

   where du dz  is the velocity gradient value in the cell center near the wall,

(20)

C

while the Reynolds shear stress at the wall τw is modeled as follows,

Jureti
c and Kozmar (2013) used the data measurements from Kozmar (2011a) and showed that the ratio between the turbulent kinetic energy and the Reynolds shear stress on the left hand side of Eq. (20) does not change significantly with height for rural, suburban and urban terrains, and this value is equal to 4.76 in all measurement points. The value of the Cμ constant is calculated using Eq. (20) and it is equal to 0.044. Richards and Hoxey (1993) noted that the inlet flow profiles derived for homogeneous flow conditions result in the Reynolds shear stress that remains constant with increasing height. The inlet velocity profile has a logarithmic profile that is maintained throughout the computational domain if the equation for the turbulence dissipation is satisfied, and it is achieved by calculating the turbulence constant σ ε as follows,

σε ¼

κ2 pffiffiffiffiffiffi ; ðC2  C1 Þ Cμ

τw uC uτ κ : ¼  ρ ln zC þz0 z0

3. Computational setup Computational modeling of the ABL in rural, suburban and urban terrain is carried out using an open-source CFD toolbox OpenFOAM®. The experimental ABL profiles in Kozmar (2011a) were used for calibration of the present CFD simulations. More experimental details on wind-tunnel geometry, flow development and uniformity in the wind-tunnel test section are reported in Kozmar (2011a, 2011b), effects of surface roughness and barrier wall in Kozmar (2008, 2012a), downscaling issues in Kozmar (2010, 2011c, 2012b) and pressure drop control in Kozmar (2005). Unlike the wind-tunnel test section with the ceiling adjustable in height (Kozmar, 2011a), the computational domain is designed with a constant height top boundary at z ¼ 1 m to demonstrate that the model is capable of resolving a homogeneous ABL in that situation, Fig. 1. The 2D simulations of the rural, suburban and urban ABL flow were performed in an empty computational domain to adjust the value of the wind-source term WSx(z). As shown in Eq. (17), the preliminary methodology to calculate the initial wind-source values is based on the experimental Reynolds shear stress data. These values were additionally adjusted during the simulations by the correction procedure implemented in the CFD code,

(21)

u3τ ; εC ¼ κðzC þ z0 Þ

old WSnew x ðzÞ ¼ WSx ðzÞ

(22)

(23)

where kC ¼ k(zC) represents the value of the turbulent kinetic energy in the adjacent cell center. The turbulence production term PC ¼ P(zC) in the wall function is set as,

Table 1 Constant values of the k-ε turbulence model used in previous works (standard constants) and the present study (model constants).

Standard constants Model constants

Cμ

C1

C2

σε

σk

0.09 0.044

1.44 1.44

1.92 1.92

1.3 1.67

1 1

τexp zx ðzÞ ; τold zx ðzÞ

exp ktop ¼ ktop ;

(27)

the turbulent viscosity is calculated using, τexp zx;top ρ

νT;top ¼  exp ; du dz

top

Table 2 Aerodynamic surface roughness length in performed CFD simulations (z0) and respective full-scale measures (z0f). Terrain exposure

z0, m

z0f, m

Rural Suburban Urban

1.70∙103 4.00∙103 15.5∙103

0.23 0.33 2.47

(26)

where WSnew x ðzÞ is the wind-source value in current iteration and WSold x ðzÞ is the wind-source value from the previous iteration. Therefore, when shear stress value from the previous iteration τold zx ðzÞbecomes equal to the experimental shear stress τexp zx ðzÞ, iterative correction is completed. The correction procedure given in Eq. (26) reduces differences between the computed shear stress and the experimental results, as the correct distribution of the shear stress is a prerequisite to obtain the appropriate distribution of the velocity and turbulent kinetic energy. Four boundary conditions are defined at the boundaries of the computational domain, Fig. 2. The top of the computational domain is pre-defined using the experimental data from Kozmar (2011a). In particular, the turbulent kinetic energy is set to the experimental value,

where z0 is aerodynamic surface roughness length adopted from Jureti
c and Kozmar (2013, 2014). The values are provided in Table 2 together with the corresponding values of the full-scale aerodynamic surface roughness length z0f (Kozmar, 2011a). The friction velocity uτ is implemented as, p ffiffiffiffiffiffipffiffiffiffiffi 4 Cμ kC ;

(25)

The velocity value in the adjacent cell center is uC ¼ uðzC Þ.

where C2 and C1 are model constants, κ is von K
arm
an constant equal to 0.41. The wall conditions presented by Richards and Hoxey (1993) were used to model the flow in the near-wall region in the present study as well. The logarithmic velocity profile is employed close to the ground, hence, the relation in Eq. (21) was used to determine the value of the σ ε constant. The same relation is proposed by Pope (2000) to reproduce the flow behavior in the log-law region. The standard values commonly used in previous studies are reported in the first row of Table 1 along with the novel model constants adopted for the present study, which are reported in the second row. The ABL flow in the near-wall region is modeled as suggested by Richards and Hoxey (1993), so dissipation rate εC ¼ ε(zC) in the adjacent cell center zC is modeled as follows,

uτ ¼

(24)

Fig. 1. Schematic view of the computational domain geometry. 292

(28)

M. Cindori et al.

Journal of Wind Engineering & Industrial Aerodynamics 173 (2018) 289–301 Table 3 Boundary conditions.

and the dissipation rate is determined using an expression provided in the standard k-ε turbulence model, 2 ktop

εtop ¼ Cμ

νT;top

Rough wall

! :

(29)

uðzÞ ¼ uref



z

α

zref

dk dz

k

The von Neumann zero gradient boundary condition is used at inflow and outflow boundaries to eliminate streamwise gradients, pressure is set to zero at the outlet and the constant shear stress is set to the top boundary as the full depth of the ABL is not simulated, with all the details of the boundary conditions adopted for the computational domain reported in Table 3. Internal field velocity profiles are initialized using the power-law, Hellman (1916), :

Top

εC ¼

p

dpeff dz

u3τ κðzC þz0 Þ

¼0

Outlet

du dx

¼0

du dx

¼0

k ¼ kexp top ε ¼ εtop

dk dx dε dx

¼0

dk dx dε dx

¼0

dpeff dz

dpeff dx

du dz

¼0

ε

Inlet  exp

u¼0

u

¼

du dz

top

¼0

¼0 ¼0

¼0

peff ¼ 0

Table 4 Characteristic parameters of the mean velocity profiles.

(30)

Terrain exposure

α, -

zref , m

uref , m/s

uτ, m/s

Rural Suburban Urban

0.16 0.20 0.37

0.202 0.202 0.202

14.97 13.48 10.14

1.11 1.10 1.43

The values of the power-law exponent α, the reference velocity uref and the friction velocity uτ for the rural, suburban and urban terrain are adopted from wind-tunnel experiments (Kozmar, 2011a), Table 4. A coupled system of the algebraic equations is solved using the SIMPLE algorithm, Patankar and Spalding (1972). Convection terms are approximated using the second-order discretization Gamma scheme (Jasak, 1996) and diffusion terms are modeled using the second-order Gaussian-based discretization scheme, e.g. OpenFOAM user guide (2016). The finite volume mesh is generated using blockMesh, the standard mesh utility in OpenFOAM®. As von Neumann boundary condition at the inlet and the outlet of the computational domain ensures homogeneous flow, only fifteen cells are generated in the streamwise direction. The resolution in the vertical direction is specified with respect to the aerodynamic surface roughness length z0. Hence, the resolution in the near-wall region differs for various types of terrain, since the distance of the near-wall grid layer centroid from the wall zC needs to be larger than z0, Blocken et al. (2007b). Thus zC for the rural terrain is equal to 1.932⋅103 m, for suburban terrain 4.5⋅103 m, and for urban terrain 1.6⋅102 m; these values are larger than the z0 values presented in Table 2. The resolution for rural, suburban and urban terrain types is the same as in Jureti
c and Kozmar (2013). In particular, 53 cells are generated in vertical direction for the rural terrain, 50 cells for the suburban terrain, and 21 cells for the urban terrain. The mesh resolutions from the bottom to the top of the computational domain for a single column of cells are presented in Fig. 3. The relative percentage error is calculated as a difference between the maximum and minimum value for each of the studied flow parameters divided by the maximum value,   jmaxðϕi Þj  jminðϕi Þj ⋅100%; ei ¼   jmaxðϕ Þj

(31)

i

Fig. 3. Mesh resolutions in the vertical direction (from the bottom to the top of the computational domain): a) rural terrain, b) suburban terrain, c) urban terrain.

and the relative mean percentage error eavg is calculated as, nP total

eavg ¼

ei

i¼1

ntotal

;

performed by calculating the absolute values of the relative error as,

(32)

ei ¼

where ntotal is the total number of sampled points. A comparison between the experimental and computational results is

jϕjmajor  jϕjminor ⋅100%; jϕjmajor

(33)

where ϕmajor (larger value) and ϕminor (smaller value) are the computational and the experimental values, respectively. The mean relative percentage error eavg is calculated using Eq. (32). 4. Results and discussion The computational results are reported up to 1 m height of the

Fig. 2. Schematic view of the boundary patches in the computational domain. 293

M. Cindori et al.

Journal of Wind Engineering & Industrial Aerodynamics 173 (2018) 289–301

Reynolds shear stress, and turbulent kinetic energy is illustrated in Fig. 5. The maximum and mean relative errors are provided in Table 5. The reported CFD results are in good agreement with the fundamental characteristics of the full-scale ABL, as reported in Stull (1988), Garratt (1992), Wyngaard (2010). The Reynolds shear stress and the turbulence intensity are largest near the ground surface. They decrease with increasing the height and eventually become nearly zero at some vertical distance from the ground surface. In general, the CFD results indicate a very good homogeneity along the computational domain, as the maximum discrepancy for all the parameters and terrain types does not exceed 1.5%. The maximum errors occur for the rural ABL, while the minimum errors are obtained for the urban ABL. These discrepancies are smallest for the mean velocity profiles and largest for the Reynolds shear stress profiles. This is likely due to a current near-wall treatment, as the mesh resolution in the vertical direction depends on the aerodynamic surface roughness length z0 and this

Fig. 4. Schematic view of the sample lines in the computational domain.

respective wind-tunnel ABL simulations. The mean velocity, Reynolds shear stress, and turbulent kinetic energy profiles are calculated and reported at the inlet (Sample line 1), center (Sample line 2), and outlet (Sample line 3) of the computational domain, as illustrated in Fig. 4. The 2D CFD simulations are analyzed with respect to the ABL homogeneity along the computational domain by comparing the computational and experimental results at the sampling lines. Homogeneity of the ABL normalized mean flow velocity, normalized

Fig. 5. 2D results for the normalized mean velocity, normalized Reynolds shear stress and turbulent kinetic energy profiles along the computational domain for the rural, suburban, and urban ABL simulations. Profiles are sampled at the sample line 1 - SL1 (o), sample line 2 – SL2 (x) and sample line 3 – SL3 (þ). 294

M. Cindori et al.

Journal of Wind Engineering & Industrial Aerodynamics 173 (2018) 289–301

computational domain (SL2) is presented in Fig. 6. Errors between the computational results with the applied wind-source term and the experimental results are reported in Table 6. The differences between the mean velocity profiles obtained with and without the wind-source term are small and those velocity profiles correspond to the respective power-law distributions. The profiles for the characteristic flow parameters obtained in simulations with and without the wind-source term are homogenous along the computational domain for all the three ABL simulations due to the von Neumann boundary condition applied at the inlet and the outlet of the computational domain. When the wind-source term is not applied, the wall function at the

Table 5 Maximum relative error and mean relative error (in brackets) of the relevant ABL parameters calculated along the computational domain. Parameter

Rural ABL

Suburban ABL

Urban ABL

uðzÞ=uref τzx ðzÞ=u2τ k(z)

0.21 (0.1) % 1.43 (0.6) % 0.72 (0.3) %

0.15 (0.1) % 1.41 (0.6) % 0.47 (0.3) %

0.19 (0.1) % 0.57 (0.3) % 0.41 (0.2) %

can affect how the computational data is fitted to the experimental data. A comparison of the computational results with and without the wind source term WSx(z) against the experimental results at the center of the

Fig. 6. Computational results obtained using the wind-source term (o), without the wind-source term (x), both captured in the center of the computational domain (SL2) in comparison with the experimental results (þ) and the respective power-law profiles () for mean velocity. 295

M. Cindori et al.

Journal of Wind Engineering & Industrial Aerodynamics 173 (2018) 289–301

the wind-source term is calculated using the standard deviation as,

Table 6 Maximum relative error and mean relative error (in brackets) between the experimental and computational results obtained with the applied WSx(z) and captured at the sample line 2 (SL2). Parameter

Rural ABL

Suburban ABL

Urban ABL

vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u N u1 X

CFD 2 σ¼t τ  τexp zx;i ; N i¼1 zx;i

uðzÞ=uref τzx ðzÞ=u2τ kðzÞ

8.2 (4.1) % 4.8 (1.6) % 4.5 (1.6) %

5.8 (3.0) % 5.0 (1.9) % 5.8 (2.9) %

5.6 (2.4) % 5.8 (2.2) % 21.1 (5.6) %

where N is the total number of the computational and experimental shear stress values with respect to the height z, τCFD zx;i is the calculated shear

(34)

stress value at the specific point i, and the τexp zx;i is the experimental shear stress value at the same point. The calculated standard deviations presented in Fig. 7 clearly indicate that the wind-source term profiles converge in all three ABL simulations. The converged wind-source profiles are reported in Fig. 8. The windsource term is calculated as a derivative of the experimental shear stress values. The obtained wind-source values are perturbed due to differentiation of the measured shear stress. Nevertheless, the obtained computational shear stress distributions still indicate that the provided windsource values can satisfactorily replicate the respective experimental results. The pressure gradient values obtained in the simulations with the applied wind-source term are constant for all three terrain types and calculated as (pinlet – poutlet)/L. pinlet is pressure value at the inlet of the computational domain, poutlet is value at the outlet of the computational domain and L is domain length. The calculated pressure gradient for the rural terrain is 1.36⋅102 Pa/m, for suburban terrain 5.1⋅102 Pa/m and

bottom surface and the constant shear stress at the domain top provide a constant shear stress and turbulent kinetic energy distribution along the height z. On the contrary, a decay of the shear stress and turbulent kinetic energy with increasing height can be noticed when the wind-source term is applied. When the wind-source term is applied, the agreement between the computational and experimental results is good, while the discrepancies are larger for the turbulent kinetic energy, which is particularly exhibited close to the wall for the urban ABL. However, the mean relative error is less than 6% for all the parameters and terrain types. The maximum relative error is less than 9% in all the tests except for the turbulent kinetic energy in the urban ABL, where it is 21%. The Cμ constant of the standard k-ε turbulence model, representing the ratio between the Reynolds shear stress and the turbulent kinetic energy, does not remain constant any more, and this forces over-estimation of the turbulent kinetic energy near the ground for the urban terrain type. During the iterative computational process, WSx(z) profiles are tuned using the correction procedure provided in Eq. (26). The wind-source profiles were corrected 24 times in each iteration. The convergence of

Fig. 9. Schematic view of the 3D computational domain with defined boundaries (all dimensions are in meters).

Fig. 7. Standard deviations calculated in the rural, suburban and urban ABL simulations.

Fig. 8. Adjusted WSx(z) profiles obtained in 2D calculations. 296

M. Cindori et al.

Journal of Wind Engineering & Industrial Aerodynamics 173 (2018) 289–301

5. Surface pressure distribution at the cubic building model

for urban terrain 0.26 Pa/m. Yaglom (1979) performed similarity and dimensional analysis for the constant-pressure and pressure-gradient turbulent wall flows using the pressure-gradient length scale δP, u2 δP ¼  τ  : 1 dp ρ dx

The 3D computation of the suburban ABL flow over the cubic building model is performed for the suburban terrain. The goal is to find out if the model ensures valid pressure distribution at cube surfaces, as well as if the wind-source term does have any significant influence on the flow in the region near the cube. Accordingly, the pressure distribution at the cubic building model is compared to the pressure distribution obtained in the simulation without the wind-source term and preliminary investigation on the wind-source term influence is presented. The distribution of the pressure coefficient at surfaces of the cubic building model is reported and discussed. The turbulence intensity Iu value of the ABL simulation with the applied wind-source term at the height h ¼ 0.2 m is approximately 20% for the suburban ABL (Kozmar, 2011a). The obtained results are qualitatively compared with the relevant previous studies, e.g. Castro and Robins (1977), Paterson and Apelt (1990), Lim et al. (2009). The problem is symmetric and hence half of the flow domain is

(35)

The pressure gradient does not have significant influence on the flow if the pressure-gradient length scale is significantly larger than the ABL thickness δ or some other characteristic vertical length scale. In our ABL simulations, the calculated pressure-gradient length scale is approximately 113 m for the rural terrain, 37 m for the suburban terrain, and 10 m for the urban terrain. The characteristic vertical length scale in the reported computational simulations of the ABL is the height of the computational domain, which is 1 m. Therefore, the pressure-gradient influence on the flow can be considered negligibly small, as the characteristic vertical length scale is much smaller than the calculated pressuregradient length scale for all tested terrain types, in agreement with Yaglom (1979).

Fig. 10. Finite volume mesh used in performed 3D computations: a) mesh generated for the suburban ABL flow over the cubic building, b) a detail of the mesh in the region near the cube.

Table 7 Boundary conditions adopted for the 3D simulation of the suburban ABL flow over the cube. Cube u

Rough surface

u¼0

u¼0

dk ¼ dy ¼ dk dz ¼ 0

dk dx

ε

Standard wall function

εC ¼

p

dpeff dx

dpeff dz

¼

dpeff dy

¼

dpeff dz

dk dz

¼0

k

¼0

u3τ κðzC þz0 Þ

¼0

Inlet

Top  exp

uðzÞ

du dz

¼

du dz

Side boundary

Symmetry plane

Outlet

Symmetry

Symmetry

du dx

¼0

dk dx dε dx

¼0

top

kðzÞ

k ¼ kexp top

Symmetry

Symmetry

εðzÞ

ε ¼ εexp

Symmetry

Symmetry

dpeff dx

dpeff dx

Symmetry

Symmetry

¼0

297

¼0

¼0

peff ¼ 0

M. Cindori et al.

Journal of Wind Engineering & Industrial Aerodynamics 173 (2018) 289–301

Fig. 11. Surfaces of the cubic building model subjected to the suburban ABL flow: a) pressure coefficient contours obtained in the simulation with WSx(z), b) pressure coefficient contours obtained in the simulation without WSx(z), c) transverse pressure coefficient distributions, d) vertical pressure coefficient distributions. CFD WSx(z) refers to the present computational results with the wind-source term, CFD WSx(z) ¼ 0 to the present computational results without the wind-source term, CR to results from Castro and Robins (1977), PA from Paterson and Apelt (1990) and L from Lim et al. (2009).

298

M. Cindori et al.

Journal of Wind Engineering & Industrial Aerodynamics 173 (2018) 289–301

(leeward) and side surfaces of the cubic building model are reported in Fig. 11 for both simulations, with and without the applied wind-source term. The computational results indicate that all of the main characteristics of the pressure distribution at the cube surfaces are captured with both computations, with and without the wind-source term. The area of the largest positive coefficient values (stagnation zone) is located near the top of the front surface, while a sudden drop occurs along with strongly negative pressure coefficients at the top surface due to flow separation. The same phenomenon occurs at the cube side surfaces. Further downstream at the side and top surfaces there is a characteristic increase in pressure coefficients due to a gradual reattachment of the flow. These results are in line with Castro and Robins (1977), Paterson and Apelt (1990) and Lim et al. (2009), Fig. 11. To preliminarily determine the impact of the wind-source term on surface pressure distribution at the cubic building model, as well as on the flow around the cube, a new scalar field ψ is defined, which is calculated as the ratio between the magnitude of the wind-source term and the pressure gradient,

computed to reduce computational time. Thus the edge of the cube has length of 0.2 m in the vertical (z) direction, 0.2 m in the horizontal (x) direction, and 0.1 m in the lateral (y) direction. The cube is placed 1 m downstream of the inlet boundary surface and 3 m upstream of the outlet boundary surface in order to prevent possible influence of the boundaries on the final solution. The geometry of the computational domain with defined boundaries is reported in Fig. 9. The basic block-structured mesh is used for spatial discretization of the computational domain. The mesh sensitivity study was performed and the presented results are mesh independent, which is in agreement with the sensitivity studies performed by Krajnovi
c and Davidson (2002) and Jureti
c (2004), as they showed that all of the relevant patterns of the flow over the cube can be captured even with coarser meshes. The mesh is generated as proposed by Franke et al. (2007) but with limited number of cells distributed in the vertical z-direction, constrained by the minimum near-wall cell center distance from the rough wall. Therefore mesh has 4 285 000 cells and it is presented in Fig. 10. The homogeneous profiles of the mean velocity, turbulent kinetic energy and dissipation rate calculated in the 2D simulation for the suburban terrain are used as the inflow condition for the 3D simulation. The von Neumann boundary condition is used at the outlet boundary surface. The aforementioned near-wall treatment proposed by Richards and Hoxey (1993) is used to establish the rough-wall conditions at the rough surface boundary with parameters defined in Table 2, while the standard wall function (Launder and Spalding, 1974) is used at surfaces of the cubic building model. The mean velocity gradient, turbulent kinetic energy and dissipation rate values are calculated using the experimental data and employed at the top boundary surface. Since the incompressible flow is assumed, the zero pressure value is defined in one cell adjacent to the outlet boundary surface. The full set of boundary conditions used in the performed 3D computations is reported in Table 7. The pressure coefficient Cp distributions at surfaces of the cubic building model are calculated from, Cp ¼

p  pref ; 1 2 ρUref 2

ψ¼

jWSx ðzÞj : 1 ρ jrpj

(37)

The ψ field is calculated in the region near the cube. The values are sampled at y-z (x ¼ 1.1 m), x-y (z ¼ 0.1 m) and x-z (y ¼ 0.1 m) planes, Fig. 12. It is noticed that the ψ field has the largest values (darkest contours) in the area far away from the cube and close to the domain boundaries, where the wind-source term dominates over the pressure gradient, which is an expected behavior. On the other hand, in the region near the cube, the local forces, such as the pressure gradient, dominate over other forces, including the wind-source term. Therefore, the ψ values are smallest in the near-cube region (white colored contours). In the x-y and x-z contour plots, larger ψ values are noticed locally. Large values of ψ near the cube are present in the separation region on the leeward side where there is a significant pressure gradient drop. In order to perform preliminary investigation on the effects of the (z), shear stress wind-source term WSx(z), the profiles of mean velocity u τzx(z) and turbulent kinetic energy k(z) are sampled at the symmetry boundary (y ¼ 1.1 m plane) far from the cubic building model. Sample line 1 is defined at the inlet boundary (x ¼ 0 m), sample line 2 at x ¼ 1.1 m downstream of the inlet boundary and sample line 3 at the

(36)

where pref and Uref are reference values obtained in the CFD simulations over suburban flat terrain in an empty 3D domain and sampled at the location of the front (windward) surface of the cubic building model at height h ¼ 0.2 m. The Cp contours recorded at the front, top (roof), rear

Fig. 12. The ψ scalar field calculated for the suburban ABL with the applied wind-source term and sampled at a) y-z (x ¼ 1.1 m), b) x-y (z ¼ 0.1 m) and c) x-z (y ¼ 0.1 m) planes. 299

M. Cindori et al.

Journal of Wind Engineering & Industrial Aerodynamics 173 (2018) 289–301

Fig. 13. Comparison between the mean velocity, shear stress and turbulent kinetic energy profiles sampled at three different positions in the domain (x ¼ 0 m – Line 1, x ¼ 1.1 m – Line 2 and x ¼ 4.2 m – Line 3) at the symmetry boundary far from the cube (y ¼ 1.1 m): a) results obtained in the ABL simulation with the applied wind-source term, b) results obtained in the ABL simulation without the applied wind-source term.

are performed for the rural, suburban and urban flat terrains and the neutral thermal stratification of the atmosphere. The results for the mean velocity, Reynolds shear stress, and turbulent kinetic energy profiles obtained in the 2D computations indicate a good agreement between the calculated profiles and the experimental data. These profiles are nearly homogeneous along the computational domain, as the mean discrepancy is 0.1% for the mean velocity profile, 0.6% for the Reynolds shear stress profile, 0.3% for the turbulent kinetic energy profile. The mean discrepancy between the experimental and computational results is 4.5% for the mean velocity profile, 2.5% for the Reynolds shear stress profile, 6% for the turbulent kinetic energy profile, when the wind-source term is applied. A more significant discrepancy of 21% is observed for the turbulent kinetic energy profile in the urban ABL. The 3D CFD simulation of the suburban ABL flow over the cubic building model was performed to assess the wind pressure distribution. It is shown that the current computational model is able to simulate the appropriate loading conditions at the cube surfaces, as good qualitative agreements with the computation without the applied wind-source term and previously reported experimental measurements and computational results are obtained. The influence of the wind-source term on the pressure field is almost negligible near the cube. The influence of the pressure gradient is negligible in the region far away from the cube where the wind-source term dominates the distribution of the ABL profiles. The wind-source term maintains the inflow profile distribution along the

outlet boundary (x ¼ 4.2 m). The results are compared to the results obtained in the simulation without the applied wind-source term and presented in Fig. 13. Fig. 13 indicates that away from the cube the wind-source term maintains the distribution of the mean velocity, Reynolds shear stress, and turbulent kinetic energy profiles throughout the computational domain. The homogeneity errors are not larger than 3%. The same homogeneity errors are achieved without the applied wind-source term, as the 2D simulation with the von Neumann boundary condition applied at the inlet and the outlet of the computational domain ensures absence of streamwise gradients. However, the absence of the wind-source term yields inability to recreate the wind-tunnel results for the boundary layer characteristics. 6. Conclusions This study presents an approach to steady Reynolds-averaged NavierStokes (RANS) simulations of the atmospheric boundary layer (ABL). The main goal is to achieve homogeneity of the relevant ABL parameters (mean velocity, Reynolds shear stress, turbulent kinetic energy) within the computational domain, which are in agreement with the respective experimental data. For that purpose, the wind-source term is implemented in the momentum equation, whereas the experimental data is used for the calibration of the computational model. The 2D simulations

300

M. Cindori et al.

Journal of Wind Engineering & Industrial Aerodynamics 173 (2018) 289–301

height of the domain and homogeneity errors are not larger than 3%. Future work will attempt to further improve the current methodology for the wind-source term calculation. Computations on wind effects on structures and the street canyons will be performed to further validate this methodology.

Kozmar, H., 2005. Scale Effects on the Structure of the Atmospheric Boundary Layer Model (In Croatian). PhD Thesis. University of Zagreb, Croatia. Kozmar, H., 2008. Influence of spacing between buildings on wind characteristics above rural and suburban areas. Wind Struct. 11 (5), 413–426. Kozmar, H., 2010. Scale effects in wind tunnel modeling of an urban atmospheric boundary layer. Theor. Appl. Climatol. 100 (1–2), 153–162. Kozmar, H., 2011a. Truncated vortex generators for part-depth wind-tunnel simulations of the atmospheric boundary layer flow. J. Wind Eng. Ind. Aerod. 99 (2–3), 130–136. Kozmar, H., 2011b. Characteristics of natural wind simulations in the TUM boundary layer wind tunnel. Theor. Appl. Climatol. 106 (1–2), 95–104. Kozmar, H., 2011c. Wind-tunnel simulations of the suburban ABL and comparison with international standards. Wind Struct. 14 (1), 15–34. Kozmar, H., 2012a. Improved experimental simulation of wind characteristics around tall buildings. J. Aerosp. Eng. 25 (4), 670–679. Kozmar, H., 2012b. Physical modeling of complex airflows developing above rural terrains. Environ. Fluid Mech. 12 (3), 209–225. Krajnovi
c, S., Davidson, L., 2002. Large-Eddy Simulation of the flow around a bluff body. AIAA J. 40 (5), 927–936. Launder, B.E., Spalding, D.B., 1974. The numerical computation of turbulent flows. Comput. Methods Appl. Mech. Eng. 3 (2), 269–289. Lim, H.C., Thomas, T.G., Castro, I.P., 2009. Flow around a cube in a turbulent boundary layer: LES and experiment. J. Wind Eng. Ind. Aerod. 97, 96–109. O'Sullivan, J.P., Archer, R.A., Flay, R.G.J., 2011. Consistent boundary conditions for flows within the atmospheric boundary layer. J. Wind Eng. Ind. Aerod. 99 (1), 65–77. Parente, A., Gorl
e, C., van Beeck, J., Benocci, C., 2011a. Improved k-ε model and wall function formulation for the RANS simulation of ABL flows. J. Wind Eng. Ind. Aerod. 99 (4), 267–278. Parente, A., Gorl
e, C., van Beeck, J., Benocci, C., 2011b. A comprehensive modelling approach for the neutral atmospheric boundary layer: consistent inflow conditions, wall function and turbulence model. Boundary-Layer Meteorol. 140, 411–428. Patankar, S.V., Spalding, D.B., 1972. A calculation procedure for heat, mass and momentum transfer in three-dimensional parabolic flows. Int. J. Heat Mass Transf. 15, 1787–1806. Paterson, D.A., Apelt, C.J., 1990. Simulation of flow past a cube in a turbulent boundary layer. J. Wind Eng. Ind. Aerod. 35, 149–176. Pontiggia, M., Derudi, M., Busini, V., Rota, R., 2009. Hazardous gas dispersion: a CFD model accounting for atmospheric stability classes. J. Hazard. Mater. 171, 739–747. Pope, S.B., 2000. Turbulent Flows. Cambridge University Press, Cambridge, UK. Richards, P.J., Hoxey, R.P., 1993. Appropriate boundary conditions for computational wind engineering models using the k-ε turbulence model. J. Wind Eng. Ind. Aerod. 46–47, 145–153. Richards, P.J., Norris, S.E., 2012. Appropriate boundary conditions for a pressure driven boundary layer. In: 18th Australasian Fluid Mechanics Conference, Launceston, Australia. Richards, P.J., Norris, S.E., 2013. Equilibrium profiles for a pressure driven boundary layer. In: 6th European and African Conference on Wind Engineering, Cambridge, UK. Richards, P.J., Norris, S.E., 2015. Appropriate boundary conditions for a pressure driven boundary layer. J. Wind Eng. Ind. Aerod. 142, 43–52. Stull, R.B., 1988. An Introduction to Boundary Layer Meteorology. Kluwer, Dordrecht, Netherlands. Wyngaard, J.C., 2010. Turbulence in the Atmosphere. Cambridge University Press, Cambridge, UK. Yaglom, A.M., 1979. Similarity laws for constant-pressure and pressure-gradient turbulent wall flows. Annu. Rev. Fluid Mech. 11, 505–540. Yang, Y., Gu, M., Chen, S., Jin, X., 2009. New inflow boundary conditions for modelling the neutral equilibrium atmospheric boundary layer in computational wind engineering. J. Wind Eng. Ind. Aerod. 97 (2), 88–95. Yang, Y., Xie, Z., Gu, M., 2017. Consistent inflow boundary conditions for modelling the neutral equilibrium atmospheric boundary layer for the SST k-ω model. Wind Struct. 24 (5), 465–480.

Acknowledgements The authors acknowledge the Croatian Science Foundation IP-201606-2017 (WESLO) support. References Balogh, M., Parente, A., Benocci, C., 2012. RANS simulation of ABL flow over complex terrains applying an Enhanced k-ε model and wall function formulation: Implementation and comparison for fluent and OpenFOAM. J. Wind Eng. Ind. Aerod. 104–106, 360–368. Balogh, M., Parente, A., 2015. Realistic boundary conditions for the simulation of atmospheric boundary layer flows using an improved k-ε model. J. Wind Eng. Ind. Aerod. 144, 183–190. Blocken, B., Stathopoulos, T., Carmeliet, J., 2007a. CFD simulation of the atmospheric boundary layer: wall function problems. Atmos. Environ. 41 (2), 238–252. Blocken, B., Carmeliet, J., Stathopoulos, T., 2007b. CFD evaluation of wind speed conditions in passages between parallel buildings — effect of wall-function roughness modifications for the atmospheric boundary layer flow. J. Wind Eng. Ind. Aerod. 95, 941–962. Boussinesq, J., 1877. Essai sur la th'eorie des eaux courantes. M
emoires pr
esent
es par divers savantsa l'Acad
emie des Sciences XXIII, pp. 1–680. Cai, X., Huo, Q., Kang, L., Song, Y., 2014. Equilibrium atmospheric boundary-layer flow: computational fluid dynamics simulation with balanced forces. Boundary-Layer Meteorol. 152 (3), 349–366. Castro, I.P., Robins, G., 1977. The flow around a surface-mounted cube in uniform and turbulent streams. J. Fluid Mech. 79 (2), 307–335. Dyrbye, C., Hansen, S., 1997. Wind Loads on Structures. John Wiley & Sons, New York. Franke, J., Hellsten, A., Schlünzen, H., Carissimo, B.E., 2007. Best Practice Guideline for the CFD Simulation of Flows in the Urban Environment. Cost action 732: quality assurance and improvement of microscale meteorological models. Garratt, J.R., 1992. The Atmospheric Boundary Layer. Cambridge University Press, Cambridge, UK. Greenshields, C.J., 2016. OpenFOAM: the OpenFOAM Foundation, User Guide Version 4.0. OpenFOAM Foundation Ltd. Hargreaves, D.M., Wright, N.G., 2007. On the use of the k-ε model in commercial CFD software to model the neutral atmospheric boundary layer. J. Wind Eng. Ind. Aerod. 95 (5), 355–369. Hellman, G., 1916. Über die Bewegung der Luft in den untersten Schichten der Atmosph€ are. Meteorol. Z. 34, 273–285. Jasak, H., 1996. Error Analysis and Estimation for the Finite Volume Method with Applications to Fluid Flow. PhD Thesis. Imperial College, University of London, UK. Jones, W.P., Launder, B.E., 1972. The prediction of laminarization with a two-equation model of turbulence. Int. J. Heat Mass Transf. 15, 301–314. Jureti
c, F., 2004. Error Analysis in Finite Volume CFD. PhD Thesis. Imperial College, University of London, UK. Jureti
c, F., Kozmar, H., 2013. Computational modeling of the neutrally stratified atmospheric boundary layer flow using the standard k-ε turbulence model. J. Wind Eng. Ind. Aerod. 115, 112–120. Jureti
c, F., Kozmar, H., 2014. Computational modeling of the atmospheric boundary layer using various two-equation turbulence models. Wind Struct. 19 (6), 687–708.

301