
Numerical Wind Potential Analysis in Urban Environments

A thesis accepted by the Faculty of Aerospace Engineering and Geodesy of
the University of Stuttgart in partial fulfilment of the requirements for the

degree of Doctor of Engineering Sciences (Dr.-Ing.)

by

Maximilian Siegfried Ulrich von der Grün

born in Rosenheim

Main referee: Prof. Dr.-Ing. Ewald Krämer
Co-referee: Prof. Dr. Ursula Voß

Date of defense: 29.09.2023

Institute of Aerodynamics and Gas Dynamics
University of Stuttgart

2024


I dedicate this thesis to my family and my friends,

none of this would have been possible without their love, patience and support.


Acknowledgements

During my PhD candidature I had the pleasure to be supported by numerous friends,

colleagues, students and my family. At this point I would like to thank these people for

making this work possible.

I wish to express my gratitude to my supervisors, Prof. Dr. Ursula Voß and Prof. Dr.

Ewald Krämer, for their incredible amount they have taught me during the course of this

thesis, their professional guidance and assistance.

My thanks also go to Pradip Zamre and Dr. Thorsten Lutz for the excellent and inspiring

discussions and they always came up with new ideas.

It was a great honor being part of the interdisciplinary Joint Graduate Research Training

Group ”Windy Cities” in corporation with the University of Stuttgart, the Hochschule

für Technik Stuttgart and the Hochschule Esslingen, meeting and working together with

people from all over the world.

I would like to thank the Stuttgart Wind Energy @ Insitute of Aircraft Design (SWE),

University of Stuttgart, in particular Yiyin Chen for conducting the Lidar measurement

campaign and providing the wind measurement data.

The computational resources were provided by the bwUniCluster funded by the Ministry

of Science, Research and Arts and the universities of Baden-Württemberg, Germany,

within the framework program bwHPC.


Abstract

Due to the global warming, the man-made CO2 emissions have to be reduced. A possible

role can play the switch of fossil energies towards electric energy. The increasing trend

to urbanism rises the energy demand in urban areas. Short power supply distances in

urban areas, where the most energy is needed, reduce the total energy costs and the power

loss in long distance grids. A local renewable power generation, such as by small wind

turbines, which are mounted in urban areas, can help to overcome these problems and

meet the increasing demand in electric energy.

It is not substantially investigated how urban areas can be used economically for small

wind turbines. The influences of environmental and atmospheric parameters, such as

atmospheric turbulence intensities or vegetational properties, on the flow field and the

energy yield of small wind turbines are still unknown. Usually, the evaluation of suitable

positions for small wind turbines is based on wind speed. For that measurement data

from metmasts or sonic measurement devices are used which do not sufficiently consider

local wind effects. But the measurement data is just valid within a short range since

local effects, for example induced by buildings, complex terrain structures or vegetation,

can strongly influence the wind field and are not represented in calculation methods of

the AEP. Furthermore, an extensive three-dimensional validation of urban wind field

simulations with on-site measurement data is still lacking. So far only punctual sonic

data are used to validate simulation data.

The objective of this work is the numerical investigation of the urban wind field at the

university campus Morgenstelle in Tübingen (southern Germany) which is chosen as a

test site for this study. This is done using CFD and highly resolved DES simulations with

ANSYS Fluent in this work. The simulations include complex terrain and vegetation

with the focus on the energy yield of small wind turbines. The building geometries used

in the simulations are based on a real 3D city model geometry. A workflow is used

to prepare the city models for urban CFD simulations from geometry optimisation to

meshing requirements. An already existing vegetation model is extended to consider local

tree heights. A sensitivity analysis is performed to study the influence of different meshing

parameters and domain sizes. Additionally, the effect of the vegetation in general and the

effects of various environmental parameters such as the atmospheric turbulence intensity,

different tree heights and foliage densities to simulate the wind flow in summer and winter

times are investigated. The numerical setup and the vegetation modelling are compared


with on-site LiDAR measurement data in Tübingen. A new method is presented how to

use planar on-site wind measurement data as an inflow boundary condition.

Since the AEP is one of the most important key parameters defining suitable locations

for wind turbines, a procedure for calculating local AEP values and local wind statistics

is presented. For that the wind flow is simulated from four evenly distributed wind

directions with and without buildings which surround the buildings of interest. For that,

a new building model is presented which models the buildings by means of a volume force

and an indicator function. That allows a faster preparation and a simplified meshing of

complex building geometries. The calculation of the AEP is based on large scale synthetic

wind statistics which have been already available. The simulations of the four wind

directions are used in the method to achieve local frequencies of wind directions and wind

speeds. The AEP is calculated for each building using a real horizontal and vertical axis

wind turbine. An already existing approach for correcting the AEP by the atmospheric

turbulence intensity and the angle of attack of the wind turbines is implemented in the

procedure to obtain a more realistic and precise prediction of the energy yield. Finally,

suitable locations for small wind turbines are evaluated and discussed.

Especially, the local wind statistics and the approach for correcting the AEP by the

atmospheric turbulence intensity lead to a significant difference in the calculated AEP

values. The study further shows a potential of small wind turbines at the university

campus Morgenstelle in Tübingen, e.g., the highest AEP value is reached on the highest

building on the campus when a horizontal axis wind turbine is used.


Kurzzusammenfassung

Aufgrund der globalen Erwärmung müssen die menschengemachten Emissionen reduziert

werden. Eine mögliche Rolle kann dabei der Umstieg von fossilen Energien auf elektrische

Energie spielen. Der zunehmende Trend zur Urbanisierung erhöht den Energieverbrauch

in urbanen Gegenden. Kurze Stromversorgungswege in städtischen Gebieten, wo die

meiste Energie benötigt wird, reduzieren die Gesamtenergiekosten und den Leistungsver-

lust in Fernnetzen. Eine dezentrale regenerative Stromerzeugung, beispielsweise durch

Kleinwindkraftanlagen, die in städtischen Gebieten aufgestellt werden, kann helfen, diese

Probleme zu überwinden und den steigenden Bedarf an elektrischer Energie zu decken.

Es ist bislang nicht ausreichend untersucht, inwiefern städtische Gebiete für Kleinwind-

kraftanlagen wirtschaftlich genutzt werden können. Die Einflüsse von Umweltparame-

tern und atmosphärischen Parametern, wie z.B. atmosphärische Turbulenzintensitäten

oder Vegetationsparameter, auf das Strömungsfeld und den Energieertrag von Kleinwind-

kraftanlagen sind bislang noch unbekannt. Normalerweise basiert die Auswertung von

geeigneten Standorten für Kleinwindkraftanlagen auf der Windgeschwindigkeit. Dafür

werden Messdaten von Messmasten oder Anomometerdaten verwendet, welche lokale

Windeffekte nicht ausreichend berücksichtigen. Aber die Messdaten gelten nur innerhalb

eines kleinen Bereichs, weil lokale Effekte, z.B. induziert durch Gebäude, komplexe Bo-

denstrukturen oder Vegetation, das Windfeld stark beeinflussen können und nicht in den

Berechnungsmethoden des Jahresenergieertrags abgebildet sind. Darüber hinaus fehlt

noch eine umfassende dreidimensionale Validierung urbaner Windfeldsimulationen mit

Messdaten vor Ort. Bisher werden nur punktuelle Anemometerdaten zur Validierung von

Simulationsdaten verwendet.

Ziel dieser Arbeit ist die numerische Untersuchung des städtischen Windfeldes auf dem

Universitätscampus Morgenstelle in Tübingen (Süddeutschland), der als Testgelände für

diese Studie ausgewählt wurde. Das soll in dieser Arbeit mit Hilfe von CFD und hoch

aufgelösten DES Simulationen mit ANSYS Fluent geschehen. Die Simulationen um-

fassen komplexes Terrain und Vegetation mit dem Fokus auf den Energieertrag von Klein-

windkraftanlagen. Die Gebäudegeometrien, die in den Simulationen verwendet werden,

basieren auf einer realen 3D-Stadtmodellgeometrie. Es wird ein Workflow verwendet, um

die Stadtmodelle für urbane CFD-Simulationen vorzubereiten, von der Geometrieopti-

mierung bis hin zu den Anforderungen an die Vernetzung. Ein bereits bestehendes Ve-

getationsmodell wird erweitert, um lokale Baumhöhen zu berücksichtigen. Es wird eine


Sensitivitätsanalyse durchgeführt, um den Einfluss verschiedener Vernetzungsparameter

und Domaingrößen zu untersuchen. Zusätzlich werden die Auswirkungen der Vegetation

im Allgemeinen und die Auswirkungen verschiedener Umweltparameter wie die atmo-

sphärische Turbulenzintensität, verschiedene Baumhöhen und Laubdichten zur Simulation

der Windströmung im Sommer und imWinter untersucht. Der numerische Aufbau und die

Vegetationsmodellierung werden mit LiDAR-Messdaten vor Ort in Tübingen verglichen.

Es wird ein neues Verfahren vorgestellt, wie man flächige Vor-Ort-Windmessdaten als

Einlassrandbedingung verwenden kann.

Da der Jahresenergieertrag einer der wichtigsten Schlüsselparameter zur Bestimmung

geeigneter Standorte für Windkraftanlagen ist, wird ein Verfahren zur Berechnung lokaler

Jahresenergieerträge und lokaler Windstatistiken vorgestellt. Dafür wird die Windströ-

mung aus vier gleichmäßig verteilten Windrichtungen mit und ohne Gebäude simuliert,

welche die Gebäude, die von Interesse sind, umgeben. Dazu wird ein neues Gebäude-

modell vorgestellt, das die Gebäude mittels Volumenkraft und einer Indikatorfunktion

modelliert. Dies ermöglicht eine schnellere Vorbereitung und eine vereinfachte Vernet-

zung von komplexen Gebäudegeometrien.

Die Berechnung des Jahresenergieertrags basiert auf großskaligen synthetischen Wind-

statistiken, die bereits verfügbar sind. Die Simulationen der vier Windrichtungen werden

im Verfahren verwendet, um lokale Häufigkeiten von Windrichtungen und Windgeschwin-

digkeiten zu erhalten. Der Jahresenergieertrag wird für jedes Gebäude mit einer realen

Horizontal- und Vertikalachsenwindkraftanlage berechnet. Ein bereits vorhandener Ansatz

zur Korrektur des Jahresenergieertrags durch die atmosphärische Turbulenzintensität und

den Anstellwinkel der Windturbinen wird in das Verfahren implementiert, um eine rea-

listischere und genauere Vorhersage des Energieertrags zu erhalten. Abschließend werden

geeignete Standorte für Kleinwindkraftanlagen diskutiert und bewertet.

Insbesondere die lokale Windstatistik und der Ansatz zur Korrektur des Jahresenergieer-

trags um die atmosphärische Turbulenzintensität führen zu einem signifikanten Unter-

schied in den berechneten Jahresenergieertragswerten. Die Studie zeigt weiterhin ein

Potential von Kleinwindkraftanlagen auf dem Universitätscampus Morgenstelle in Tübin-

gen, z.B. der höchste Jahresenergieertrag auf dem höchsten Gebäude auf dem Campus

erreicht, wenn eine Horizontalachsenwindkraftanlage verwendet wird.


Contents

List of Symbols IV

List of Subscripts XI

List of Abbreviations XII

List of Figures XV

List of Tables XXI

1 Introduction 1

