atmosphere

Article
Atmospheric Wind Field Modelling with OpenFOAM for
Near-Ground Gas Dispersion
Sebastian Schalau 1, *, Abdelkarim Habib 1 and Simon Michel 2

1 Bundesanstalt für Materialforschung und—Prüfung (BAM), 12205 Berlin, Germany; karim.habib@bam.de

2 Fakultät für Mathematik, Informatik und Naturwissenschaften, Meteorologisches Institut,
Universität Hamburg, 20146 Hamburg, Germany; simon.michel@uni-hamburg.de
* Correspondence: sebastian.schalau@bam.de

Abstract: CFD simulations of near-ground gas dispersion depend significantly on the accuracy of the
wind field. When simulating wind fields with conventional RANS turbulence models, the velocity
and turbulence profiles specified as inlet boundary conditions change rapidly in the approach flow
region. As a result, when hazardous materials are released, the extent of hazardous areas is calculated
based on an approach flow that differs significantly from the boundary conditions defined. To solve
this problem, a turbulence model with consistent boundary conditions was developed to ensure
a horizontally homogeneous approach flow. Instead of the logarithmic vertical velocity profile, a
power law is used to overcome the problem that with the logarithmic profile, negative velocities
would be calculated for heights within the roughness length. With this, the problem that the distance
of the wall-adjacent cell midpoint has to be higher than the roughness length is solved, so that a
high grid resolution can be ensured even in the near-ground region which is required to simulate
gas dispersion. The evaluation of the developed CFD model using the German guideline VDI



3783/9 and wind tunnel experiments with realistic obstacle configurations showed a good agreement

 between the calculated and the measured values and the ability to achieve a horizontally homogenous
Citation: Schalau, S.; Habib, A.; approach flow.
Michel, S. Atmospheric Wind Field
Modelling with OpenFOAM for Keywords: atmospheric boundary layer; OpenFOAM; gas dispersion; CFD; turbulence model;
Near-Ground Gas Dispersion. hazard assessment; horizontal homogeneity; wind field
Atmosphere 2021, 12, 933. https://
doi.org/10.3390/atmos12080933

Academic Editor: Patrick Armand 1. Introduction

When assessing hazards from industrial plant applications, calculating the gas disper-
Received: 2 July 2021
Accepted: 15 July 2021
sion is one of the main tasks for determining safety distances. Simple models for calculating
Published: 21 July 2021
the gas dispersion are generally not able to account for obstacles such as buildings or the
topography of the dispersion area. To simulate dispersion scenarios in complex areas,
Publisher’s Note: MDPI stays neutral
computational fluid dynamics (CFD) codes can be used. By solving the full Navier–Stokes
with regard to jurisdictional claims in
equations, CFD codes can simultaneously model the wind field and the dispersion of
published maps and institutional affil- pollutants, taking into account obstacles, topography, and thermal stratification, as well
iations. as the type of release (with or without momentum). Simulating gas dispersion is mainly
influenced by the wind field. The main problem when using CFD codes is the cost–benefit
ratio, as the computational time is exceedingly high compared to simpler models. To reduce
the computational time, relying on RANS (Reynolds-averaged Navier–Stokes) equations is
Copyright: © 2021 by the authors.
a common practice. To close the RANS equations, the k-ε turbulence model is used in a
Licensee MDPI, Basel, Switzerland.
wide number of publications for calculating the atmospheric wind field.
This article is an open access article
When simulating near-ground atmospheric boundary layer phenomena, fully de-
distributed under the terms and veloped wind and turbulence profiles are specified as inlet boundary conditions. Due
conditions of the Creative Commons to numerical reasons, every CFD simulation requires a computational domain of spatial
Attribution (CC BY) license (https:// extensions much above the dimensions of the area of interest. Therefore, an obstacle free
creativecommons.org/licenses/by/ approach flow region is required between the inlet and the first obstacle [1]. To ensure
4.0/). that the posed inlet boundary conditions reach the area of interest, where the release

Atmosphere 2021, 12, 933. https://doi.org/10.3390/atmos12080933 https://www.mdpi.com/journal/atmosphere

Atmosphere 2021, 12, 933 2 of 14