1.1 Previous studies about urban wind field simulations . . . . . . . . . . . . . 2

1.2 Objective . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9

2 Fundamental Equations and Numerical Approach 11

2.1 Governing equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11

2.2 Methodology of simulations . . . . . . . . . . . . . . . . . . . . . . . . . . 12

2.3 Generation and injection of atmospheric turbulence . . . . . . . . . . . . . 15

3 Testcase Validation: Flow around a High-rise Building Structure 18

3.1 Description of the experiment . . . . . . . . . . . . . . . . . . . . . . . . . 18

3.2 Numerical setup . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19

3.3 Results and validation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21

3.3.1 Turbulence model and shielding function . . . . . . . . . . . . . . . 22

3.3.2 Turbulence generation method . . . . . . . . . . . . . . . . . . . . . 25

3.4 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26

4 Geometry Preparation and Meshing for Urban Simulations 28

4.1 Site information . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29

4.2 Building preparation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30

4.3 Terrain modelling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32

4.4 Vegetation modelling details . . . . . . . . . . . . . . . . . . . . . . . . . . 32

4.5 Modelling of surrounding buildings . . . . . . . . . . . . . . . . . . . . . . 36

4.6 Meshing of the urban domain . . . . . . . . . . . . . . . . . . . . . . . . . 40

I


5 Evaluation: Sensitivity Analysis 42

5.1 Mesh convergence study . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42

5.1.1 Simulation setup . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42

5.1.2 Influence of meshing parameters . . . . . . . . . . . . . . . . . . . . 43

5.2 Domain size analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47

5.2.1 Simulation setup . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47

5.2.2 Domain height . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48

5.2.3 Domain width . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52

5.3 Influence of environmental conditions . . . . . . . . . . . . . . . . . . . . . 54

5.3.1 Simulation setup . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54

5.3.2 Investigation of different atmospheric turbulence intensities . . . . . 57

5.3.3 Investigation of different vegetation modelling . . . . . . . . . . . . 61

5.3.4 Summer - winter comparison . . . . . . . . . . . . . . . . . . . . . . 63

5.3.5 Investigation of different tree heights . . . . . . . . . . . . . . . . . 67

5.4 Simulation of surrounding buildings . . . . . . . . . . . . . . . . . . . . . . 71

5.4.1 Site information . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71

5.4.2 Simulation setup . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72

5.4.3 Wake flow analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . 73

5.5 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78

6 Evaluation against LiDAR Measurement Data 80

6.1 Preparation and filtering of LiDAR measurement data . . . . . . . . . . . 80

6.2 Numerical setup and inflow conditions . . . . . . . . . . . . . . . . . . . . 82

6.3 Comparison with on-site LiDAR data . . . . . . . . . . . . . . . . . . . . . 83

6.4 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85

7 Windfield Evaluation and local Annual Energy Production 87

7.1 Local site information . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87

7.2 Numerical setup and inflow conditions . . . . . . . . . . . . . . . . . . . . 88

7.3 Flow analysis of the simulated wind directions . . . . . . . . . . . . . . . . 89

7.4 Procedure for calculating the local AEP . . . . . . . . . . . . . . . . . . . 96

7.4.1 Interpolation of the synthetic wind statistics onto the grid . . . . . 97

7.4.2 Calculation of local wind statistics . . . . . . . . . . . . . . . . . . 99

7.4.3 Calculation of the local AEP . . . . . . . . . . . . . . . . . . . . . . 101

7.4.4 Correction of the local AEP . . . . . . . . . . . . . . . . . . . . . . 102

7.5 Local wind statistics and local AEP . . . . . . . . . . . . . . . . . . . . . . 104

7.6 Influence of surrounding buildings on local AEP . . . . . . . . . . . . . . . 110

7.7 Evaluation of locations for small wind turbines . . . . . . . . . . . . . . . . 112

7.8 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 118

8 Conclusion and Outlook 120

II


Bibliography 124

A Appendix 137

A.1 Data sheet of horizontal axis wind turbine [5] . . . . . . . . . . . . . . . . 137

A.2 Data sheet of vertical axis wind turbine [6] . . . . . . . . . . . . . . . . . . 138

III


List of Symbols

Symbol Description Unit

Ain Area of inflow plane [m2]

Acs Cross-sectional area of the domain [m2]

BWS Wind speed bins in synthetic wind statistics
[m
s
]

Bsim Wind speed bins in simulations
[m
s
]

CDES Constants in the IDDES model [−]

CDES1 Constants in the IDDES model [−]

CDES2 Constants in the IDDES model [−]

Cw Constant in equ. 2.17 [−]

Cµ Constant in the k-production term [−]

Co Courant number [−]

Dw Diameter of wind turbine [m]

F1 Blending function in the k-ω SST model [−]

F2 Blending function in the k-ω SST model [−]

Fb,i Building drag force
[
kg
ms2

]

Fw,i Forest drag force
[

kg
m2s2

]

G (x, x′) LES filter operation [−]

Hb Building height [m]

Hc Height of cube [m]

Ht Tree height [m]

Ht,av Average tree height [m]

IV


Hw Height of wind turbine [m]

I Turbulence intensity [−]

L Turbulent length scale [m]

L11 Longitudinal integral length scale [m]

LAI Leaf area index [−]

Lb,x Length of building in x-direction [m]

LD Domain length [m]

Ld Downwind length of the domain [m]

Lu Upwind length of the domain [m]

Lγ Constant in equ. 7.10 [m]

M(i) Measurement signal [−]

N Number of vortex points [−]

Nil Number of inflation layers [−]

P Power of wind turbine [W ]

P10%TI Wind turbine power with 10% turbulence intensity [W]

Pw Wind power [W]

Pcorr Corrected wind turbine power [W]

PK k production term in equ. 2.14
[
kg
s3·m

]

P̃K k production term in equ. 2.7
[
kg
s3·m

]

Pl,i Wind turbine power according to the

lower wind speed bin limit [W]

Pmax Rated power of wind turbine [W]

PsimTI Wind turbine power with simulated I [W]

Pu,i Wind turbine power according to the

upper wind speed bin limit [W]

Re Reynolds number [−]

V


Reλ Taylor-scale Reynolds number [−]

Ri Richardson number [−]

S Strain rate magnitude
[
1
s
]

Sij Strain rate tensor
[
1
s
]

Vc Cell volume [m3]

a Local foliage density
[

1
m
]

a1 Constant in equ. 2.16 [−]

c Threshold value in median filter [−]

cd Drag coefficient [−]

cp Pressure coefficient [−]

cp,b Pressure coefficient for buildings [−]

cr Resistance coefficient [−]

cVM Constant in the Vortex Method [−]

d Side length of cube [m]

dw Wall distance [m]

f̃d Blending function in equ. 2.17 [−]

fe Elevating function in equ. 2.17 [−]

fG10% Gaussian distribution function (TI= 10%) [−]

fGsim Gaussian distribution function (TI from the simulations) [−]

g Gravitational acceleration
[m
s2
]

g(x) Value in median filter [−]

g′(x) Median value in median filter [−]

hmax Maximum edge length of a cell [m]

k Turbulent kinetic energy
[
m2

s2
]

l Filtered length scale [m]

lg Local grid scale [m]

VI


lEI Demarcation length scale between the energy-

containing range of eddies (l > lEI) and

smaller eddies (l < lEI) [m]

lIDDES IDDES length scale [m]

lLES LES length scale [m]

lRANS RANS length scale [m]

p Pressure [Pa]

p∞ Atmospheric pressure [Pa]

q Induction factor of HAWT [−]

r Radius in building model [m]

s Skewness of a numerical cell [−]

t Time [s]

u Wind speed
[m
s
]

u′ Fluctuating velocity field
[m
s
]

ui Velocity component (ui = ūi + u′i)
[m
s
]

u′i Fluctuating velocity component
[m
s
]

u′∗i Fluctuating velocity in the Vortex Method
[m
s
]

ūi Time-averaged velocity component
[m
s
]

uref Reference velocity in velocity profile
[m
s
]

ūsim Time-averaged and wind-direction averaged

simulated wind speed
[m
s
]

ūWS Yearly-averaged wind speed in synthetic wind

statistics
[m
s
]

u60 Wind speed from 60 ◦ wind direction
[m
s
]

u90 Wind speed from 90 ◦ wind direction
[m
s
]

u150 Wind speed from 150 ◦ wind direction
[m
s
]

VII


w Raster width in building model [m]

wm Weighting factor in median filter [−]

w1, w2 Weightings in the bilinear interpolation method [−]

x x-coordinate of the domain [−]

y y-coordinate of the domain [−]

y+ Normalized wall distance [−]

z z-coordinate of the domain [−]

zhom Demarcation height between lateral homogeneous

velocity profile (z > zhom) and lateral

inhomogeneous velocity profile (z < zhom) [m]

zref Reference height in velocity profile [m]

Γ Circulation
[
m2

s

]

∆ Cell size [m]

∆be Cell size of building edges [m]

∆bs Cell size of building surfaces [m]

∆t Cell size terrain surfaces [m]

Θ Virtual potential temperature [K]

Φsim Frequency based on the simulations [−]

ΦWS Frequency in the synthetic wind statistics [−]

Ψ Load factor [−]

α Exponent in velocity profile [−]

αω Constant in ω-equation 2.9 [−]

β Constant in ω-equation 2.9 [−]

β∗ Constant in k-equation 2.7 [−]

βγ, βγ0 Constants in equ. 7.10 [−]

γ Angle of attack [◦]

VIII


δ Angle of local wind direction [◦]

ϵ Dissipation rate
[
m2

s3
]

ζ Blocking ratio [−]

η Spatial distribution of vortex points [−]

θe Angle for an equiangular face or cell [◦]

θmax Largest angle in a face or cell [◦]

θmin Smallest angle in a face or cell [◦]

κ Wavenumber
[

1
m
]

κc Filter cutoff wavenumber
[

1
m
]

λ2 Second eigenvalue of
(
S2 +Ω2

) [
1
m
]

µ Dynamic viscosity
[
kg
m s

]

µt Modelled turbulent viscosity
[
kg
m s

]

ν Kinematic viscosity
[
m2

s

]

ξ Vortex size [m]

ρ Density
[
kg
m3

]

σ Standard deviation [−]

σk Constant in k-equation 2.7 [−]

σω, σω2 Constants in ω-equation 2.9 [−]

τij Shear stress tensor
[

kg
m s2

]

τ sgsij Subgrid stress tensor
[

kg
m s2

]

ϕ Scalar quantity [−]

ϕ̄ Filtered quantity [−]

ψ Ratio of ūsim to ūWS [−]

ω Specific dissipation rate
[
1
s
]

C⃗c Cell coordinates [m]

R⃗1−4 Coordinates of synthetic wind statistics [m]

IX


n⃗ Normal vector of the gradient of the streamwise

mean velocity in equ. 2.21 [m]

n⃗nord Vector of the northern wind direction [m]

n⃗x Vector in x-direction [m]

u⃗ Velocity vector
[m
s
]

u⃗h Horizontal velocity vector
[m
s
]

u⃗xz Velocity vector in x-z-plane
[m
s
]

v⃗′ Planar fluctuating velocity field
[m
s
]

x⃗ Location vector [m]

z⃗ Unit vector in equ. 2.21 [-]

S Symmetric tensor of the velocity gradient
[
1
s
]

Ω Asymmetric tensor of the velocity gradient
[
1
s
]

X


List of Subscripts

Subscript Description

av average

max maximum

rms root mean square

x x-coordinate

y y-coordinate

z z-coordinate

∞ freestream

0 origin

XI


List of Abbreviations

Abreviation Description

ABL Atmospheric Boundary Layer

AEP Annual Energy Production

CFD Computational Fluid Dynamics

DDES Delayed Detached Eddy Simulation

DES Detached Eddy Simulation

DNS Direct Numerical Simulation

HAWT Horizontal Axis Wind Turbine

IDDES Improved Delayed Detached Eddy Simulation

LAI Leaf Area Index

LES Large Eddy Simulation

LGL Landesamt für Geoinformation und Landentwicklung

(State Agency for Spatial Information and Rural Development)

LiDAR Light Detection And Ranging

LKM Linear Kinematic Model

LoD Level of Detail

PANS Partially-Averaged Navier-Stokes

pdf probability density function

RANS Reynolds-Averaged Navier-Stokes

SDES Shielded Detached Eddy Simulation

XII


sgs subgrid stress

SpS Spectral Synthesizer

SST Shear Stress Tensor

TI Turbulence Intensity

TIN Triangulated Irregular Network

URANS Unsteady Reynolds-Averaged Navier-Stokes

VAWT Vertical Axis Wind Turbine

VM Vortex Method

XIII


List of Figures

3.1 Computational domain of the validation case [104] (modified). . . . . . . . 20

3.2 Time-averaged first- and second-order velocity profiles in x-direction and

elements of the Reynolds stress tensor . . . . . . . . . . . . . . . . . . . . . 23

3.3 Iso-surface of λ2 for the DDES shielding function. . . . . . . . . . . . . . . 24

3.4 Iso-surface of λ2 for the IDDES shielding function. . . . . . . . . . . . . . . 24

3.5 Time-averaged first- and second-order velocity profiles in x-direction and

elements of the Reynolds stress tensor . . . . . . . . . . . . . . . . . . . . . 26

4.1 Model of the campus Morgenstelle of the University of Tübingen (Ger-

many). Relative height (z) of the complex terrain with all buildings included. 30

4.2 Optimized building geometries (blue) and original geometries (green). . . . 31

4.3 Profiles of a over the dimensionless tree height Ht for different LAI, based

on [115]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33

4.4 Seasonal LAI-values for pedunculate oaks in different vegetation times [23]

(modified). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33

4.5 Forested zones (dark green) in the computational model. . . . . . . . . . . 34

4.6 Local dimensionless tree heights in the domain. . . . . . . . . . . . . . . . 34

4.7 Probability density functions and tree height distributions for σ/Ht,av =

0.05, 0.1, 0.15 and the local forested zone in Tübingen. . . . . . . . . . . . . 35

4.8 Influence of the local tree height standard deviation σ on |ū| for different
distances (z/Ht,av) to the ground. . . . . . . . . . . . . . . . . . . . . . . . 35

4.9 Influence of the local tree height standard deviation σ on k for different

distances (z/Ht,av) to the ground. . . . . . . . . . . . . . . . . . . . . . . . 35

4.10 Schematic view of a mesh to check, if the cell is inside the building: A red

cross must be captured by the green circle with radius r. . . . . . . . . . . 37

4.11 Time-averaged streamwise velocity with different cr. . . . . . . . . . . . . . 38

4.12 Time-averaged streamwise rms velocity ux,rms with different cr. . . . . . . . 38

4.13 Streamwise velocity with different raster widths w and corresponding r. . . 39

4.14 Velocity cross section: the flow infiltrates the cube at the black arrows

because r is chosen too small. . . . . . . . . . . . . . . . . . . . . . . . . . 39

4.15 Meshed geometry with buildings, terrain and forested zones (light areas). . 40

XV


5.1 Computational domain of the simulations used in the mesh convergence

study. Relative height (z) of the complex terrain with the buildings. . . . . 43

5.2 Iso-surfaces of the λ2 criterion, obtained from the simulations with different

cell sizes ∆ and coloured by the velocity magnitude from blue (0m/s) to red

(11m/s)). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44

5.3 cp,b values averaged over all buildings for different meshing parameters. . . 45

5.4 Computational domain used in the domain height study. The vertical eval-

uation lines are highlighted in red. Relative height (z) of the complex

terrain with all buildings included. . . . . . . . . . . . . . . . . . . . . . . 48

5.5 Blockage ratio ζ of the computational domain with a height of 6Hb,max,

12Hb,max and 18Hb,max over the domain length LD. . . . . . . . . . . . . . 49

5.6 Computational domain used in the domain height study with the horizontal

evaluation line which is highlighted in blue. Largest cross-sectional area of

the domain Acs,max is indicated by the dashed line. . . . . . . . . . . . . . 50

5.7 Time-averaged velocity magnitude |ū| evaluated at the horizontal line 4

indicated as blue line in fig. 5.6 which is perpendicular to the main wind

direction. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51

5.8 Time-averaged velocity magnitude |ū| evaluated at vertical lines 1 - 3 using

different domain heights, wherein line 1 is above building A and line 2 is

above building B. Blue: 6Hb,max, red: 12Hb,max, green: 18Hb,max. . . . . . 51

5.9 Computational domain with a width of 10Hb,max. The vertical evaluation

lines are highlighted in red. Relative height (z) of the complex terrain with

all buildings included. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52

5.10 Computational domain with the horizontal evaluation lines highlighted in

blue. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53

5.11 Time-averaged velocity magnitude evaluated at vertical lines 1 - 3 using

different domain widths, wherein line 1 is above building A, line 2 is above

building B and line 3 is above building C. Blue: 10Hb,max and red: 16Hb,max 53

5.12 Time-averaged velocity magnitude evaluated at horizontal lines in y-direction

as shown by the blue lines in fig. 5.9 in a height of 70m. Blue: 10Hb,max

and red: 16Hb,max . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54

5.13 Computational domain for the simulations of different environmental con-

ditions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55

5.14 Vertical evaluation lines in front of building B at y=80m (blue) and in

front of building A at y=-4.5m (red). . . . . . . . . . . . . . . . . . . . . . 57

5.15 Time-averaged velocity magnitude evaluated in front of building B at y=80m

with an inlet turbulence intensity of 10% (green), 20% (blue) and 30% (red). 58

5.16 Turbulent kinetic energy evaluated in front of building B at y=80m with

an inlet turbulence intensity of 10% (green), 20% (blue) and 30% (red). . 59

XVI


5.17 Time-averaged velocity magnitude evaluated in front of building A at y=-

4.5m with an inlet turbulence intensity of 10% (green), 20% (blue) and

30% (red). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59

5.18 Turbulent kinetic energy evaluated in front of building A at y=-4.5m with

an inlet turbulence intensity of 10% (green), 20% (blue) and 30% (red). . 60

5.19 Time-averaged streamwise velocities ux plotted over height with different

types of forests: volume force modelled forest (black), solid body modelled

forest (green) and without forest (red), based on [50]. . . . . . . . . . . . . 62

5.20 Time averaged velocity magnitude of the wind field with summer foliation

in the forested zone (LAI=3.5). . . . . . . . . . . . . . . . . . . . . . . . . 63

5.21 Time averaged velocity magnitude of the wind field with winter foliation

in the forested zone (LAI=0.5). . . . . . . . . . . . . . . . . . . . . . . . . 63

5.22 Time-averaged velocity magnitude evaluated in front of building B at y=80m

with summer foliation (blue) and winter foliation (red). . . . . . . . . . . . 64

5.23 Turbulent kinetic energy evaluated in front of building B at y=80m with

summer foliation (blue) and winter foliation (red). . . . . . . . . . . . . . . 65

5.24 Time-averaged velocity magnitude evaluated in front of building A at y=-

4.5m with summer foliation (blue) and winter foliation (red). . . . . . . . . 65

5.25 Turbulent kinetic energy evaluated in front of building A at y=-4.5m with

summer foliation (blue) and winter foliation (red). . . . . . . . . . . . . . . 66

5.26 Time-averaged velocity magnitude of the wind field with 50% of the tree

height. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67

5.27 Time-averaged velocity magnitude of the wind field with 150% of the tree

height. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67

5.28 Time-averaged velocity magnitude in front of building B at y=80m with

tree heights of 50% (blue), 100% (green) and 150% (red) related to the

original tree heights. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69

5.29 Turbulent kinetic energy evaluated in front of building B at y=80m with

tree heights of 50% (blue), 100% (green) and 150% (red) related to the

original tree heights. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69

5.30 Time-averaged velocity magnitude in front of building A at y=-4.5m with

tree heights of 50% (blue), 100% (green) and 150% (red) related to the

original tree heights. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70

5.31 Turbulent kinetic energy evaluated in front of building A at y=-4.5m with

tree heights of 50% (blue), 100% (green) and 150% (red) related to the

original tree heights. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70

5.32 University campus Tübingen with modelled surrounding buildings (left)

and surrounding buildings with solid geometries (right). . . . . . . . . . . . 72

XVII


5.33 Time-averaged velocity magnitude with a simulated wind direction of 240 ◦,

evaluated at z=15m with surrounding buildings included (left) and without

surrounding buildings (right). . . . . . . . . . . . . . . . . . . . . . . . . . 74

5.34 Time-averaged velocity magnitude with a simulated wind direction of 330 ◦,

evaluated at z=15m with surrounding buildings included (left) and without

surrounding buildings (right). . . . . . . . . . . . . . . . . . . . . . . . . . 75

5.35 Time-averaged velocity magnitude on building A with surrounding build-

ings (red) and without surrounding buildings (blue). . . . . . . . . . . . . . 77

5.36 Time-averaged velocity magnitude on building B with surrounding build-

ings (red) and without surrounding buildings (blue). . . . . . . . . . . . . . 78

6.1 Velocity in streamwise direction u over time - unfiltered and filtered by the

original median filter, von der Grün et al. [50]. . . . . . . . . . . . . . . . . 81

6.2 Velocity in streamwise direction u over time - unfiltered and filtered by the

extended version of the median filter, von der Grün et al. [50]. . . . . . . . 81

6.3 Height(z) of the complex terrain with forested zones, buildings and instan-

taneous streamwise velocity ux from the simulation, von der Grün et al.

[50]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83

6.4 Simulated time-averaged velocity magnitude (red) compared with LiDAR

data (black), see also [50]. . . . . . . . . . . . . . . . . . . . . . . . . . . . 84

6.5 Section of the domain showing the locations of the evaluated lines, see also

[50]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84

7.1 Computational domain with forested zones and buildings used for the calcu-

lation of the annual energy production. The four simulated wind directions

are illustrated by the arrows. . . . . . . . . . . . . . . . . . . . . . . . . . . 88

7.2 Time-averaged velocity magnitude (left) and turbulent kinetic energy (right)

with a simulated wind direction of 60 ◦ and evaluated at z=15m. . . . . . . 89

7.3 Time-averaged velocity magnitude (left) and turbulent kinetic energy (right)

with a simulated wind direction of 150 ◦ and evaluated at z=15m. . . . . . 91

7.4 Time-averaged velocity magnitude (left) and turbulent kinetic energy (right)

with a simulated wind direction of 240 ◦ and evaluated at z=15m. . . . . . 93

7.5 Time-averaged velocity magnitude (left) and turbulent kinetic energy (right)