and dispersion occurs, a CFD model should be chosen that is able to conserve the inlet
conditions in the approach flow region, resulting in a horizontally homogenous boundary
layer flow. Otherwise, the turbulence model of the CFD code will modify the approach
flow depending on the domain length and grid resolution, resulting in different extents of
hazardous areas for the same release scenario [2].
Several publications [3–6] state that reaching a horizontally homogenous boundary
layer flow with RANS CFD models is not evident. Richards and Hoxey [7] state that a
horizontally homogenous boundary layer flow can only be reached if there is a consistency
between the turbulence model and the boundary conditions, such as the atmospheric inlet
condition and wall functions. They propose inlet conditions on the basis of a logarithmic
wind profile, as well as necessary modifications to the k-ε turbulence model. Several other
authors [8–10] present comparable solutions for the problem of horizontal homogeneity.
Whilst all aforementioned publications try to achieve a solution by modifying the constants
of the k-ε model, [11] they achieved the horizontal homogeneity by introducing an addi-
tional source term to the dissipation equation [12] and by introducing it in the turbulent
kinetic energy equation.
To reduce the computational cost, wall functions are used for describing the near wall
area using a roughness factor, to avoid modelling even the smallest details (e.g., vegetation,
fences, cars). Different types of roughness parameters and their suitability for modelling
horizontally homogenous boundary layer flows were discussed [13–15]. However, for all it
was stated that the general formulations of wall functions for rough walls (based on the
equivalent sand grain roughness k s or the roughness length z0 ) introduce a dependency of
the wall-adjacent cell size from the defined roughness value. This results in wall-adjacent
cell sizes that are much too big for simulating near-ground or near-obstacle effects with
acceptable accuracy. If a near-wall cell size smaller than the requirements due to the defined
roughness is chosen, the simulation will not result in a horizontally homogenous boundary
layer flow [13].
For near-ground gas dispersion scenarios, a coarse wall-adjacent mesh will lead to an
unrealistic prediction of the gas cloud sizes, as the concentration gradients near ground
cannot be resolved with sufficient accuracy.
In this work, a CFD model based on the open source CFD code OpenFOAM is
presented, which was developed with the aim to compute a horizontally homogenous
boundary layer flow and to overcome the wall resolution restrictions due to wall roughness.
In this model the wind profile is described by a power law instead of a log wind profile. The
used flow solver is a transient, compressible flow solver, to also be able to simulate in the
future gas dispersion, gas cloud explosions, or gas cloud combustion. In this work we will
focus on the wind field modelling. The model will be evaluated systematically using the
German Guideline VDI 3783/9 [16,17]. This guideline provides test cases in which generic
obstacle configurations are investigated and objective evaluation criteria are provided. To
evaluate the model performance for a more realistic obstacle configuration comparable to
an industrial area, wind tunnel measurements are used. In opposite to free field trials, they
show a high reproducibility, and the average boundary conditions can be kept constant over
a long time period to produce a statistically sound database. The wind tunnel data used
were generated at the University of Hamburg in the context of an actual research project.
The quality of the presented model to simulate a horizontally homogenous boundary layer
approach flow and the flow around obstacles will be demonstrated and evaluated with the
VDI Guideline 3783/9 and the aforementioned wind tunnel experiments.

2. Turbulence Model and Boundary Conditions for Neutral Stratification

The presented CFD model is based on the transient and compressible solver rhoRe-
actingBuoyantFoam [18,19], as well as on the k-e-turbulence model of Launder and Spald-
ing [20] and El Tahry [21], implemented in OpenFOAM 5.0. The differential equations
Atmosphere 2021, 12, 933 3 of 14

describing the turbulent kinetic energy k (1) and the dissipation ε (2) are coupled to the
Navier–Stokes equations through the turbulent viscosity νt (3):
→ →
       
∂(ρk ) νt 2
= ∇· ρ + ν ∇k − ∇ · ρkU + Pk − ρ ∇ · U k − ρe , (1)
∂t σk 3

→ → e2
       
∂(ρe) νt C1 e 2
= ∇· ρ + ν ∇e − ∇ · ρeU + Pk − C1 ρ ∇ · U e − C2 ρ + Se , (2)
∂t σe k 3 k

k2
νt = Cµ . (3)
e
The values of the model constants C1 , C2 , Cµ , σk and σε in Equations (1)–(3) are based
→
on [22]. The variable ρ is the density of the fluid and U is the velocity vector, Pk describes
the production of the turbulent kinetic energy and Sε is an additional source term, which
will be derived in the following.
With the idealized assumption of a horizontally homogenous boundary layer flow
and with the additional assumptions that:
- the flow is stationary,
- the flow fields are free of horizontal gradients,
- the vertical velocity is zero,
Equations (1) and (2) can be simplified to:
 
∂ νt ∂k
0= ρ + Pk − ρe (4)
∂z σk ∂z

e2
 
∂ νt ∂e C e
0= ρ + 1 Pk − C2 ρ + Se (5)
∂z σe ∂z k k
The dynamic viscosity is assumed to be negligible (νt  ν).
Assuming a two-dimensional horizontally homogenous boundary layer flow, the
turbulence production is defined as
 2
∂Ux
Pk = ρνt (6)
∂z

Based on Equations (4) to (6), models have been developed to achieve a horizontally
homogenous boundary layer solution (e.g., [7,10,12]). All these models are based on a
logarithmic wind velocity profile. This assumption leads to a dependency of the grid
resolution near ground from the roughness length [13]. To overcome this restriction, in this
work a power-law wind profile is used as inlet boundary condition
!m
z
Ux = Ure f (7)
zre f