with a simulated wind direction of 330 ◦ and evaluated at z=15m. . . . . . 95

7.6 One of the 14 used synthetic large scale wind statistics with frequencies

(%) for wind direction sectors (◦) and wind speeds (m/s) at the university

campus in Tübingen [3]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98

7.7 Large scale synthetic wind statistics at the university campus in Tübingen

[3] (dashed lines) and the computational domain (black lines). . . . . . . . 99

7.8 Pressure coefficient for the HAWT and the VAWT in dependence of γ. . . 103

7.9 Power loss of the HAWT with increasing γ. . . . . . . . . . . . . . . . . . . 103

XVIII


7.10 HAWT power curve with a pdf with I=20% and 70%. . . . . . . . . . . . 104

7.11 HAWT power curve affected by turbulence intensity. . . . . . . . . . . . . 104

7.12 Frequencies Φ of local wind directions based on the simulations and syn-

thetic wind statistics on building 12 with surrounding buildings in the

simulations. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105

7.13 Frequencies Φ of local wind speeds based on the simulations and synthetic

wind statistics on building 12 with surrounding buildings in the simulations.106

7.14 Calculated AEP associated to buildings at the university campus Morgen-

stelle in Tübingen for VAWT. . . . . . . . . . . . . . . . . . . . . . . . . . 107

7.15 Calculated AEP associated to buildings at the university campus Morgen-

stelle in Tübingen for HAWT. . . . . . . . . . . . . . . . . . . . . . . . . . 108

7.16 Calculated AEP associated to buildings at the university campus Morgen-

stelle in Tübingen for VAWT when the surrounding buildings are included

in the simulations. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111

7.17 Calculated AEP associated to buildings at the university campus Morgen-

stelle in Tübingen for HAWT when the surrounding buildings are included

in the simulations. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112

7.18 Local wind statistics including frequencies of wind directions and wind

speeds compared to synthetic wind statistics for building 12. . . . . . . . . 113

7.19 Local wind statistics including frequencies of wind directions and wind

speeds compared to synthetic wind statistics for building 13. . . . . . . . . 114

7.20 Local wind statistics including frequencies of wind directions and wind

speeds compared to synthetic wind statistics for building 19. . . . . . . . . 115

7.21 Local wind statistics including frequencies of wind directions and wind

speeds compared to synthetic wind statistics for building 21. . . . . . . . . 117

XIX


List of Tables

5.1 Mesh properties of the simulations with various cell sizes. . . . . . . . . . . 44

5.2 Mesh properties of the simulations with various cell resolutions of the ter-

rain and building surfaces and building edges. . . . . . . . . . . . . . . . . 46

5.3 Mesh properties of the simulations with various numbers of inflation layers. 46

7.1 Comparison of the horizontal and vertical wind turbines used in the calcu-

lation of the local AEP. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101

7.2 Comparison of the corrected and uncorrected AEP values using the VAWT. 109

XXI


1 Introduction

The human-induced climate change is one of the most important issues the humanity has

to face. Compared to the preindustrial time the global temperature is increased by 1.2 ◦C

until now and the annual rise even accelerated. To limit the increase in temperature to

1.5 ◦C, the CO2 emissions have to be reduced by 50% by 2030 compared to 1990 and

eliminated by 2050 [99]. In the Paris Agreement all nations commonly agreed to limit

the global warming well below 2 ◦C. An increase in temperature of 2 ◦C could be already

a threshold value which leads to irreversible changes in the ecosystem, stronger increase

in global warming than anticipated and impacts on human societies and economies [124].

The European Union set the objective to limit global warming to 1.5 ◦C in his climate

legislation. That means that 40 percent of the emissions have to be reduced by 2030

compared to 1990 and eliminated by 2050 which is transferred in national laws.

The world population increased from 2015 to 2020 by 78 million people per year and 9.7

billion people are expected to live on the planet by 2050, projected by the UN Depart-

ment of Economics and Social Affairs [8]. Worldwide 56.2 percent of the population live

in urban areas, in the US this ratio already increased from 65% in 1950 to 83% until

now. This trend will continue in the future, so the design and the supply of urban areas

become more important. Due to the increase in urban population and industry, the ur-

ban energy consumption will rise in future. The climate legislation makes it necessary to

switch the primary energy to electric energy. Short power supply distances in urban areas,

where the most energy is needed, reduce the total energy costs and the power loss in long

distance grids. The additional lack of people’s acceptance delays the implementations of

small wind turbines. A local renewable power generation, such as by small wind turbines,

which are mounted in urban areas, can reduce the total energy costs and the power loss

in long distance grids and meet the increasing demand in electric energy.

There is a presumed potential for small wind turbines in cities. For example, in Munich

the technical potential for small wind turbines for buildings with flat roofs and a height

above 20m is estimated to be 27.000MWh/year corresponding to 6750 households [4]. To

estimate the yield potential, a geo-information system based approach was used, which

combines wind data, exemplary system characteristics and building-specific characteristics

such as roof shape and height. The wind data were taken from the German Wheather Ser-

vice based on test reference years, while the building-specific parameters are taken from a

3D model of the city of Munich. Based on the performance curve of a wind turbine, yields

1


were calculated at various hub heights for each square kilometre of the test reference years

and assigned to the respective buildings or roof areas. In Surakarta, the yearly available

energy of buildings with a height above 50m is assumed to be 1221 kWh/m2/year [127].

In this study, the assessment of the wind energy potential was analyzed using a Weibull

distribution in a period from 2011 to 2015 on the height of 50 m. The wind speed data

were taken from a meteorological station nearby and extrapolated to the intended location

and height. In a next step, the average wind velocity and the standard deviation of the

extrapolating result of the wind data were calculated.

Varying wind directions are not considered when calculating the annual energy yield.

Another weakness of the studies regarding the wind potential in Munich and Surakarta

is that local flow effects around single buildings and the vegetation are not taken into

account. However, the latter are considered to be of great importance for an accurate

prediction of the annual energy yield and the potential for small wind turbines [41].

This work is part of the interdisciplinary Joint Graduate Research Training Group ”Windy

Cities” which investigate the economic use of small wind turbines in urban areas. That

also includes the research about new web-based visualisation, new energy storage tech-

nologies and intelligent load management of small wind turbines [7].

1.1 Previous studies about urban wind field

simulations

Computational Fluid Dynamics (CFD) is a widely used and appropriate tool for wind field

simulations with increasing significance [24]. For instance, the aerodynamics and loads of

wind turbine blades [45], the wake [101], the performance [119] or the layout of wind farms

[67] can be studied with CFD. In the last years CFD has established as a common tool in

architecture and urban planning [24, 133, 136, 128]. Blocken [24] gave an overview of the

development of CFD for computational wind engineering and pointed out the advantage

of CFD compared to on-site measurements. CFD provides detailed information on the

relevant flow variables in the whole computational domain, especially when using scale

resolving simulations. But CFD still needs high-quality measurements for solution veri-

fication and validation studies. In addition, the results of CFD simulations can be very

sensitive to the computational parameters. Valger et al. [136] provide studies, in which

influences of turbulence models, thermal stratification regimes or concentration of gaseous

emissions are investigated with CFD. Moreover, CFD is used to optimize building-roof

shapes for the wind energy exploitation and to simulate and optimize real geometries of

vertical axis wind turbines on building roofs [128].

An urban wind field is mainly affected by buildings, vegetation and urban terrain [41]. A

study using generic buildings and an artificial vegetation layout shows a large impact of

2


even small vegetation obstacles like hedges [87]. This effect is also confirmed by studies in

built-up areas and open terrain simulations [73]. Recent studies with isolated mountain

islands in the sea show topographic speed-up effects when complex terrain is considered

in the simulations [57, 135]. These effects have to be taken into account for a precise wind

flow prediction.

Urban wind field simulations are not only performed with focus on wind energy but also

include a wide range of urban related topics. For example, Yoshie et al. [146] did some

validation studies with generic and real urban building geometries compared to wind tun-

nel experiments to investigate the wind flow at pedestrian level and proposed that the

results can be improved with scale-resolved simulations. The simulation of the Amster-

dam ArenA stadium by van Hooff et al. [56] showed large differences in the air change rate

of up to 42% depending on the wind direction. This rate can be improved by increasing

the size of openings near the roof of the stadium.

In new urban areas and new buildings, wind safety and wind comfort studies for pedes-

trians are required by some authorities. Thus, Blocken et al. [26] investigated the wind

flow around the Eindhoven University campus at pedestrian level and concluded that the

local wind speed and the wind direction can be totally different compared to the inflow

parameters and have to be simulated individually for each location. The wind in passages

between high buildings may be perceived as uncomfortable by pedestrians. Blocken et

al. [27] developed a control system, controlled by local wind measurements, with sliding

doors to modify the wind climate in passages.

Due to the climate change and rising temperatures, the effect of urban heat islands has

come more into the focus of the research in the recent years. Toparlar et al. [133] in-

vestigated the design of climate adaptive urban areas in a case study of a quarter in

Rotterdam. By means of urban greening and evaporative cooling the effect of a heat

storage mechanism and wind pattern on the temperature field is analysed. Besides veg-

etational evaporation, radiation plays a big role in an urban microclimate. Thus, Qu

et al. [97] showed the importance of a 3D radiation model to consider the non-uniform

radiation by the sun position and shadows which causes asymmetrical wall heating of

buildings. The microclimate in urban areas can also be influenced by the layout and

the morphologies of the buildings [18]. The temperatures of building facades depend on

the distance to other buildings with non-uniform heights. Near complex geometries the

facade temperatures are reduced due to increased shadowing effects. A higher ventilation

between the buildings cool down the facades.

Water surfaces can also lower the temperatures in urban areas as shown by Tominaga et

al. [131]. At pedestrian level the maximum temperature is decreased by the evaporation

of the water surface by approximately 2 ◦C. The wind propagates the cooling effect to an

unobstructed distance of 100m. The transpirational cooling effect of different types of

vegetation during a heat wave in the Netherlands was investigated by Gromke et al. [48].

The vegetation lowers the urban temperature by 1.6 ◦C and a facade greening by 0.3 ◦C

3


while a roof greening does not change the temperature in the street canyons.

Urban air pollution comes into the focus of research in the last recent years. Traffic is one

of the main contributors to air pollution. The polluted air is also distributed to streets

with lower traffic, not just limited to busy roads. The degree of pollution is locally in-

fluenced by buildings and vegetation, as shown by a concentration map for Madrid [108].

However, tree alleys which can be mostly found in every city, significantly changes the

flow in street canyons and can have an adverse effect on the pollution since they increase

the concentration near the walls of buildings [38]. Jeanjean et al. [61] investigated the

effectiveness of trees to disperse road traffic emissions in Leicester City Centre. At pedes-

trian level the pollution is decreased by 7% by tree alleys due to increased turbulence

and mixing effects. Trees show the biggest effect on reduction in emissions when they are

planted in open areas.

Another topic related to urban wind field simulations is the investigation of gas or other

chemical dispersions in populated areas. When gas is released due to terrorist attacks

or accidental events, the simulation could be used as a decision support system which

indicates the affected areas. Pontiggia et al. [95] simulated a real gas event and found

out that the presence of buildings leads to a different cloud footprint than in a free field

dispersion. Hanna et al. [53] investigated the release of chlorine gas from a railroad tank