The value of the exponent m of the wind profile is set in accordance to the roughness
length z0 , which is now only indirectly linked to the governing model equations, allowing
us to choose the grid resolution independently of the ground roughness.
The vertical profile of the turbulent kinetic energy of the incident flow, as described
by, e.g., Richards and Hoxey [7], is dependent on the wall shear stress velocity which can
be defined as nearly constant with height for a neutrally layered atmosphere

u2
k = p∗ (8)
Cµ
Atmosphere 2021, 12, 933 4 of 14

Ure f κ
u∗ =  z +z  (9)
re f 0
ln z0

To achieve a horizontally homogenous atmospheric boundary layer flow, it is manda-

tory that the flow profiles at the inlet boundary satisfy the equations of the turbulence
model. This is achieved for the turbulent kinetic energy Equation (4) with the condition
that ∂k/∂z = 0 for
!
ρe = Pk (10)
From Equations (3), (6)–(8), and (10), the vertical profile of the turbulent dissipation ε
at the inlet boundary is   m
z
u2∗ · m · Ure f zre f
e= (11)
z
The inlet boundary conditions (7), (8), and (11) only satisfy the equation for the
turbulent dissipation (5) with

νt ∂2 e C1 e ∂Ux 2 e2
 
∂ρ νt ∂e ρ ∂νt ∂e
Se = − − −ρ 2
− ρνt + C2 ρ (12)
∂z σe ∂z σe ∂z ∂z σe ∂z k ∂z k

With this modification, horizontally homogenous boundary layer flows can be simu-
lated for near-ground grid resolutions that are not dependent on the roughness length. This
source term is evaluated solely based on the vertical inlet profiles at the beginning of the
simulation and is therefore time independent. As this source term is relevant and designed
to achieve a horizontally homogenous flow without obstacles, this term is only active in
the upstream region when simulating obstacles in the domain. In the vicinity of obstacles,
the turbulence model (Equations (1)–(3)) is used with Sε = 0. The term Sε compensates
the expected dissipation of turbulent kinetic energy with the standard k--model in an
obstacle-free flow domain.
Besides the described modifications to the turbulence model, the simulation of a
horizontally homogenous boundary layer flow also requires a formulation for wall func-
tions consistent to the modifications made. Wall functions available in OpenFOAM de-
scribe the turbulent energy dissipation and production for hydraulically smooth walls
(epsilonWallFunction), as well as the turbulent viscosity based on the roughness length
(nutkAtmRoughWallFunction). To avoid reintroducing here the dependency of the possible
grid resolution on the roughness length, and at the same time to satisfy Equation (11), the
turbulent dissipation at the wall-adjacent cell center is defined as
 m
zp
u2∗ · m · ure f zre f
ep = (13)
zp

Coupling of the wall function with the flow field occurs in analogy to the OpenFOAM
epsilonWallFunction [18], by replacing the friction velocity by the wall adjacent turbulent
kinetic energy through Equation (8)
 m
p zp
Cµ · k p · m · ure f zre f
ep = (14)
zp

In accordance with Equation (6), the turbulence production in the wall adjacent cell is
 2
∂u x ∂u x ∂u x
Pk,p = µt = µt (15)
∂z ∂z ∂z

In analogy to the OpenFOAM nutkAtmRoughWallFunction [18], one of the two

velocity gradients of the wall function is replaced by the analytical derivative of the power
Atmosphere 2021, 12, 933 5 of 14

law profile and the other by its spatial discretization. Thus, the turbulence production can
be rewritten as  m
zp

∂u x
2 m · u re f zre f ∆u x
Pk,p = µt ≈ µt (16)
∂z zp ∆z
The wall sheer stress τ0 , considered as constant within the Prandtl layer, is

∂u x
τo = ρu2∗ = µt (17)
∂z
and by integration with respect to the height and assuming a no slip wall, Equation (17)
becomes
up
ρu2∗ = µt (18)
zp
resulting in the wall function for the turbulent viscosity, by substituting up with Equation (7):

u2∗ z p
νt =  m (19)
zp
ure f zre f

u∗ = Cµ0.25 k0.5 (20)

For the turbulent kinetic energy, a zero gradient is specified at the wall, in accordance
with the inlet boundary condition.
Studies [13,23] observed that choosing fixed shear stress as a boundary condition on
the top of the domain instead of a symmetry boundary condition leads to a less inhomo-
geneous horizontal boundary layer. This was also confirmed by [24] who noticed a better
agreement with experimental values when setting a fixed shear stress as the top boundary
condition. Therefore, in the presented model, constant values for the wind speed and the
turbulent values based on Equations (7)–(9) are set as the top boundary conditions. All
other values are set to zero gradient.

3. Results
Several investigations have been carried out to assess the accuracy of the presented
model. Beginning with a study of the horizontal homogeneity of the boundary layer
flow when using the original and the modified code over a benchmark with the VDI
3783/9 [16,17], to a comparison with wind tunnel experiments for complex terrain, the
quality of the new model is shown.

3.1. Horizontally Homogeneous Boundary Layer

By using a power law to describe the vertical velocity profile and introducing the
additional source term Sε , it should be possible to simulate a horizontally homogenous
boundary layer flow with a grid resolution not depending on the roughness length. Figure 1
shows the wind velocity in the main direction Ux calculated by the k-e turbulence model
without modification (standard model) and with the model presented here (modified
model) at a height of 1 m above the ground from the inlet up to a distance of 200 m. A
power law with a reference velocity of 3 m/s at a reference height of 10 m was specified
as the inlet boundary condition and a roughness length of 1.2 m was assumed. For the
standard model, a wall function based on the roughness length is chosen (OpenFOAM
nutkAtmRoughWallFunction), whilst the roughness length is approximated by the ex-
ponent of the power law in the modified model. The simulations were carried out in a
400 m × 8 m × 100 m (l × w × h) pseudo 2D domain without obstacles. In Figure 1, it can
be clearly seen that the standard model, as expected, does not conserve the wind profile
over the distance, whereas the modified model shows insignificant deviations from the
given inlet velocity.
Atmosphere 2021, 12, 933 6 of 14

Figure 1. Development of the boundary layer for the standard model and the model presented.

To ensure the presented model being able to carry out simulations of a horizontally
homogenous boundary layer flow where the grid resolution can be chosen independently
of the surface roughness, which is not possible for standard turbulence models in CFD [13],
further simulations were carried out in the same domain. In these simulations, the wind
speed at 10 m in height over the ground was varied in a range of 1 m/s up to 5 m/s, and
for wind profiles corresponding to roughness lengths between 0.1 m and 1.2 m, covering
a typical range for hazard assessment scenarios. The wind profiles were investigated
after a distance of 300 m, a typical distance from the inlet according to the best practice
guidelines [1], where distances from the inlet of two to eight times the obstacle height are
recommended as the approach flow distance. The grid resolution ranged from 0.5 m down
to 0.125 m for the wall-adjacent cell, corresponding to cell midpoint heights z p of 0.25 m
down to 0.0625 m.
Figure 2 shows the relative deviation between the calculated vertical wind profile in
the main direction (Ux ) and the inlet boundary condition. As no significant changes were
observed over the range of wind speeds, the figure shows exemplarily the results for a
wind speed of 3 m/s at 10 m over the ground. The wind profiles after the investigated
distance of 300 m do not differ from the inlet profile by more than +/−5% close to the
ground. At a height of 1 m or more from the ground, the deviation between the inlet profile
and the computed values is negligible. Similar observations can be made for the eddy
viscosity profiles (Figure 3). The maximum deviation reached is around +/−10%, which
seems small enough to be qualified as acceptable.
As these observations for the wind velocity and eddy viscosity profiles are valid for
any of the investigated combinations of roughness length and grid resolution, it can be
concluded that the modified turbulence model and wall functions presented here are able
to simulate a horizontally homogenous boundary layer flow for near-wall grid resolutions
that can be chosen independently to the roughness length.

3.2. VDI 3783/9 Benchmark

The German Guideline VDI 3783/9 [16,17] is used to evaluate the quality of wind
field models for a built-up environment. In the guideline, a number of test cases are
defined to evaluate different characteristics of the wind field model. In some test cases,
comparisons of the model have to be performed with wind tunnel experiments that were
carried out for generic obstacle configurations (scale 1:200) by the meteorological institute
of the University of Hamburg [25]. As the aim of this work is to provide a wind field
model for, e.g., dispersion scenarios in built up environment, the evaluation of the model
presented here is carried out with these test cases.
Atmosphere 2021, 12, 933 7 of 14

Figure 2. Wind profiles at 300 m distance from the inlet for varying surface roughness and grid resolution.

Figure 3. Eddy viscosity profiles at 300 m distance from the inlet for varying surface roughness and
grid resolution.

The evaluation method of the guideline VDI 3783/9 gives a quantitative estimate of
the model quality by determining a hit ratio q

n ∑ ni
q= = (21)
N N
Atmosphere 2021, 12, 933 8 of 14

It indicates in percent the proportion of the total of correctly predicted values n in the
total number of comparison values N. For successful validation in comparison with wind
tunnel experiments, it is required that q > 66%. The number of correctly predicted values ni
results from the comparison of the normalized results Pi and the normalized comparison
values Oi [17]: (   )
 −Oi 
1, i f  Pi O ≤ D or | Pi − Oi | ≤ W
ni = (22)
i

0, otherwise
In Equation (22), W is the absolute permitted difference and D the permitted difference,
which are defined separately for each case in [16,17].
The focus is laid on the test cases c3, c4 and c6 of the 2005 version of the guideline [16],
as a simple generic case for a multiple building scenario (test case c6) is available. Figure 4
shows the geometrical representation of the test cases, and the corresponding parameters
are given in Table 1. Parameters of each test case. For the test cases c3 and c4 (single
obstacle), the guideline requires evaluating the wind field in the near field (grey dotted
region in Figure 4a,b), as well as in the whole computational domain. For test case c6, only
the whole domain has to be considered.