car. The released chlorine cloud could initially extend a hundred meters or more in the

upwind and crosswind directions and follows the terrain drainage, e.g. river channels or

valleys. The buildings slow down the gas distribution and can locally increase the gas

concentration.

The three-dimensional building geometry used in this study is obtained from a CityGML

model. The availability of these city models, which are provided for whole Germany [49]

strongly has increased over the last years. The CityGML models are categorized into

different Levels of Details (LoD) from 0 to 4: LoD0 is actually just a 2D geometry con-

taining the ground area of buildings, LoD1 is the extruded geometry of LoD0 by the mean

building height. In LoD2, the simplified roof shapes, such as a flat roof or a gabled roof,

are included. The LoD3 contains more exterior details, e.g. balconies, dormers, windows

etc.. The highest Level of Detail LoD4 further includes interior rooms which can be used

for indoor ventilation studies. For urban CFD simulations usually LoD1 or LoD2 build-

ings are sufficient. With CityGML models different urban topics can be investigated. For

instance, Jang et al. [60] used these geometries to simulate the flooding status, the degree

of flooding and the level of building damage after heavy rainfall in a case study in Ko-

rea. Furthermore, the loss of human life and property damage is estimated in this study.

Trometer et al. [134] simulated scenarios in which detected unexploded bombs have to

be defused or detonated. This can be used for a more precise prediction of which areas

must be evacuated for the deactivation. One of the most important scenarios happening

4


in complex and tall buildings is the evacuation of people. Atila et al. [19] combined

CityGML buildings with an individual evacuation model to find the shortest path with

safety in case of disasters such as fire. Besides the heat demand forecast of single buildings

or city districts [88], CityGML models are applied in photovoltaic potential estimation.

Rodriguez et al. [100] determined their potential in urban and regional scale and also

the fraction of electricity demand that can be covered. Lu et al. [77] investigated a 3D

noise propagation on building facade level with different heights by using a 3D city model

and integrating all the noise coming from individual traffic such as cars or motorcycles as

well as planes and railroad based vehicles. Furthermore, studies based on 3D city models

or architectural proposal embedded in an existing city model can be used for decision

support in cities and communities.

Basically, two different approximate forms of the Navier-Stokes-equations can be used

for urban wind simulations: the Large Eddy Simulation (LES) in which the larger ed-

dies are solved directly and the smaller ones are modelled and the Reynolds Averaged

Navier-Stokes (RANS) simulation in which all eddies are modelled. In Detached Eddy

Simulations (DES), also called hybrid methods, the wall boundary layers in a domain are

modelled with RANS and the remaining flow field with LES. Blocken et al. [25] compared

both methods for the use of outdoor and indoor building simulations. In general, LES

shows more precise results, but more computational resources and higher requirements

for the grid generation are needed. Due to the high simulation complexity of LES a lack

of knowledge to set up such a simulation can yield to less accurate and less reliable results

than those by RANS [25]. Depending on the problem to be solved, RANS simulations

can be sufficiently accurate. Back to the basic bluff body flow simulation, Bourdreau et

al. [28] compared unsteady RANS (URANS) and DES simulations by performing simu-

lations of the wake of a bluff body. In that, the DES approach tends to overestimate the

time-averaged streamwise velocity component, especially in the near wake. The velocity

fluctations agree better with experimental data and better than the URANS results. For

a more detailed overview about bluff body flow it is referred to chapter 3.1. Krajnovic

et al. [68] also investigated the bluff body flow with Partially Averaged Navier Stokes

(PANS) simulations. They found out that the simulation results with LES and DES show

more details and fluctuations than PANS and RANS. Moreover, LES and DES tend to

be more accurate and more suitable for bluff body flow simulations. Liu et al. [75] also

confirmed the superiority of LES and DDES over RANS models when they simulated an

isolated high-rise building. LES and DDES predicted similar results in the wake region,

but the DDES approach has a lower overall mesh requirement. Thus, it is recommended

to use DDES for building simulations, since instantaneous wind characteristics can be

useful for a more accurate analysis of wind comfort [75].

5


The wind field in the atmospheric boundary layer is affected by the structure of the ground

surface which is complex terrain, vegetation and buildings. The effects of complex terrain

on the wind field under various wind directions were investigated by Huang et al. [57].

A hilly island in the sea is taken as a test site and the simulations are validated using

wind tunnel experiments. The island strongly influences the wind flow in terms of wind

speed and higher turbulent kinetic energy. Another study of topographic speed-up effects

in complex terrain using a hilly island in the sea was performed by Uchida et al. [135].

The wind field with different locations for high wind turbines is investigated. Behind the

island a recirculation area of eight times the height of the island is formed. The wind flow

is mainly affected by the terrain instead of the inflow shear boundary condition. Letzgus

et al. [73] performed highly resolved DDES simulations including complex terrain and

forested zones. A real inflow boundary condition with velocity fluctuations is used pro-

vided by a previous mesoscale simulation. The comparison with on-site met mast data

shows a good agreement regarding the power spectra and the mean velocity values.

For a precise prediction of the wind flow atmospheric wind properties have to be considered

in the simulations. The properties have to be applied as inflow conditions. A synthetic

turbulence generator is crucial for scale resolving simulations like DES. Li et al. [74]

compared two different inflow conditions: the first is a classical logarithmic/exponential

law to define the velocity profile without generating atmospheric turbulence at the in-

let. The second is to interpolate the velocities and turbulence properties provided by a

previous simulation of the upstream region. The last boundary condition shows a better

prediction, but also needs more computational resources. Schulz et al. [114] investi-

gated the power response of a wind turbine to inflow turbulence and terrain effects, and

confirmed the importance of applying atmospheric turbulence at the inflow plane. That

leads to an increase in load and power fluctuations and to a decrease of the tower blockage.

An approach for a more universal predication of urban wind field studies is to categorize

buildings and city quarters. Therefore, morphological indicators can be defined which po-

tentially correlate with the wind environment. The parameters such as building density,

building coverage in terms of the aspect ratio, variation level of the building volumes or

porosity of the city (street canyons, parks, etc.) have a large correlation with the wind

potential [140]. For example, the impact of the building length on the wind potential on

roofs is much more important than the building width [139].

The urban windfield is influenced by the arrangement of street canyons, the arrangements

and the shape of buildings. Wang et al. [141] defined various urban morphologies and cat-

egorized urban quarters. They defined seven typical urban forms from the city of Beijing

and made a cross analysis of the wind potential over the rooftops of the highest build-

6


ings in each urban form, wherein just the local dominant wind direction is simulated and

the wind energy density is evaluated. The determination of suitable locations for wind

turbines requires the use of real-scale urban geometries. Juan et al. [64] investigated the

wind resource assessment around high-rise buildings in real urban areas based on real ur-

ban geometries. The domain includes roof geometries, upstream obstacles, arrangements

of integrated building complexes and parallel high-rise buildings. Possible locations have

been evaluated depending on distances from rooftop sidewalls or lowest mounting heights

above rooftops. These zones indicate high wind power densities and acceptable turbulence

intensities for wind turbines.

A common approach evaluating wind turbine locations is calculating the mean wind speed,

for example in the above-mentioned studies fromWang et al. [141], Sunderland et al. [126]

and Balduzzi et al. [21]. Additionally, the turbulence intensity or the wind power density

can be considered, as for example in the studies of Ledo et al. [69] and Juan et al. [64].

Toja-Silva et al. [129] used the threshold value of the turbulent intensity of 0.15 according

to the European Wind Turbine Standards II [2] to define heights for installing HAWTs.

Another approach evaluating wind turbine locations is to calculate local Weibull distri-

butions of the wind field based on wind roses as annual wind speed distributions can be

approximated by a Weibull distribution [126]. To represent wind turbines in the evalua-

tion approach, the wind turbine power curve can be included in the approach [32, 43].

Kalmikov et al. [65] calculated the probability distribution of an urban wind flow which

is described by the parameters of the Weibull distribution. They used the MIT campus

to study and validate the mean wind speed and wind power density by integrating local

wind measurements and observations of near reference sites in their simulations. The

micro-climate around the campus was analysed and optimal locations for small wind tur-

bines were studied. Sunderland et al. [126] used Weibull distributions for a more accurate

power prediction of wind turbines than current wind turbine power output measurements

which are based on average wind speeds over an observation period. Two models were

developed which are based on the normal and Weibull distributions, respectively. Both

models predicted the mean wind speed, the standard deviation within a 10min time in-

terval and the turbulence intensity. The validation shows that the Weibull based model

shows a more accurate prediction of the energy production under realistic conditions.

However, Toja-Silva et al. claims that Weibull distributions of wind statistics would be

incorrect in urban settings [128].

The installation site has a great influence on the potential energy yield. Balduzzi et al.

[21] did a study with generic buildings which vary in shapes, geometrical proportions

and various arrangements. The skew angle of the flow, which may change near the roof

top, is considered in the wind turbine power curve. As a result, a vertical axis wind

turbine (VAWT) shows better performances with the skewed flow than a horizontal axis

wind turbine (HAWT), so that a VAWT may be more recommended for urban wind flow.

Gagliano et al. [43] assessed the feasibility of building integrated micro wind turbines and

7


calculated the wind distribution within urban areas and the yearly energy yield. That

allows a possible owner of a wind turbine to evaluate its effective potentiality of wind

energy generation in urban areas. A large number of building integrated micro wind tur-

bines, will make a large contribution to local energy production [43].

Dadioti did an extensive research about urban wind field simulation with the focus on

small wind turbines [33, 32]. In this study, the university campus in Leicester (UK) is

used as a test site for the DES simulations. After detailed validation studies with bluff

body flows and local anemometer data, the optimum location for micro wind turbines

installation was investigated based on the calculation of the annual energy production.

This study concludes and recommends that DES offers robustness and accuracy over a

range of wind conditions and is well suited to the analysis of wind energy potentials in

complex urban environments.

But one of the most interesting questions regarding small wind turbines in urban areas

is the question if it is worth to mount a small wind turbine on a specific building. To

answer this question, different micrositing processes were applied in the literature. For

example Wang et al. [142] performed CFD simulations of wind flow in built environ-

ment under the urban atmospheric boundary layer inflow condition. Their simulations

have been compared to wind lidar measurements. They concluded that the wind turbine

should be installed at the height of 1.3 to 1.5 of a building height. An overview of existing

studies gives the review paper of Toja-Silva et al. [128]. The studies considering several

real buildings and mentioned therein use mostly RANS simulations, a real-scale geometry

and isolated or several generic buildings. Toja-Silva et al. [128] emphasized the influence

and the importance of surrounding buildings and figured out that it is not legitimate

to conclude from the flow around one isolated building to the wind field around several

buildings. The few studies mentioned in the review paper did not consider a combination

of vegetation, the flow from more than one wind direction, and any kind of evaluation of

the suitable positions for small wind turbines based on AEP values on buildings or local

wind statistics as a whole in one study.

But only the minority of the studies used the power curve of a wind turbine to determine

the possible energy output such as Gagliano et al. [43] or Dadioti [32]. As Toja-Silva et

al. mentioned in their review paper there is a lack of studies which determine the local

frequencies of the wind speed and show the expected great potential of LES for accurate

turbine performance evaluation [128].

The energy yield of a wind turbine strongly depends on the wind flow around its position.

This wind flow is affected by the building itself on which the wind turbine is mounted. As