Figure 4. Geometrical representation of the test cases (a) test case c3 cube with 90◦ approach flow,
(b) test case c4 cube with 270◦ approach flow, (c) test case c6 7 × 3 array of buildings.

Table 1. Parameters of each test case.

Test Case C3 Test Case C4 Test Case C6

Domain in m (l × w × h) 500 × 200 × 200 500 × 200 × 200 520 × 400 × 200
Nr. of Buildings 1 Building 1 Building 3 × 7 Buildings
Building Size in m (l × w × h) 25 × 25 × 25 25 × 25 × 25 30 × 20 × 25
Mesh size 1.60 Mio cells 1.80 Mio cells 7.30 Mio cells
z0 = 0.10 m z0 = 0.10 m z0 = 0.10 m
m = 0.22 m = 0.21 m = 0.21
Atmospheric Inlet Conditions Ure f = 5.04 m/s Ure f = 4.43 m/s Ure f = 6.30 m/s
zre f = 75 m zre f = 75 m zre f = 132 m
u∗ = 0.31 m/s u∗ = 0.27 m/s u∗ = 0.36 m/s
upwind of the building: 25 m
downwind of the building: 50 m
Definition of near fields -
each side transverse to the direction of flow: 12.5 m
above the building: 12.5 m
Permitted difference D 25% 25% 25%
Absolute permitted
0.06 0.06 0.10
difference W
Atmosphere 2021, 12, 933 9 of 14

During a grid sensitivity study, the grid was refined from near-wall mesh sizes of
1.0 m down to 0.25 m, showing that the hit ratio is not very sensitive to the grid refinement,
as only negligible differences for all refinements can be observed. In the following, only
results for a mesh with a near-wall mesh size of 0.5 m will be discussed, although the
observations made are also valid for all investigated meshes.
Table 2 Hit ratio in % for all three test cases. shows the hit ratio q for all three test cases,
subdivided into hit ratios for qUx —corresponding to the main wind direction, qUy —the
crosswind direction, and qUz —the vertical to the main wind direction. For all test cases
and all three wind direction components, the required minimum hit ratio of 66% is reached
or exceeded when evaluating the whole computational domain.

Table 2. Hit ratio in % for all three test cases.

qUx qUy qUz

Whole Whole Whole
Near Field Near Field Near Field
Domain Domain Domain
Test case c3 83 72 92 88 91 85
Test case c4 81 70 69 62 77 61
Test case c6 78 - 87 - 80 -

While the near-field values for test case c3 meet the requirements of the guideline,
test case c4 does not reach the required minimum value for the hit ratio for some velocity
components. Due to the high velocity gradients occurring at obstacles where a flow
detachment or recirculation can be observed, the choice of the turbulence model has a
significant influence on the quality of the result, whilst, in general, a higher deviation
between computed values and experimental data is to be expected in these cases. Although
the minimum criterion for the hit ratio is reached for test case c3, in test case c4 with a
45◦ approach flow, this can only be achieved for the hit ratio of the main wind direction
component. Study [26] stated that the approach flow in the experiments for test case c4
was not exactly 45◦ but showed a slight deviation of 2◦ . It is still to be examined whether
correcting the simulation values for the 2◦ deviation might lead to better hit ratios.

3.3. Wind Tunnel Experiments

After having compared the model results to such generic cases as the ones defined
in [16,17], the model is still to be validated for a more realistic obstacle configuration.
Due to the high reproducibility and to a lack of field experiments, wind tunnel data
generated at the University of Hamburg in the context of an actual research project are
used. The experiments were carried out at the EWTL large boundary layer wind tunnel
of the University of Hamburg in a scale of 1:100. The 25 m long wind tunnel has an 18 m
long, 4 m wide and approx. 3 m high test section. The investigated model area (Figure 5)
was planned with aerial photos and building data from real industrial areas. In real scale,
it covers about 300 m × 300 m. The model consists of 20 cuboid shaped buildings which
have varying side lengths of 20 m to 40 m. Each building has a height of 15 m. The
spaces between the buildings range between 10 m and 80 m. In the middle of the model
area, a large building is located. This main building has a base area of 30 m × 60 m. The
main building has three openings on the sides and one opening on the roof with a size of
4 m × 4 m. In the future, it is planned to carry out gas dispersion experiments with this
setup where the mentioned openings on the buildings will be used as gas inlets.
Flow measurements were performed with Laser-Doppler Anemometry (LDA). The
LDA measures two components of the turbulent wind vector simultaneously with a spatial
resolution of approx. 110 mm in full scale. To receive statistical representative time series
in a turbulent flow, measurements are performed for several minutes. In real scale, this
represents a measurement taken for 30 h under stationary meteorological conditions.
Atmosphere 2021, 12, 933 10 of 14