compared to typical rectangular roofs, a rounded roof design produces a lower turbulence

intensity and higher power density which is increased up to 86.5% [145]. Furthermore,

high-rise buildings tend to block the incoming wind in the upstream direction and induce

8


higher turbulence intensities [145]. Another approach to improve the wind flow towards

building integrated wind turbines are ducted openings. Their aerodynamic design is op-

timized for wind energy harvesting by Ruiz et al. [16]. Rounded openings suppress the

flow separation, enhance the magnitude and uniformity in the wind, and tend to lower

the turbulent kinetic energy. With an optimized design the wind speed can be increased

up to 78%. Ledo et al. [69] found out that the power density above flat roofs is increased

compared to pyramidal roofs, independently of the wind direction. In general, the wind

velocity decelerates above pitched and pyramidal roofs when the inclination angle of the

roof exceeds 55 ◦ in the pitched case and 67 ◦ in the pyramidal roof. With a domed, vaulted

and a wedged roof shape an increase of energy yield of more than 45% can be reached

[9]. The lowest increase of power is obtained with a pyramidal roof [9]. Moreover, Lu et

al. [78] mentioned that the concentration effect of buildings and the heights of buildings

could enhance wind power utilization by increasing the wind speed by a factor of 1.5 –

2. Thus, not only the type of roof affects the wind flow but also the interaction between

adjacent buildings.

It is well known that numerical simulations still rely on experimental data sets to validate

numerical setups, turbulence models, or any kind of new implementation in the source

code. For urban wind simulations many on-site measurement devices are needed to cap-

ture the wind flow in urban domains. Due to the large influence of local conditions the

measured wind data is only valid within a short range. Wind tunnel experiments can help

to overcome these problems. For simple validation purposes, wind tunnel data with sim-

ple building geometries are preferred. Yoshie et al. [146] from the Architectural Institute

of Japan proposed a workflow with different validation steps based on comparative stud-

ies. In the first step, the flow around two types of single high-rise buildings is simulated.

The second step includes a high-rise building surrounded by an array of equal simplified

buildings. Just the last step is based on actual urban areas with two types of building

complexes. But in all stages the complex terrain and the vegetation is not included [146].

To include complex terrain, vegetation and real scale buildings in the validation process,

on-site measurements are crucial but still lacking for validation purposes.

1.2 Objective

The objective of this work is to investigate the urban wind field at the university cam-

pus Morgenstelle in Tübingen (southern Germany) with the focus on the energy yield

of small wind turbines. Highly resolved DES simulations are preformed which include

building geometries, complex terrain and vegetation using the ANSYS Fluent software

[1]. The building geometries are based on a real 3D city model geometry. A method for

using 3D city models for urban CFD simulations from geometry optimisation to meshing

9


requirements is presented. An already existing vegetation model from Shaw and Schu-

mann [115] is extended to consider local tree heights and seasonal effects of the vegetation.

A sensitivity analysis is performed to study the influences of different meshing parameters

and domain sizes. Additionally, the effect of various environmental parameters such as

turbulence intensity at the inlet, different tree heights and foliage densities to simulate

a summer and winter wind flow are investigated. The wind flow is evaluated with and

without buildings which surround the target zone. The numerical setup and the vege-

tation model are validated with on-site LiDAR measurement data in Tübingen. A new

method for using planar on-site wind measurement data as an inflow boundary condition

is introduced.

To investigate if it is worth to install small wind turbines on one of the buildings, the

wind fields from four different wind directions are simulated and evaluated. A new method

for calculating local wind statistics for each numerical cell is introduced using large-scale

synthetic wind statistics. The local wind statistics are used to determine the maximal

annual energy production associated to each building using a real vertical and horizon-

tal axis wind turbine. The annual energy production values are corrected according to

local turbulence intensities and angles of attack. Finally, the annual energy production

is taken as the main parameter to evaluate suitable locations on the buildings for small

wind turbines.

10


2 Fundamental Equations and

Numerical Approach

The simulations in this work were performed with the commercial CFD program ANSYS

Fluent in the versions 18.1 [12] and 2019R3 [13]. For the highly resolved Detached Eddy

Simulation (DES) an incompressible solver is used with the Improved Delayed DES (ID-

DES) shielding function combined with the k-ω SST turbulence model. At the inflow

plane the atmospheric turbulence is synthetically generated. Their governing equations

and the methodologies are described in the following chapter.

2.1 Governing equations

A Newtonian fluid can be described by the Navier-Stokes equations which include the

conservation equations for mass, momentum and energy. The simulations in this work

are performed assuming incompressible flow only. Hence, the continuity equation [111] is

described by
∂ui
∂xi

= 0 (2.1)

where i = 1, 2 and 3 which corresponds to the x, y, z-direction in the computational

domain and ui is the velocity in i−direction. The momentum equation is defined by [123]

ρ
Dui
Dt

= ρgi −
∂p

∂xi
+ µ

∂2ui
∂xj∂xj

, (2.2)

where j and i = 1, 2 and 3, ρ is the density, p is the pressure, t is the time, gi is the

gravitational acceleration in i−direction and µ is the dynamic viscosity. According to

Pope [96] the material derivative is defined by

D

Dt
≡ ∂

∂t
+ uj

∂

∂xj
. (2.3)

Since the air temperature is not considered in this work, the energy conservation equation

is not listed here. The Navier-Stokes equations are solved numerically which is described

in the following.

11


2.2 Methodology of simulations

CFD simulations can be basically divided into three major categories according to the

degree of modelling: Reynolds-Averaged Navier–Stokes simulations(RANS), Large Eddy

Simulation (LES) and Direct Numerical Simulations (DNS). In RANS simulations the

variables are time-averaged and no information about fluctuating velocities is available.

The entire flow and all turbulent eddies are modelled. Due to low computational costs

RANS simulations are widely used in industrial applications. In LES small scale eddies are

filtered out by a filter operation in order to be modelled by a turbulence model. However,

the remaining large scale eddies are resolved directly. Since the large scale eddies are

the most energy containing eddies and are mainly affected by the geometry [96], e.g. in

urban areas which are mostly affected by buildings and vegetation zones [40, 96], LES

simulation results are assumed to be more accurate [25]. In DNS simulations 100% of the

turbulence is resolved directly, down to the smallest scales, also known as the Kolmogorov

scales [96]. That makes DNS simulations the most expensive type of simulation and is

especially used for validation purposes. In general, the degree of modelling increases the

computational costs, but also the accuracy of the results [96].

In a hybrid simulation, both RANS and LES are involved simultaneously, also called

Detached Eddy Simulation (DES). The basic idea behind DES is that the near-wall region

is simulated with the RANS method and the outer, i.e. detached flow, with the LES

method. The transition from RANS to LES is coupled by the length scale and the

viscosity. When the boundary layer is modelled, the cell layers inside the boundary

layer can be thicker compared to when the boundary layer is solved directly because the

small turbulent eddies within a boundary layer do not have to be resolved in LES mode.

That increases temporal resolution in the detached regions and reduces computational

costs [122]. Hence, DES can be seen as a preferable compromise for urban wind field

simulations, since the focus is on the directly resolved outer region where wind turbines

are mounted on buildings and unsteady turbulence information can be provided with

reasonable costs.

Before the LES method and the transition between RANS and LES is explained, the

RANS approximation of Navier-Stokes equations is described briefly. A velocity u is

decomposed as the sum of a time-averaged velocity ū and a fluctuating velocity u′. After

inserting u = ū+ u′ in equ. 2.2 and time-averaging this equation yields to

ρ
∂ūi
∂t

+ ρūj
∂ūi
∂xj

= − ∂p̄

∂xi
+

∂

∂xj

(
µ
∂ūi
∂xj

− ρu′iu
′
j

)
, (2.4)

where the term −ρu′iu′j is the unknown Reynolds stress tensor. Defining a conservation

equation for−ρu′iu′j would lead to further unknown terms in this equation. More equations

would lead again to more unknown terms resulting in an endless loop which cannot be

closed [96].

12


The starting point for solving the closure problem is Boussinesq’s analogy to the law of

viscosity. The turbulent shear stress is directly proportional to the mean deformation

rate. In tensor notation, this analogy results

τij = −ρuiuj = µt

(
∂ūi
∂xj

+
∂ūj
∂xi

)
− 2

3
ρδijk, (2.5)

where τij is the shear stress tensor and µt the modelled turbulent viscosity. Thus, the

modelling of the turbulent stresses is reduced to the modelling of µt. In this work, the

k-ω SST approach by Menter et al. [83] is used as the RANS model for DES simulation

[10, 11]. Accordingly, the turbulent viscosity is modelled by the following term

µt = Cµ
k

ω
(2.6)

Since µt is a function of k and ω, the Reynolds stress tensor is modelled by defining two

additional transport equations, one for the turbulent kinetic energy k and one for the

specific dissipation rate ω. The original k-ω SST approach, developed by Menter [84],

uses the k-ω model by Wilcox [143] in the sub- and log-layer and gradually switches to

the k-ε model in the outer part of the boundary layer [84], wherein ϵ is the dissipation

rate. The governing equation for k implemented in ANSYS according to the formulation

of Menter et al. [83] is

∂(ρk)

∂t
+
∂(ρuik)

∂xi
= P̃k − β∗ρkω +

∂

∂xj

[(
µ+

µt

σk

)
∂k

∂xj

]
(2.7)

with the limited production term P̃k

P̃k = min

(
τij
∂ui
∂xj

, 10 · β∗ρkω

)
, (2.8)

where µ is the dynamic viscosity and β∗ and σk are constants. The transport equation

for ω implemented in ANSYS is

∂(ρω)

∂t
+
∂(ρuiω)

∂xi
= αωρS

2−βρω2+
∂

∂xj

[(
µ+

µt

σω

)
∂ω

∂xj

]
+2 (1− F1)·

ρσω2
ω

∂k

∂xi

∂ω

∂xi
(2.9)

with the strain rate magnitude S and the constants αω, β, σω and σω2 [83]. The blending

function F1 ensures that the k-ω model gradually switches to the k-εmodel with increasing

wall distance.

The outer region in DES is simulated with LES. Herein, the eddies are filtered in the

computational domain so that large scale eddies can be resolved directly, while smaller

scales can be represented by simple models [96]. The general filtering operation G(x, x′)

13


applied to a quantity ϕ(x′) is introduced by Leonard [72] and is defined by

ϕ̄(x) =

∫
G (x, x′)ϕ (x′) dx′ (2.10)

to obtain the filtered quantity ϕ̄(x) [96]. The filter width determines the filtered length

scale l which is not necessarily the cell size ∆ but the condition l > ∆ is always valid

[40]. In ANSYS Fluent, the filtered length scale is based on the local grid scale lg and

is calculated by lg = V
1/3
c , where Vc is the volume of the computational cell [10, 11]. By

introducing the subgrid stress (sgs) tensor τ sgsij [96] based on the filtered velocities in equ.

2.10

τ sgsij = −ρu′iu′j, (2.11)

the filtered momentum equation for LES yields [40]:

∂ (ρūi)

∂t
+
∂ (ρūiūj)

∂xj
=

∂

∂xj

[
µ

(
∂ūi
∂xj

+
∂ūj
∂xi

)]
−
∂τ

sgs
ij

∂xj
− ∂p̄

∂xi
. (2.12)

Similar to the Reynolds stress tensor, the sgs tensor has to be modelled to close the