Figure 5. Wind tunnel and obstacle configuration.

Two different boundary layer approach flows are modelled at a scale of 1:100 to
investigate the influence of the flow on the measurements inside the model area. One
boundary layer flow represents the conditions of a “moderately rough” boundary layer.
This type of boundary layer flow develops over grasslands or farmlands. The second
wind boundary layer represents a flow that develops over suburban area and is described
by the definition of a “rough” boundary layer. The roughness classes are defined by
VDI 3783/12 [27]. The similarity between wind tunnel experiments and real scale has
been ensured by checking that the turbulence intensity, the wind profiles, the spectral
distribution of turbulent energy and the lateral homogeneity are similar to nature-given
profiles and distributions. In addition, the functional relationship between the roughness
length z0 and the wind profile exponent m has been checked, according to [28].
In the following, the simulated wind fields are compared to the profiles measured in
the wind tunnel for the “moderately rough” case. For the considered case, the roughness
length for the “moderately rough” case is z0 = 0.06 m, corresponding to an exponent of
the wind profile power law of m = 0.17. The same approach flow (Figure 6a) was used
for two different wind directions, case M1 and case M2 (with respect to the obstacles,
Figure 6b). The horizontal flow components (Ux , Uy ) were measured around the obstacles
at 327 positions for the direction M1 and at 383 positions for M2, at a height over the
ground varying from 2 m up to 20 m. Figure 6a shows a good agreement of the measured
and calculated normalized vertical wind velocity profile of the approach flow. For each
horizontal, homogeneity was verified.

Figure 6. (a) Comparison of the boundary layer flow measured in the wind tunnel and calculated
with the presented model; (b) schematic representation of the obstacle configuration and the approach
flows M1 and M2.
Atmosphere 2021, 12, 933 11 of 14

For the following investigations of the model performance, an objective evaluation

method comparable to the VDI Guideline 3783/9 is desirable. In the VDI Guideline 3783/9,
in the 2017 version [17], a test case is defined for a real built-up area. The evaluation
criteria for this case (D = 25% and W = 0.08, threshold for q = 66%) are applied to the cases
investigated here, as they also represent realistic topographies. It should be noted that the
evaluation criteria of the VDI Guideline are defined as case-specific and therefore might
not be fully adapted to the cases M1 and M2. Nevertheless, as the permitted differences
(D and W) and the hit ratio q do not change significantly for all cases of the guideline, this
approach appears to be permissible.
For the cases M1 and M2, grid refinements were carried out to investigate if the hit
ratio might be independent of the resolution, as this was the case in Section 3.2. For both
cases, it is noticed that, with each refinement of the grid, the hit ratio increased up to a
point where a grid resolution of 0.25 m adjacent to the obstacles finally led to a successful
“validation” in terms of reaching the hit ratio threshold. These simulations already showed
a total number of cells > 20 Mio. As further refinement would have led to domains with at
least 80 Mio. cells, no further refinement was investigated due to the resulting extremely
long simulation time.
Figures 7 and 8 show the repartition of the measuring points at heights of 2 m
(Figures 7a and 8a) and 8 m over the ground (Figures 7b and 8b), respectively, for cases M1
and M2. Whilst in the VDI Guideline the experimental data used to define the thresholds
are evenly spread over the whole domain, the measuring points in cases M1 and M2 are
concentrated around buildings and close to the ground, making it harder to reach the
threshold. Both figures represent the investigated area in such way that the approach
flow is oriented parallel to the x-axis coming from the left. Full black bullets show a point
where the calculated hit ratios qUx and qUy reach values higher than the threshold of 66%,
half-filled bullets show points where the hit ratio of one velocity component did not reach
the threshold, and for the blank bullets, no component of the horizontal velocity reached
the threshold.

Figure 7. Illustration of the hit ratio for case M1 (a) at 2 m height over ground, (b) at 8 m height
over ground.

With increasing height over the ground, the number of points which show a hit ratio
higher than the threshold increases in both cases. This is due to the fact that the closer to
the ground, the more complex the turbulence-influenced effects are on the flow. To resolve
these effects in a satisfactory manner, at least a LES (Large-Eddy Simulation) would be
required. Table 3 Hit ratio for the whole domain of case M1 and M2. shows that the hit
ratio for the whole domain exceeds the threshold value of 66% for both cases. For the aim
of providing an adequate wind field for hazard assessment purposes, the reached accuracy
is therefore sufficient.
Atmosphere 2021, 12, 933 12 of 14

Figure 8. Illustration of the hit ratio for case M2 (a) at 2 m height over ground, (b) at 8 m height
over ground.

Table 3. Hit ratio for the whole domain of case M1 and M2.

qUx qUy
M1 71 89
M2 67 82