equations, similar to the RANS formulation. In ANSYS Fluent, the Boussinesq hypothesis

is applied to the subgrid models to compute the subgrid-scale stresses with the strain-rate

tensor Sij [10, 11]

τ sgsij − 1

3
τ sgsbb δij = µt

(
∂ūi
∂xj

+
∂ūj
∂xi

)
= 2µtSij, (2.13)

where b, i and j = 1, 2 and 3. Finally, the turbulent viscosity µt is the only left quantity

which has to be modelled.

In the newer version of DES by Spalart et al [121], also called the Delayed DES (DDES),

a shielding function is introduced to protect the outer LES region to fall back into the

RANS modelled boundary layer. When an undesired LES-fallback into the boundary

layer occurs, velocity fluctuations cannot be resolved properly because the grid of the

boundary layer is too coarse. But providing velocity fluctuations for the LES mode is

crucial. It is achieved by Menter et al [83] by using the F1 and F2 functions of the k-ω

SST RANS model to identify the boundary layer and the transition to LES since the

thickness of the boundary layer is unknown during the grid generation process.

Since the simulations in this work are performed with the IDDES k-ω SST model, this

model is described briefly. An overview of all tested DES turbulence models and their

different results is given in the next chapter. The IDDES model in ANSYS Fluent is based

on the k-ω SST model [84] and the IDDES approach by Shur et al [118] including the

modifications proposed by Gritskevich et al. [47]. This model provides shielding against

the LES-fallback and includes wall-modelled LES capabilities in the attached flow regions

[47]. The following governing equations are extracted from the publication of Gritskevich

14


et al. [47]. The equation for k is defined as

∂(ρk)

∂t
+
∂(ρuik)

∂xi
=

∂

∂xj

[(
µ+

µt

σk

)
∂k

∂xj

]
+ Pk − ρ

√
k3

lIDDES

, (2.14)

where lIDDES is the IDDES length scale, and with the limited production term Pk

Pk = min
(
µ2
tS

2, 10 · Cµρkω
)

(2.15)

and the constant Cµ. The turbulent viscosity is modelled as

µt = ρ
a1 · k

max (a1 · ω, F2 · S)
(2.16)

with the constant a1 and the blending function F2. While the ω-equation 2.9 and the

blending functions F1 and F2 remain unmodified, the sink term of the k-equation is

modified using the IDDES length scale lIDDES which is defined as follows

lIDDES =f̃d · (1 + fe) · lRANS +
(
1− f̃d

)
· lLES

lLES =CDES ·min {Cw ·max [dw, hmax] , hmax}

lRANS =

√
k

Cµω

CDES =CDES1 · F1 + CDES2 · (1− F1)

(2.17)

with the constants CDES, CDES1, CDES2, Cw, the wall distance dw, the LES length scale

lLES and the RANS length scale lRANS. The variable hmax is the maximum edge length of

a cell and its definition includes only hexahedral cells [47]. This formulation is adapted

to other cell types in ANSYS Fluent [10, 11]. For the empiric blending function f̃d and

the elevating function fe it is referred to the work of Gritskevich et al. [47].

2.3 Generation and injection of atmospheric

turbulence

In order to simulate the atmospheric turbulence in the wind field, the turbulent structures

have to be generated synthetically at the inflow plane of the domain in scale-resolving

simulations. ANSYS Fluent provides two different methods to generate turbulent velocity

fluctuations: the Vortex Method (VM) and the Spectral Synthesizer (SpS) [10, 11]. In

most of the simulations of this work the Vortex Method is used, thus only the algorithm of

this method is described in the following chapter which is mainly related to the ANSYS

Fluent Theory Guide [10, 11]. This approach is originally published and validated in

several cases by Mathey et al. [81, 80] and is a relatively inexpensive and precise way to

15


generate random fluctuations at the inlet since it is temporally and spatially correlated

[81].

To consider perturbations of a turbulent inlet velocity profile, a time-dependent 2D vortex

method is applied. The fluctuating velocities are added to the mean velocity profile at

the inflow plane. A vorticity field is created randomly on the two-dimensional boundary

inlet plane. This approach is based on the Lagrangian form of the 2D evolution equation

of the vorticity written in the Biot-Savart law form [11, 123]. This equation is solved by a

particle discretization representing the vortex points. Taking the number of vortex points

N and the area of the inflow plane Ain into account, the amount of vorticity carried by a

given particle i is represented by the circulation Γi(x, y)

Γi(x, y) = 4

√
πAink(x, y)

3N [2ln(3)− 3ln(2)]
(2.18)

and a spatial distribution η depending on the location vector x⃗

η (x⃗) =
1

2πξ2

(
2e−|x|2/2ξ2 − 1

)
2e−|x|2/2ξ2 . (2.19)

The quantity ξ specifies the size of a vortex by a turbulent mixing length hypothesis to

ensure a wide applicability. Using k and the dissipation rate ε it yields

ξ =
cVMk

3/2

2ε
, (2.20)

where the constant cVM amounts to 0.16. The discretization of the velocity field is given

by

u⃗ (x⃗) =
1

2π

N∑

i=1

Γi

((x⃗i − x⃗)× z⃗)
(
1− e|x⃗−x⃗′|2/2ξ2

)

|x⃗− x⃗′i|2
(2.21)

with the unit vector z⃗ in streamwise direction. The sign of the circulation of each vortex

changes randomly every characteristic time scale which is the time for a vortex to travel

100 times the mean vortex size ξ with the velocity in streamwise direction [81]. The

minimum size of the vortex is bounded by the cell width which ensures that it belongs to

resolved turbulent scales which is important for partially resolved simulations as LES or

DES.

A simplified Linear Kinematic Model (LKM) is used based on the work of Mathey et al.

[10, 11, 81]. This model imitates the influences of the vortices in the streamwise mean

velocity field. The fluctuating velocity field u′ in streamwise direction is given by

u′ = −v⃗′ · n⃗, (2.22)

16


where v⃗′ is the planar fluctuating velocity field and n⃗ is the normal vector of the gradient

of the streamwise mean velocity. The perturbations are equally distributed among the

velocity components. In the case normal fluctuations are known, e.g. from previous

simulations or defining the Reynolds stresses for isotropic turbulence, they can be applied

to the flow field in order to fulfil the normal statistic fluctuations uiui [10, 11]. The

resulting velocity fluctuations are then calculated by

u′∗i = u′i

√
uiui√
2/3k

. (2.23)

17


3 Testcase Validation: Flow around

a High-rise Building Structure

In this chapter different setups and turbulence models with various shielding functions

were tested and validated. Since no real-scale measurements for buildings are available,

a surface mounted cube in a boundary flow profile was chosen which is similar to a single

high-rise building. The results are validated with wind tunnel experiments by Bourgeois

et al. [29] and DNS simulations by Saeedi et al. [104].

3.1 Description of the experiment

Bluff body flows are very common for testing simulation setups and turbulence models.

Both can be validated by wind tunnel experiments and DNS simulations. For instance,

Yakhot et al.[144] performed a DNS simulation of a wall mounted cube in a fully developed

channel flow and analysed the wake and the vortex structure. Also Saeedi et al. [104]

used DNS simulations to investigate the wake of a surface mounted slim cube and studied

the Reynolds stresses and velocity fluctuations in the wake. Valger et al.[137] simulated

the flow around a surface mounted cube and discovered an over-prediction of k in RANS

with k-ϵ and k-ω turbulence models. Elkhoury et al. [37] did a bluff body flow study

with different RANS models. They came to the conclusion that the Spalart-Allmaras

turbulence model and the k-ω Scale Adaptive Simulation best reproduce the separation

region in front of the cube and the Spalart-Allmaras model predicts the best velocity pro-

file of the roof of the cube. But the length of the separation region behind the bluff body

is overpredicted by all models used in the study. Due to the complex flow around bluff

bodies a scale resolving simulation is highly recommended and shows a higher accuracy in

the results than using RANS [28, 75, 68]. Roy et al.[102] showed in their simulations that

RANS over-predicted the length of the recirculation zone in comparison to LES. Paik et

al.[92] simulated the flow around two wall-mounted cubes in tandem. Herein, URANS

fails to capture key features such as single vortices while good agreements with various

DES and DDES models are achieved. Robertson et al.[98] observed similar results testing

the DDES and the IDDES shielding function.

The flow around an array of several cubes was simulated using LES by Stoesser et al.[125]

18


who achieved a good agreement with wind tunnel experiments. Hanna et al.[54] simulated

various arrangements of several obstacles using LES to investigate street canyon effects.

Both Hang et al.[52] and Santiago et al.[107] investigated the flow around cube arrays,

validated their results with wind tunnel data and obtained a good agreement. The impor-

tance of a synthetic turbulence generator to capture the atmospheric turbulence in LES

simulations is confirmed by Shi et al.[116].

Also on the experimental side there are some studies which are worth to mention. Ha-

jimirzaie et al.[51] used particle image velocimetry in the wind tunnel to investigate the

wake behind two different ellipsoid bodies and observed counter-rotating distributions of

vorticity inducing downwash (tip structures), upwash (base structures), and horseshoe

vortices in the wake. Bourgeois et al.[29] investigated the large scale structures behind a

wall mounted square cylinder. They observed a strong interaction between the free end

and the wall and a resulting deformation of the vortex structures with a principal core

behind the cube. Two zones in the flow around a surface mounted cube are detected by

Sattari et al.[109]. While the first and windward zone is dominated by alternate forma-

tion and vortex shedding, the second and lee side zone is characterized by two co-existing

vortices throughout the shedding cycle.

The wind tunnel experiments by Bourgeois et al. [29] and Sattari et al. [109], who inves-

tigated a small-scale high-rise building in a wind tunnel, and DNS simulations by Saeedi

et al. [104] are taken to validate the numerical setup which is used for the simulations

in this work. The building with an aspect ratio of 4:1 in height to width is placed in a

boundary layer flow. A schematic view of this geometry is shown in fig. 3.1. The side

length d and the height of the cube Hc is 0.0127m and 0.0508m. The computational

domain is 0.31496m long, 0.17018m wide and 0.1143m high. Saeedi et al. [104] also

mentioned differences in the blockage ratio and in the domain size between the wind tun-

nel experiments and their DNS simulation. Since the results in this work are compared

in particular with the results of Saeedi et al. [104], the size of the domain is according

to the DNS domain of Saeedi et al. [104]. The cube is located 4 d downstream behind

the inflow plane and in the centre of the z-direction as shown in fig. 3.1. The boundary

layer flow enters the domain on the left with a free stream velocity u∞ of 15m/s. The

Reynolds number is 12.000, based on u∞ and d as the characteristic length scale. The

centre of origin is located at the ground area of the cube in its centre as illustrated in fig.

3.1.

3.2 Numerical setup

In order to exactly match the thickness of the boundary layer of 0.18 4 when it hits the

obstacle, the same formulation for the mean inlet velocity u(y) is used as Saeedi et al.

19


9 d

u∞

Free surface

13.4 d
4 d 19.8 dd

4 d x

y

z

Figure 3.1: Computational domain of the validation case [104] (modified).

[104]

u(y) =




u∞ ·

(
y

0.15Hc

)0.16

for y < 0.15Hc

u∞ for y ≥ 0.15Hc

(3.1)

A synthetic turbulence generator is applied to the inlet plane using a turbulent inten-

sity of 0.8% and a zero mean of the velocity fluctuations to reproduce equal turbulence

structures as in the DNS simulations. The velocity fluctuations are added to the mean

velocity profile of u(y) at the inlet plane, as mentioned in chapter 2.3. In ANSYS Fluent

two methods are implemented modelling fluctuating velocities at the inlet, the Vortex

Method and the Spectral Synthesizer [10, 11]. Both of them are tested to study their

influences on the results. The free surfaces are modelled as zero gradient boundary condi-

tions, and a no-slip boundary condition is applied to all solid surfaces. The outlet plane

is modelled as a pressure-outlet boundary condition. In the vortex method (VM) the

number of vortices has to be defined a priori, but is not published by Saeedi et al. [104]

and by Bourgeois et al. [29]. Thus, the value of 200 as the number of vortices is taken

from a similar bluff body flow simulation performed by Mathey et al. [81]. In order to

investigate the sensitivity of the vortex number, also 150 and 250 vortices considered.

In contrast, the Spectral synthesizer method does neither require any additional input

quantities nor a specific vortex number.

Pope [96] proposed that 80% of the turbulent scales in a flow have to be resolved in

scale-resolving simulations. The cell size has to be determined carefully since it limits

the minimal size of turbulent scales [40]. It determines the degree of modelled turbulent

structures and influences the simulation accuracy [96].

The following approach to calculate the cell width in an LES grid is proposed by Pope

[96]. According to Pope, 80% of the energy is contained in motions of length scales

1/6L11 < l < 6L11 where L11 is the longitudinal integral length scale. Regarding a

20


cumulative turbulent kinetic energy spectra, approximately 80% of the energy is resolved

in the range until lEI = 1/6L11 where lEI is the demarcation length scale between the

energy containing range of eddies (l > lEI) and smaller eddies. Considering the Reynolds

number of 12.000 in the validation case, the Taylor-scale Reynolds number Reλ is 155,

estimated from the relation

Reλ ≈
√
2Re. (3.2)

The ratio of L11/L, wherein L is the turbulent length scale L = k3/2/ε, depends on the

Reynolds number of the model spectrum Reλ, which is 155 in this case, and L11/L is 0.51

[96]. For 80% of the resolved turbulent kinetic energy, the product κL11 is 15 [96] where

κ is the wave number. Considering the definition of the cut-off wave number

κc =
π

∆
, (3.3)

the cell width is

∆0.80 =
π

κc
=
π · L11

15
=
π · 0.51L

15
= 0.107L. (3.4)

With L = 0.009m obtained in a previous RANS simulation the cell width is 0.001m for an

80% resolution. In the simulation of the surface mounted cube, only tetragonal cells are

used. This strongly simplifies the meshing process of urban domains including buildings,

vegetation and terrain. The grid has a total number of 6.235.529 tetragonal cells. In

comparison to that, the grid, which Saeedi et al. [104] used in their simulations, has 35.5

million cells. The boundary layer consists of 25 inflation layers with the first cell layer

height of 0.0000175m, fulfilling the requirement y+ = 1 for an entirely resolved boundary

layer at the wall, where y+ is the normalized wall distance. The growth rate of the cells

in the boundary layer towards the outer flow is 1.2. For a better comparison, the total

runtime of the simulations is eight flow-through times which consists of two and a half

times for stabilizing the flow and five and a half times for sampling unsteady statistics

according to the times used by Saeedi et al. [104]. The time step is chosen according to

the stability criteria of the Courant number Co < 1 [40, 111].

3.3 Results and validation

For the validation of the numerical setup, the flow around a surface mounted cube is sim-

ulated with different turbulence generators and DES turbulence models combined with

various shielding functions. Regarding synthetic turbulence generation, the Spectral Syn-

thesizer and the Vortex Method with 150, 200 and 250 vortices are applied. In the case

no turbulence generation method is explicitly mentioned, the Spectral Synthesizer is used

in the simulations. On the turbulence model side, the k-ω SST [84], the realizable k-ϵ

[117] and the one-equation Spalart-Allmaras [120] turbulence model are used for the in-

21


vestigation. To study the effect of the shielding functions on the flow field, the following

shielding functions are used in the validation study: DDES [121], IDDES [47, 118] and

the Shielded Detached Eddy Simulation (SDES) [10, 11].

3.3.1 Turbulence model and shielding function

In fig. 3.2a the dimensionless time-averaged velocity in streamwise direction, located at

x/d = 3.5 and y/d = 3 is plotted against the dimensionless z-axis. All evaluation lines

are located in the wake behind the cube. The best agreement with the DNS simulations

is achieved with the SDES and the IDDES model, while the velocities obtained with the

first model are closer to the measured velocity at z = 0. The other turbulence models

with the DDES shielding function, show an insufficient performance and underestimate

the measured velocity by 0.3 at z=0. The models with the DDES shielding function

predicted a wider and longer wake behind the cube which leads to a higher blockage and

consequently to higher velocity peaks at the edge of the wake. Roy et al. [102] made a

similar observation in their study. Herein, good agreements between LES simulations and

experiments were achieved, while with RANS the recirculation zone was predicted too

long. A thicker RANS boundary layer around the cube in the DDES shielding function

could be an explanation for the longer recirculation zone. The deviation between the

numerics and the experiments in the outer regions of the domain can be traced back to

the slightly different domains and the higher blockage ratio in the simulations [104], as

mentioned in the description of the experiment.

The u′xu
′
z element of the Reynolds stress tensor, evaluated along a line at x/d = 3.5 and

y/d = 3, is shown in fig. 3.2b. The best agreement is obtained with the IDDES model.

In the left peak in fig. 3.2b, the u′xu
′
z-stresses are almost identical to the DNS data. The

DDES k-ω SST and the SDES models show an overshoot in the results, while the DDES

k-ϵ realizable model basically struggles to resolve the turbulent fluctuations in the wake.

The right peak in fig. 3.2b is over-predicted by almost all turbulence models except the

DDES k-ϵ realizable model. With the Spalart-Allmaras and the k-ϵ realizable models the

peak is narrower and is located more outside in spanwise direction.

In fig. 3.2c the time-averaged velocity ux is plotted in downstream direction at y/d = 1

and z = 0. The best results compared to the DNS simulation are obtained with the

DDES Spalart-Almaras and the IDDES model. Both turbulence models fail to predict

the location of the negative peak velocity correctly. The recirculation zone is predicted

to be shorter with the IDDES model, but longer with the DDES Spalart-Almaras model.

However, with increasing distance to the cube the simulated velocity with the IDDES

model is closer to the velocity in the DNS simulations than the velocity with the DDES

Spalart-Almaras model.

22


-6 -4 -2 0 2 4 6

z/d [-]

-0.6

-0.4

-0.2

0

0.2

0.4

0.6

0.8

1

1.2

1.4

u
x
/u

[-
]

Exp

DDES k- SST

DDES k- realizable

DDES Spalart-Almaras

SDES k- SST

IDDES k- SST

DNS

(a)

-6 -4 -2 0 2 4 6

z/d [-]

-0.1

-0.08

-0.06

-0.04

-0.02

0

0.02

0.04

0.06

0.08

0.1

u
' x

u
' z

/u
2

[-
]

Exp

DDES k- SST

DDES k- realizable

DDES Spalart-Almaras

SDES k- SST

IDDES k- SST

DNS

(b)

0 2 4 6 8 10 12 14

x/d [-]

-0.5

0

0.5

1

u
x
/u

[-
]

Exp

DDES k- SST

DDES k- realizable

DDES Spalart-Almaras

SDES k- SST

IDDES k- SST

DNS

(c)

0 2 4 6 8 10 12 14

x/d [-]

-0.05

0

0.05

0.1

0.15

0.2

0.25

0.3

0.35

u
x

,r
m

s
/u

[-
]

Exp

DDES k- SST

DDES k- realizable

DDES Spalart-Almaras

SDES k- SST

IDDES k- SST

DNS

(d)

Figure 3.2: Time-averaged first- and second-order velocity profiles in x-direction and el-
ements of the Reynolds stress tensor: (a) ux/u∞ at x/d = 3.5 and y/d = 3,
(b) u′xu

′
z/u

2
∞ at x/d = 3.5 and y/d = 3, (c) ux/u∞ at y/d = 1 and z = 0, (d)

ux,rms/u∞ at y/d = 1 and z = 0 in comparison to the wind tunnel experiments
by Bourgeois et al. [29] and the DNS simulations by Saeedi et al. [104].

Fig. 3.2d shows the dimensionless root mean square velocity in streamwise direction

ux,rms. The best qualitative and quantitative agreement with the DNS data is obtained

with the IDDES model. The ux,rms values almost reach the maximal peak value and further

decrease slightly towards the outlet with the same slope as in the DNS simulations. The

DDES k-ϵ realizable model completely fails to predict the initial slope of ux,rms behind

the cube currently due to the over-prediction of the length of the recirculation zone which

is also confirmed by the streamwise velocity distribution in fig. 3.2c. In the simulation

with the DDES k-ω SST model, the maximum value of ux,rms is over-estimated by 50%

and by 25% with the SDES and the DDES Spalart-Allmaras model. The peak of ux,rms

in fig. 3.2d indicates the reunion of the flow which is separated by the cube after the

recirculation zone. The location of the peak is moved downstream using the turbulence

23


models with the DDES shielding functions. The value of ux,rms at the outlet is quite

similar independently of the used turbulence models.

To analyse the entire vortex structure in the domain, the vortices are visualized using iso-

surfaces of the λ2-criterion for the DDES (in fig. 3.3) and the IDDES shielding function

(in fig. 3.4) with the same k-ω SST turbulence model. The λ2-criterion was developed

by Jeong and Hussain [63] and is a method to visualize vortex structures from a three-

dimensional velocity field. The variable λ2 defines the second eigenvalue of the tensor

S2 + Ω2, where S and Ω are the symmetric and antisymmetric tensor of the velocity

gradient, respectively.

In the figs. 3.3 and 3.4, the horseshoe vortex is a characteristic vortex which is formed

u/u∞
0 0.2 0.4 0.6 0.9 1.1 1.3 1.5

u

Figure 3.3: Iso-surface of λ2 for the DDES
shielding function.

u/u∞
0 0.2 0.4 0.6 0.9 1.1 1.3 1.5

u

Figure 3.4: Iso-surface of λ2 for the IDDES
shielding function.

near a surface mounted cube at the ground wall [58]. It is generated by the backflow

of the flow which hits the wall pointing against the streamwise direction and the flow

which passes the cube at its sidewalls. The flow in front of the cylinder is subjected to

a pressure increase. As a result, the boundary layer in front of the cylinder detaches

immediately from the bottom and rolls up into a vortex, which wraps around the cylinder

like a horseshoe [58]. In the IDDES simulation, the horseshoe vortex is developed more

precisely.

The IDDES model is more capable to resolve finer turbulent structures which is especially

visible in the downstream part of the flow and in the boundary layer near the ground wall.

The vortices in the separation zones near the cube show more detailed structures and finer

resolved eddies in the IDDES simulation, since the IDDES model includes wall-modelled

LES capabilities in the near wall regions and provides better shielding against the LES

fallback into the RANS modelled boundary layer. For more detailed information about

the IDDES model it is referred to the literature [47, 118].

24


3.3.2 Turbulence generation method

In this subchapter the bluff body