4. Discussion
The transient, compressible flow solver rhoReactingBuoyantFoam of OpenFOAM 5.0
in combination with the standard k-ε turbulence model is promising for simulating the
dispersion process of pollutants in the atmosphere. Due to inconsistencies between the
required boundary conditions for modelling atmospheric boundary layer flow and the
formulation of the turbulence model, simulating a horizontally homogenous boundary
layer flow is not possible.
As the horizontally homogenous approach flow is mandatory for reliable estimations
of the safety distances for hazard assessment purposes, when a pollutant is released near
ground in the atmosphere, a modification to the turbulence model has been made and
was presented in this work. Introducing an additional source term in the equation for the
turbulent dissipation and replacing the log law profile by a power law profile for the wind
field, in combination with a new formulation of the wall functions, showed the ability to
simulate a horizontally homogenous boundary layer flow. Additionally, the grid resolution
adjacent to the walls can now be chosen independently of the roughness length defined
on the ground. For gas dispersion purposes in a built-up domain, the grid resolution in
the vicinity of buildings is crucial to achieve a satisfactory resolution of the concentration
gradients, leading to a reliable and reproducible prediction of the safety distances. The
model’s performance for this application is currently under investigation.
The evaluation of the presented model occurred in three steps. In the first step, it was
checked against the original model formulation for its performance to produce horizontally
homogenous wind and turbulence profiles, and showed a clear advantage compared to
the original formulation by achieving the homogeneity. In the second step, the model was
evaluated considering the rules defined by the VDI Guideline 3783/9 [16] where generic
obstacle configurations are presented, and an objective statistical method for evaluation
using a so-called hit ratio is defined. The model was evaluated against three test cases of
the guideline with a single cube in a 90◦ approach flow and a 45◦ approach flow, as well
as for an array of 7 × 3 obstacles. Considering the whole domain, the model performed
well in all cases and reached hit ratios clearly fulfilling the defined threshold. For the two
single cube cases, an evaluation of the “near-field” had to be carried out. The hit ratio
showed lower values than for the whole domain, due to restrictions of the k-ε turbulence
model in the vicinity of obstacles. The third evaluation step consisted of comparing the
Atmosphere 2021, 12, 933 13 of 14

presented model to wind tunnel measurements and using the evaluation criteria of the
VDI Guideline 3783/9 [17] to obtain an objective measure of the performance, even though
the criteria might not fully match the considered case. Nevertheless, the presented model
showed again a good performance by reaching the required thresholds.

Author Contributions: Conceptualization, S.S. and A.H.; model development and validation, S.S.;
wind tunnel measurements, S.M.; All authors were actively involved in writing and revising the
manuscript. All authors have read and agreed to the published version of the manuscript.
Funding: The IGF-Project No.: 20719 BG of the Research Association Society for Chemical Engi-
neering and Biotechnology (DECHEMA), Theodor-Heuss-Allee 25, 60486 Frankfurt am Main, was
funded by the German Federation of Industrial Research Associations (AiF) within the framework
of the Industrial Collective Research (IGF) support program by the Federal Ministry for Economic
Affairs and Energy due to a decision of the German Bundestag.
Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.
Conflicts of Interest: The authors declare no conflict of interest.

References
1. Franke, J.; Hellsten, A.; Schlünzen, H.; Carissimo, B. Best Practice Guideline for the CFD Simulation of Flows in the Urban Environment:
COST Action 732 Quality Assurance and Improvement of Microscale Meteorological Models; Meteorological Institute University of
Hamburg: Hamburg, Germany, 2007.
2. Batt, R.; Gant, S.E.; Lacome, J.-M.; Truchot, B. Modelling of stably-stratified atmospheric boundary layers with commercial CFD
software for use in risk assessment. Chem. Engineer. Trans. 2016, 48, 61–66.
3. Richards, P.; Younis, B. Comments on “prediction of the wind-generated pressure distribution around buildings” by E.H. Mathews.
J. Wind. Eng. Ind. Aerodyn. 1990, 34, 107–110. [CrossRef]
4. Quinn, A.D.; Wilson, M.; Reynolds, A.M.; Couling, S.B.; Hoxey, R.P. Modelling the dispersion of aerial pollutants from agricul-
tural buildings—An evaluation of computational fluid dynamics (CFD). Comput. Electron. Agric. 2001, 30, 219–235. [CrossRef]
5. Zhang, C.X. Numerical predictions of turbulent recirculating flows with a k–ε model. J. Wind. Eng. Ind. Aerod. 1994, 51, 177–201.
[CrossRef]
6. Riddle, A.; Carruthers, D.; Sharpe, A.; McHugh, C.; Stocker, J. Comparisons between FLUENT and ADMS for atmospheric
dispersion modelling. Atmos. Environ. 2004, 38, 1029–1038. [CrossRef]
7. Richards, P.J.; Hoxey, R.P. Appropriate boundary conditions for computational wind engineering models using the k-ε mod-el. J.
Wind. Eng. Ind. Aerod. 1993, 46–47, 145–153. [CrossRef]
8. Alinot, C.; Masson, C. k–ε model for the atmospheric boundary layer under various thermal stratifications. J. Sol. Energy Eng.
2005, 127, 438–443. [CrossRef]
9. Yang, W.; Quan, Y.; Jin, X.; Tamura, Y.; Gu, M. Influences of equilibrium atmosphere boundary layer and turbulence param-eter
on wind loads of low-rise buildings. J. Wind. Eng. Ind. Aerodyn. 2008, 96, 2080–2092. [CrossRef]
10. Yang, Y.; Gu, M.; Chen, S.; Jin, X. New inflow boundary conditions for modelling the neutral equilibrium atmospheric boundary
layer in computational wind engineering. J. Wind. Eng. Ind. Aerodyn. 2009, 97, 88–95. [CrossRef]
11. Pontiggia, M.; Derudi, M.; Busini, V.; Rota, R. Hazardous gas dispersion: A CFD model accounting for atmospheric stability
classes. J. Hazard. Mater. 2009, 171, 739–747. [CrossRef] [PubMed]
12. Van Der Laan, M.P.; Kelly, M.C.; Sørensen, N.N. A new k-epsilon model consistent with Monin-Obukhov similarity theory. Wind.
Energy 2016, 20, 479–489. [CrossRef]
13. Blocken, B.; Stathopoulos, T.; Carmeliet, J. CFD simulation of the atmospheric boundary layer: Wall function problems. Atmos.
Environ. 2007, 41, 238–252. [CrossRef]
14. Zhang, X. CFD Simulation of Neutral ABL Flows; Danmarks Tekniske Universitet, Risø Nationallaboratoriet for Bæredygtig Energi:
Kongens Lyngby, Denmark, 2009.
15. Parente, A.; Gorlé, C.; van Beeck, J.; Benocci, C. Improved k-ε model and wall function formulation for the RANS simula-tion of
ABL flows. J. Wind. Eng. Ind. Aerodyn. 2011, 34, 107–110.
16. Environmental Meteorology—Prognostic Microscale Windfield Models—Evaluation for Flow around Buildings and Obstacles; VDI
Guideline 3783, Part 9; Beuth Verlag: Berlin, Germany, 2005.
17. Environmental Meteorology—Prognostic Microscale Windfield Models—Evaluation for Flow around Buildings and Obstacles; VDI
Guideline 3783, Part 9; Beuth Verlag: Berlin, Germany, 2017.
18. OpenFOAM. User Guide. Version 7. OpenFOAM—The OpenFOAM Foundation. 2019. Available online: https://openfoam.org
(accessed on 24 July 2020).
19. Weller, H.G.; Tabor, G.; Jasak, H.; Fureby, C. A tensorial approach to computational continuum mechanics using ob-ject-oriented
techniques. Comput. Phys. 1998, 12, 620–631. [CrossRef]
Atmosphere 2021, 12, 933 14 of 14

20. Launder, B.E.; Spalding, D.B. The numerical computation of turbulent flows. Comput. Method Appl. Mech. Eng. 1974, 3, 269–289.
[CrossRef]
21. El Tahry, S.H. k-epsilon equation for compressible reciprocating engine flows. J. Energy 1983, 7, 345–353. [CrossRef]
22. Jones, W.; Launder, B. The prediction of laminarization with a two-equation model of turbulence. Int. J. Heat Mass Transf. 1972, 15,
301–314. [CrossRef]
23. Hargreaves, D.; Wright, N. On the use of the k–ε model in commercial CFD software to model the neutral atmospheric boundary
layer. J. Wind. Eng. Ind. Aerodyn. 2007, 95, 355–369. [CrossRef]
24. Franke, J.; Sturm, M.; Kalmbach, C. Validation of OpenFOAM 1.6.x with the German VDI guideline for obstacle resolving
micro-scale models. J. Wind. Eng. Ind. Aerodyn. 2012, 104–106, 350–359. [CrossRef]
25. Leitl, B. Compilation of Experimental Data for Validation of Micro Scale Dispersion Models (CEDVAL). Available online:
https://mi-pub.cen.uni-hamburg.de/index.php?id=433 (accessed on 1 May 2021).
26. Eichhorn, J. MISKAM Version 5.02—Evaluierung nach VDI-RL 3783/9. Available online: http://www.lohmeyer.de/ (accessed
on 1 May 2021).
27. Environmental Meteorology—Physical Modelling of Flow and Dispersion Processes in the Atmospheric Boundary Layer—Application of
Wind Tunnels; VDI Guideline 3783, Part 12; Beuth Verlag: Düsselsdorf, Germany, 2000.
28. Counihan, J. Adiabatic atmospheric boundary layers: A review and analysis of data from the period 1880–1972. Atmos. Environ.
1975, 9, 871–905. [CrossRef]