# A source-term formulation for injecting wind gusts in CFD simulations

Guillaume de Nayer, Michael Breuer

► To cite this version:

Guillaume de Nayer, Michael Breuer. A source-term formulation for injecting wind gusts in CFD simulations. Journal of Wind Engineering and Industrial Aerodynamics, 2020, 207, pp.104405. 10.1016/j.jweia.2020.104405. hal-02990346

HAL Id: hal-02990346

https://hal.science/hal-02990346

Submitted on 5 Nov 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# A source-term formulation for injecting wind gusts in CFD simulations

G. De Nayer, M. Breuer*

Professur für Strömungsmechanik, Helmut-Schmidt-Universität Hamburg, D-22043, Hamburg, Germany

<figure><img src="image_3.png" /></figure>

## ARTICLE INFO

**Keywords:**

Wind gust

Source term

Injection

Large-eddy simulation

ECG

EOG

## ABSTRACT

The objective of the present paper is to develop a methodology to inject strong wind gusts into the computational domain in order to efficiently simulate their effect on the fluid flow. The design of the methodology based on a source-term formulation takes the feedback effect of the resulting turbulent flow (and, if present, the impacted structure) on the wind gust itself into account. Since the injection of the wind gusts can be carried out close to the region of main interest, CPU-time intensive methods to ensure a proper transport of the gust through the flow field can be avoided. The methodology is mainly intended for the application within eddy-resolving simulations (e.g., LES), but it is not restricted to this class of simulation approaches. For the description of the gusts classical shape functions such as the Extreme Coherent Gust (ECG) and the Extreme Operating Gust (EOG) as well as a newly derived $C^2$-"1-cosine" shape are applied. Two scenarios are taken into account to assess the proposed gust injection technique. On the one hand a (laminar) undisturbed flow field is considered and the effect of different time and length scales of the gusts on their evolution and propagation through the flow field is studied in detail. On the other hand a turbulent background flow is assumed demonstrating that the methodology suggested is also applicable for practically relevant turbulent flows.

## 1. Introduction

Singular wind events such as tornados, downbursts or wind gusts are not only of interest for meteorologists who are trying to analyze, understand and foresee such extreme phenomena (e.g., Feser et al., 2015). They have also garnered the attention of engineers due to their devastating consequences on structures. An impressive example of extreme structural damages on wind turbines observed at a wind farm in 2011 was investigated in details by Hawbecker et al. (2017). A wind turbine tower was buckled and blades of several others turbines were stripped away during a storm. Simulations including transient wind gusts were performed and compared to available wind datasets. Especially constructions made of lightweight thin structures (e.g., large umbrellas, tents, stadium roofs) are vulnerable to wind events. The National Institute of Standards and Technology (2020) (NIST, USA) gathers disaster studies related to windstorms in the USA since 1969. A variety of similar reports are available. For example, Yokel et al. (1976) present the wind damages on buildings in the metropolitan area of Washington DC during a storm. The report compiles the different destructions but also the meteorological data. More recently in 2009, the Dallas Cowboys indoor practice facility, a membrane-covered frame structure, collapsed during a severe thunderstorm. The facility including the structural design and

analysis, its destruction and the wind environment during the storm are described in detail in Gross et al. (2013). Massive wind destructions also occur in nature, particularly forests leading to typical damage patterns. Consequently, wind tunnel measurements on the flow and the gust dynamics as well as studies on the modeling of these phenomena can be found in the literature (see, e.g., Gromke and Ruck, 2018; Kamimura et al., 2019; Tischmacher and Ruck, 2013).

In order to prevent structural damages due to the wind, design standards are prescribed. For example, the engineering handbook of Frost et al. (1978), the IEC-Standard 61400-21 IEC-Standard (2002) and the wind energy handbook of Burton et al. (2001) provide guidelines for the construction of wind turbines. These are often based on measurements and estimations such as in the work of Kasperski (2007). However, the assumptions made in the guidelines are rather simplistic. Therefore, thorough investigations of the proposed design submitted to different wind conditions either in experimental facilities and/or by the application of sophisticated highly-resolved numerical simulations is a must for a technological advancement in the future.

Due to the extremely wide field the following considerations are restricted on the one hand to numerical simulations and on the other hand to wind gusts occurring in the micro scale, e.g., from meters to one hundred meters (Steyn et al., 1981). Methods to describe these wind

<footnote>* Corresponding author.</footnote>

<footnote>E-mail address: breuer@hsu-hh.de (M. Breuer).</footnote>
gusts can be classified into *deterministic* and *stochastic* models.

In the *deterministic* model-based method the mean shape of the gust is defined by a mathematical function. Two well-known shapes were proposed in the IEC-Standard (2002) and prescribed by the Federal Aviation Regulation: The Extreme Operating Gust (EOG) also called Mexican Hat and the Extreme Coherent Gust (ECG) also denoted 1-Cosine shape. A comparison of selected real gust events with the idealized EOG is available in Emeis (2018). A mean shape in form of a Gaussian distribution was also recently applied by Yawar et al. (2019). Another more realistic (one- and two-dimensional) model was proposed by Knigge and Raasch (2016) derived from turbulent flow fields computed by large-eddy simulations (LES). In general, the *deterministic* representation of a gust by a prescribed mathematical function is easy to realize computationally, but also simplistic. Thus, it often does not fully represent the reality (Bos, 2017; Branlard, 2009).

The *stochastic* model-based method tries to represent the gust in a more physical manner, see, e.g., Bierbooms (2005, 2006, 2007); Bierbooms and Cheng (2002); Bierbooms et al. (1999); Bos (2017); Larsen et al. (2003) and Nielsen et al. (2003). Here the gust is not generated by superposition of a deterministic mean shape to the velocity signal, but the velocity is constrained to a prescribed gust velocity at a pre-defined time, so that the whole signal fulfills the correlations of the isotropic turbulence theory. The resulting constrained stochastic gust model was compared with classical EOG and successfully applied to simulations for the flow around wind turbines.

Since the main objective of the present contribution is to describe a method how wind gusts can be imposed on the fluid field, this aspect is described in more detail in the following. The classical method denoted *far-field boundary condition* is to inject wind gusts as inlet boundary conditions. Since the gust superimposed at the inlet boundary induces a strong velocity perturbation, for an incompressible fluid it is necessary to dynamically adapt the velocity at the inlet boundary far away from the gust, so that the global inlet mass flow remains constant. Norris et al. (2010) proposed a spatial and temporal correction of the inlet velocity profile to fulfill this requirement in their LES carried out to investigate the three-dimensional flow and turbulence structures in the wake of a wind turbine impacted by an EOG. Later the same group studied the turbulent flow around three wind turbines introducing ECG at the inlet boundary (Storey et al., 2013, 2014). Besides LES a variety of simulations relying on URANS exist requiring lower resolutions in space and time. For example, the flow around a vertical axis wind turbine undergoing an EOG was computed in 2D by Onol and Yesilyurt (2017). Various EOG were introduced at the inlet and different amplitudes were considered. Similarly, relying on 3D-URANS simulations, Menegozzo et al. (2018) studied the response of a horizontal axis wind turbine undergoing an EOG. As before the gust was introduced at the inlet boundary and freely propagated through the domain. The CFD simulation was coupled with the rotor dynamics, so that the effects induced by the gusts on important parameters of the problem such as the load coefficient were observable.

Besides the challenge of satisfying the global mass conservation mentioned above, the *far-field boundary conditions* possess two more disadvantages. First, the gust imposed at the inlet has to travel through the entire computational domain before reaching the region of interest, i.e., the building or the wind turbine. Since the grid in the far field is often rather coarse in order to save grid points for the region of interest, the effect of numerical damping of the gust during the approaching phase can be significant. This negative effect can be drastically reduced by a gust-transport mesh, i.e., a secondary fine mesh moving with the gust within the domain using the overset-grid technique (Heinrich, 2014), which is, however, rather costly and not straightforward to implement in general-purpose CFD codes. Second, the entire time interval of this pre-lude has to be simulated, which in case of a fine grid applied to reduce the damping effect can be very CPU-time consuming.

In order to avoid these drawbacks a low-cost approach called field velocity method (FVM) was proposed in the context of aeronautics. Note that wind gusts are also of interest for the aircraft design, since these events can trigger devastating phenomena such as flutter and stall

(Radespiel et al., 2013). The concept of FVM was originally introduced into an unsteady Euler solver by Singh and Baeder (1997) and Parameswaran and Baeder (1997). The objective was to compute the step response of a wing to a sudden change in the angle of attack. The underlying idea of FVM is to alter the convective fluxes by superposition of a pre-defined additional velocity varying in space and time following a gust profile. Since this additional gust velocity is prescribed, the gust can be arranged in the computational domain such that it hits the structure at the first time step of the simulation, avoiding additional CPU costs for transporting the gust from the inlet to the zone of interest. However, the fluid flow and the structure cannot interact with the prescribed wind gust and alter it. Thus, the feedback effect is not taken into account. Despite this important drawback the field velocity method is widely used, see, e.g., Weishäupl and Laschka (2001), Ghoreyshi et al. (2018a, b), Heinrich (2014), Förster and Breitsamter (2015) and Kelleners and Heinrich (2016), where in the last three references it is denoted disturbance velocity approach (DVA).

URANS results obtained by FVM/DVA were compared with the results obtained by the injection of the gust at the inlet studying 2D (airfoils) and 3D geometries (realistic aircraft shapes) in Ghoreyshi et al. (2018a,b); Heinrich (2014); Kelleners and Heinrich (2016). For wind events with short durations deviations of integral quantities such as the lift force appear, whereas the method delivers similar results as long as the gust duration is long enough.

In order to improve the field velocity method (or DVA) Wales et al. (2014) rearranged the Euler equations similar to the field velocity method, but without the assumption that the airfoil has no feedback effect on the gust. The idea lead to the so-called split velocity method (SVM), which introduces additional gust related source terms into the Euler equations. These gust related source terms correspond to the effect of the body on the gust. The SVM was systematically compared with the FVM and with a simulation injecting the gust at the inlet boundary. The SVM (similar to FVM) offers the advantage of reduced CPU costs in comparison to the injection of the gusts at the inlet. However, compared to FVM the split velocity method leads to more reasonable results for gusts with shorter wavelengths.

The objective of the present study is to develop a methodology to efficiently inject strong (turbulent) wind gusts into the flow domain with the intention to be applied later on to fluid-structure interaction simulations. The intended methodology should generate individual wind gusts or a series of such events by an injection within the computational domain and close to the object of interest. Similar to the formulation used for injecting artificial turbulent fluctuations (Breuer, 2018; De Nayer et al., 2018b; Schmidt and Breuer, 2017; Wood et al., 2016) it is based on a source-term formulation and especially considers the feedback effect of the structure and the resulting turbulent flow on the wind gust itself, which is of major importance for highly turbulent flows and their prediction via LES.

Note that the present paper does not intend to predict realistic wind gusts typical for atmospheric conditions. That is beyond the scope of the contribution focusing on the injection methodology and thus considering solely the injection of gusts into either laminar or turbulent free flows. The application towards atmospheric boundary layer flows is intended for subsequent investigations.

The paper is organized as follows: Section 2 describes the general concept and the derivation of the source-term formulation. In Section 3 different shape functions for the wind gusts are presented and a slightly modified C²-1-cosine shape is proposed. The finite-volume flow solver which is used in the present study is briefly described in Section 4. The test cases can be found in Section 5 followed by the results and conclusion sections.

## 2. Momentum source-term formulation for a gust

The general integral momentum conservation equation for an incompressible fluid is considered in the Cartesian basis $\mathcal{B}_0 = (O, e_x, e_y, e_z)$:
$$ \int_V \frac{\partial(\rho \mathbf{u})}{\partial t} dV + \int_V (\mathbf{u} \cdot \nabla)(\rho \mathbf{u}) dV + \int_V \nabla \times [\mu(\nabla \times \mathbf{u})] dV + \int_V \nabla p dV = \mathbf{S} \quad (1) $$

where $\mathbf{u}$ is the velocity vector and $p$ the pressure. $\rho$ denotes the fluid density, $\mu$ the dynamic viscosity and $\mathbf{S}$ the source term not depending on gusts, e.g., the gravity force. $V$ is the volume of the cell. $t$ denotes the physical time.

A new momentum source term $\mathbf{S}^{\text{gust}}$ (with unit: kg m s⁻²) is added to the previous conservation equation to model the generation of a gust inside the computational domain. The gust induced by this source term is only driven by momentum. No additional mass is injected into the domain. Thus, the mass conservation equation remains unchanged and does not include an additional source term. Accordingly, Eq. (1) is simply extended to:

$$ \underbrace{\int_V \frac{\partial(\rho \mathbf{u})}{\partial t} dV}_{\text{Variation in time}} + \underbrace{\int_V (\mathbf{u} \cdot \nabla)(\rho \mathbf{u}) dV}_{\text{Convection}} + \underbrace{\int_V \nabla \times [\mu(\nabla \times \mathbf{u})] dV}_{\text{Diffusion}} + \underbrace{\int_V \nabla p dV}_{\text{Pressure force}} = \mathbf{S} + \mathbf{S}^{\text{gust}} \quad (2) $$

First, this gust momentum source term is formulated in the local basis $\mathcal{B}_1 = (O, \mathbf{e}_\xi, \mathbf{e}_\eta, \mathbf{e}_\zeta)$. The direction of the gust is given by the first vector of the local basis, i.e., $\mathbf{e}_\xi$. This is a user-defined parameter. $\mathbf{e}_\eta$ and $\mathbf{e}_\zeta$ are defined in such a manner that $\mathcal{B}_1$ forms an orthonormal basis. Since convection in a gust plays the dominant role, the gust momentum source term is expressed analog to the convection term of Eq. (2), but with the total velocity of the gust denoted $\mathbf{u}_{\text{gust}}$ instead of the velocity of the fluid without the gust. That means:

$$ \mathbf{S}^{\text{gust}}|_{\mathcal{B}_1} = \int_V (\mathbf{u}^{\text{gust}}|_{\mathcal{B}_1} \cdot \nabla) (\rho \mathbf{u}^{\text{gust}}|_{\mathcal{B}_1}) dV. \quad (3) $$

The vector of the total gust velocity $\mathbf{u}_{\text{gust}}$ is defined as follows:

$$ \mathbf{u}^{\text{gust}} = \mathbf{u}_{\text{base}} + \mathbf{u}_{\text{deterministic}}^{\text{gust}}. \quad (4) $$

$\mathbf{u}_{\text{base}}$ is the vector corresponding to the velocity of the surrounding fluid without the gust. Using the explicit time-marching scheme described in Section 4, this velocity is given by the local velocity of the previous time step $\mathbf{u}^n$. $\mathbf{u}_{\text{deterministic}}^{\text{gust}}$ is the vector, which only contains the velocity profile of the gust without the surrounding fluid. This gust shape is deterministically defined in the local basis $\mathcal{B}_1 = (O, \mathbf{e}_\xi, \mathbf{e}_\eta, \mathbf{e}_\zeta)$ depicted in Fig. 1. Its definition is based on a user-defined amplitude of the gust $A_g$ and user-defined analytic functions representing the shape of the gust in time and space. As mentioned above, the direction of the gust is the first vector of the local basis $\mathbf{e}_\zeta$. Thus, the velocity of the gust in the local basis reads:

$$ \mathbf{u}_{\text{deterministic}}^{\text{gust}}|_{\mathcal{B}_1}(t, \xi, \eta, \zeta) = A_g f_1(t) f_\xi(\xi) f_\eta(\eta) f_\zeta(\zeta) \mathbf{e}_\xi. \quad (5) $$

The typically used functions for the spatial $f_\xi(\xi)$, $f_\eta(\eta)$, $f_\zeta(\zeta)$ as well as the temporal $f_1(t)$ distributions related to a gust are described in Section 3.

The momentum source term (3) can now be simplified. Since the

<figure><img src="image_4.png" /><figcaption>Fig. 1. Definition of the gust in the local basis $\mathcal{B}_1 = (O, \mathbf{e}_\xi, \mathbf{e}_\eta, \mathbf{e}_\zeta)$.</figcaption></figure>

second and third components of $\mathbf{u}^{\text{gust}}$ written in $\mathcal{B}_1$, $u_\eta^{\text{gust}}|_{\mathcal{B}_1}$ and $u_\xi^{\text{gust}}|_{\mathcal{B}_1}$, are negligible compared with its first component $u_\xi^{\text{gust}}|_{\mathcal{B}_1}$ (the $\eta$ and $\zeta$-components are solely composed of the components of $\mathbf{u}_{\text{base}}$ in $\mathcal{B}_1$), only the first component of $(\mathbf{u}^{\text{gust}}|_{\mathcal{B}_1} \cdot \nabla)$ remains. Thus, Eq. (3) reads:

$$ \mathbf{S}^{\text{gust}}|_{\mathcal{B}_1} \approx \int_V u_\xi^{\text{gust}}|_{\mathcal{B}_1} \frac{\partial}{\partial \xi} (\rho \mathbf{u}^{\text{gust}}|_{\mathcal{B}_1}) dV. \quad (6) $$

Furthermore, it is assumed that the gradient in the gust direction of the velocity without gust $\partial \mathbf{u}_{\text{base}}|_{\mathcal{B}_1}/\partial \xi$ is negligible compared with $\partial \mathbf{u}_{\text{deterministic}}^{\text{gust}}|_{\mathcal{B}_1}/\partial \xi$. Accordingly, Eq. (6) can be approximated:

$$ \begin{aligned} \mathbf{S}^{\text{gust}}|_{\mathcal{B}_1} &\approx \int_V u_\xi^{\text{gust}}|_{\mathcal{B}_1} \left( \frac{\partial}{\partial \xi} (\rho \mathbf{u}_{\text{base}}|_{\mathcal{B}_1}) + \frac{\partial}{\partial \xi} (\rho \mathbf{u}_{\text{deterministic}}^{\text{gust}}|_{\mathcal{B}_1}) \right) dV \\ &\approx \int_V u_\xi^{\text{gust}}|_{\mathcal{B}_1} \frac{\partial}{\partial \xi} (\rho \mathbf{u}_{\text{deterministic}}^{\text{gust}}|_{\mathcal{B}_1}) dV. \end{aligned} \quad (7) $$

As usual in a finite-volume discretization, the midpoint rule is applied to approximate the volume integral:

$$ \begin{aligned} \mathbf{S}^{\text{gust}}|_{\mathcal{B}_1} &\approx u_\xi^{\text{gust}}|_{\mathcal{B}_1, \text{cc}} \frac{\partial}{\partial \xi} (\rho \mathbf{u}_{\text{deterministic}}^{\text{gust}}|_{\mathcal{B}_1})|_{\text{cc}} V \\ &\approx u_\xi^{\text{gust}}|_{\mathcal{B}_1, \text{cc}} \frac{\partial}{\partial \xi} (\rho \mathbf{u}_{\text{deterministic}}^{\text{gust}}, \xi|_{\mathcal{B}_1} \mathbf{e}_\xi)|_{\text{cc}} V. \end{aligned} \quad (8) $$

Here the expression cc abbreviates the notation cell center. The vector $\mathbf{e}_\xi$ present under the derivative is constant and can be extracted from this derivative. Combining it with $u_\xi^{\text{gust}}|_{\mathcal{B}_1, \text{cc}}$, the velocity $\mathbf{u}^{\text{gust}}|_{\mathcal{B}_1, \text{cc}}$ appears again and the expression reads:

$$ \mathbf{S}^{\text{gust}}|_{\mathcal{B}_1} \approx \rho \underbrace{\frac{\partial}{\partial \xi} (u_{\text{deterministic}}^{\text{gust}}, \xi|_{\mathcal{B}_1})|_{\text{cc}}}_{\text{analog to a mass flow rate}} V \mathbf{u}^{\text{gust}}|_{\mathcal{B}_1, \text{cc}}. \quad (9) $$

Using the Gauss theorem Eq. (9) can be interpreted as a mass flow rate in $\xi$-direction multiplied by the total gust velocity.

$\mathbf{S}^{\text{gust}}|_{\mathcal{B}_1}$ has now to be written in the Cartesian basis $\mathcal{B}_0 = (O, \mathbf{e}_x, \mathbf{e}_y, \mathbf{e}_z)$. The vector $\mathbf{e}_\xi$ is a user-defined parameter in $\mathcal{B}_0$. Since it defines the direction of the gust, it will be denoted $\mathbf{n}^\xi$ in the following. The transfer matrix $\mathbf{T}_{\mathcal{B}_0 \to \mathcal{B}_1}$ from basis $\mathcal{B}_0$ to basis $\mathcal{B}_1$ is known:

$$ \mathbf{T}_{\mathcal{B}_0 \to \mathcal{B}_1} = \begin{bmatrix} \xi_x & \eta_x & \zeta_x \\ \xi_y & \eta_y & \zeta_y \\ \xi_z & \eta_z & \zeta_z \end{bmatrix}, \quad (10) $$

where $\xi_x, \xi_y$ and $\xi_z$ are the $x$, $y$- and $z$-components of the vector $\mathbf{e}_\xi$. The same notations are used for the components of the vectors $\mathbf{e}_\eta$ and $\mathbf{e}_\zeta$. Since $\mathcal{B}_0$ and $\mathcal{B}_1$ are orthonormal bases, $\mathbf{T}_{\mathcal{B}_0 \to \mathcal{B}_1}$ is an orthogonal matrix. Thus, the transfer matrix $\mathbf{T}_{\mathcal{B}_1 \to \mathcal{B}_0}$ from basis $\mathcal{B}_1$ to basis $\mathcal{B}_0$ can be easily deduced as the transposed matrix of $\mathbf{T}_{\mathcal{B}_0 \to \mathcal{B}_1}$:

$$ \mathbf{T}_{\mathcal{B}_1 \to \mathcal{B}_0} = \mathbf{T}_{\mathcal{B}_0 \to \mathcal{B}_1}^{-1} = \mathbf{T}_{\mathcal{B}_0 \to \mathcal{B}_1}^{\text{T}} = \begin{bmatrix} \xi_x & \xi_y & \xi_z \\ \eta_x & \eta_y & \eta_z \\ \xi_x & \xi_y & \xi_z \end{bmatrix}. \quad (11) $$

The local coordinates $\xi$, $\eta$ and $\zeta$ can be also expressed depending on the Cartesian coordinates:

$$ \xi = t_\xi(x, y, z) = \xi_x x + \xi_y y + \xi_z z, $$

$$ \eta = t_\eta(x, y, z) = \eta_x x + \eta_y y + \eta_z z, $$

$$ \zeta = t_\zeta(x, y, z) = \zeta_x x + \zeta_y y + \zeta_z z. $$

It leads to the following expression for the gust velocity:

$$ \mathbf{u}_{\text{deterministic}}^{\text{gust}}|_{\mathcal{B}_0}(t, x, y, z) = A_g f_1(t) f_\xi(t_\xi(x, y, z)) f_\eta(t_\eta(x, y, z)) f_\zeta(t_\zeta(x, y, z)) \mathbf{n}^\xi|_{\mathcal{B}_0}. \quad (12) $$
<figure><img src="image_5.png" /><figcaption>Fig. 2. Typical shapes for the deterministic gust models for a given scale $L_g^\phi$ and a given central value $\phi_g$.</figcaption></figure>

Consequently, the source-term formulation in the Cartesian basis $\mathcal{B}_0$ reads:

$$ S_{\mathcal{B}_0}^{\text{gust}} \approx \rho A_g f_t(t) f_\eta(t_\eta(x_{cc}, y_{cc}, z_{cc})) f_\zeta(t_\zeta(x_{cc}, y_{cc}, z_{cc})) \left. \frac{\partial f_\zeta}{\partial \zeta} \right|_{cc} V \mathbf{u}_{\mathcal{B}_0,cc}^{\text{gust}} . \quad (13) $$

To use the same terminology as in the source-term formulation of Pasquetti and Peres (2015), who injected a micro-jet inside a fluid domain, an expulsion phase and an ingestion phase have to be distinguished during the gust injection. The gust is only formed during the expulsion phase when the velocity rises. The ingestion phase does not beneficially contribute to the generation of the gust. Thus, in order to consider only the expulsion phase in the present gust source term, the idea has to be transformed from the temporal to the spatial coordinate. Taking one of the gust shapes defined below in Section 3 for $f_\zeta$, only the left part of these functions (based on the notation of Section 3 that means $\zeta < \xi_g$, where $\xi_g$ is the center of the gust injection in $\xi$-direction) leads to a positive mass flow rate as denoted in Eq. (9) and thus to a beneficial source term supporting the gust generation process. In summary, the final

### 3. Gust shapes

In this section the typically used deterministic functions required by Eq. (14) for the spatial and temporal distributions related to a gust are described. Thus, the variable $\phi$ might be the coordinates $(\phi = \{\xi, \eta, \zeta\})$ or the time $t$. The constant $\phi_g$ is the central value of the gust distribution for the corresponding $\phi$ and $L_g^\phi$ its length or time scale. As mentioned in the introduction the Extreme Coherent Gust (ECG) and the Extreme Operating Gust (EOG) are defined in the IEC-Standard (2002). Another option is a Gaussian distribution.

* Adapted “1-cosine” shape (Extreme Coherent Gust, ECG):

The original “1-cosine” shape found in the literature (IEC-Standard, 2002) is expressed as follows:

$$ S_{\mathcal{B}_0}^{\text{gust}} \approx \begin{cases} \rho A_g f_t(t) f_\eta(t_\eta(x_{cc}, y_{cc}, z_{cc})) f_\zeta(t_\zeta(x_{cc}, y_{cc}, z_{cc})) \left. \frac{\partial f_\zeta}{\partial \zeta} \right|_{cc} V \mathbf{u}_{\mathcal{B}_0,cc}^{\text{gust}} & \text{for } \zeta < \xi_g \\ 0 & \text{else} \end{cases} \quad (14) $$

source-term formulation reads:

Regarding this new source-term formulation the following features have to be clearly mentioned:

* Contrary to the micro-jet model by Pasquetti and Peres (2015) using a punctual source, the present gust model relies on a volumetric distribution of the source term.

* Since the underlying physics of micro-jets and gusts significantly differ, the model for the micro-jets relies on a momentum and an additional mass source term, whereas gusts can be generated without a mass source term. Although in lab-scale experiments gusts are often generated by a supply of an additional mass flux (see, e.g., Volpe et al., 2013; Tischmacher and Ruck, 2013; Gromke and Ruck, 2018), in nature that is not the case. Instead the existing fluid of the surrounding is locally accelerated. Consequently, the pure source-term formulation for the momentum equation as suggested here is the more natural choice in case of gusts.

* Owing to the fact that no additional mass flux is included in the source-term formulation, it is not necessary to relax the divergence-free constraint in the region of the injection as done in the case of the micro-jets (Pasquetti and Peres, 2015).

$$ f_\phi(\phi) = \begin{cases} \frac{1}{2} \left( 1 - \cos \left( \frac{\pi \phi}{L_g^\phi} \right) \right) & \text{for } \phi \in [0, L_g^\phi] \\ 0 & \text{else} \end{cases} \quad (15) $$

In order to achieve more control, the original shape is adapted introducing the central value $\phi_g$ (De Nayer et al., 2019):

$$ f_\phi(\phi) = \begin{cases} \frac{1}{2} \left( 1 + \cos \left( \frac{2\pi(\phi - \phi_g)}{L_g^\phi} \right) \right) & \text{for } (\phi - \phi_g) \in \left[ -\frac{L_g^\phi}{2}, \frac{L_g^\phi}{2} \right] \\ 0 & \text{else} \end{cases} \quad (16) $$

* Adapted “Mexican hat” shape (Extreme Operating Gust, EOG):

In the same manner the original “Mexican hat” shape (IEC-Standard, 2002) is adapted introducing $\phi_g$ (De Nayer et al., 2019):

$$ f_\phi(\phi) = \begin{cases} 0.37\cos\left(\frac{3\pi(\phi - \phi_g)}{L_g^\phi}\right) \left(1 + \cos\left(\frac{2\pi(\phi - \phi_g)}{L_g^\phi}\right)\right) & \text{for } (\phi - \phi_g) \in \left[-\frac{L_g^\phi}{2}, \frac{L_g^\phi}{2}\right] \\ 0 & \text{else} \end{cases} \quad (17) $$
<figure><img src="image_6.png" /><figcaption>Fig. 3. Derivatives of the three shape functions for the deterministic gust models for a given scale $L_g^\phi$ and a given central value $\phi_g$.</figcaption></figure>

### Gaussian distribution:

$$f_{\phi}(\phi) = \exp\left(-18\left(\frac{\phi - \phi_g}{L_g^{\phi}}\right)^2\right). \qquad (18)$$

This Gaussian distribution is similar to the distribution used in the digital filter method applied to generate synthetic turbulent velocity components (Klein et al., 2003). The coefficient 18 is chosen so that the value of the function is less than 0.01% of the peak value at $\phi = -L_g^\phi/2$ and $\phi = L_g^\phi/2$.

Fig. 2 illustrates the three typical distributions for a given scale $L_g^\phi$ and a given central value $\phi_g$.

Since the source-term formulation proposed in Eq. (14) also relies on the derivative $\partial f_\xi / \partial \xi$ of the shape function, the three corresponding derivatives are given and plotted in Fig. 3.

### Adapted “1-cosine” shape (Extreme Coherent Gust, ECG):

$$\frac{\partial f_{\phi}(\phi)}{\partial \phi} = \begin{cases} \frac{\pi}{L_g^{\phi}} \sin\left(\frac{2\pi(\phi - \phi_g)}{L_g^{\phi}}\right) & \text{for } (\phi - \phi_g) \in \left[-\frac{L_g^{\phi}}{2}, \frac{L_g^{\phi}}{2}\right] \\ 0 & \text{else.} \end{cases} \qquad (19)$$

### Adapted “Mexican hat” shape (Extreme Operating Gust, EOG):

$$\frac{\partial f_{\phi}(\phi)}{\partial \phi} = \begin{cases} \frac{0.37 \pi}{L_g^{\phi}} \left[ 3 \sin\left(\frac{3\pi(\phi - \phi_g)}{L_g^{\phi}}\right) \left(1 + \cos\left(\frac{2\pi(\phi - \phi_g)}{L_g^{\phi}}\right)\right) + 2 \cos\left(\frac{3\pi(\phi - \phi_g)}{L_g^{\phi}}\right) \sin\left(\frac{2\pi(\phi - \phi_g)}{L_g^{\phi}}\right) \right] & \text{for } (\phi - \phi_g) \in \left[-\frac{L_g^{\phi}}{2}, \frac{L_g^{\phi}}{2}\right] \\ 0 & \text{else.} \end{cases} \qquad (20)$$

### Gaussian distribution:

$$\frac{\partial f_{\phi}(\phi)}{\partial \phi} = - \frac{36}{\left(L_g^{\phi}\right)^2} (\phi - \phi_g) \exp\left(-18\left(\frac{\phi - \phi_g}{L_g^{\phi}}\right)^2\right). \qquad (21)$$

In the current source-term formulation $S^{\text{gust}}|_{\partial\Omega_0}$ has non-zero values only for $\xi < \xi_g$ (see Eq. (14)). Applying this cut on the above presented derivatives leads to a kink in the derivative of each function at $\xi = \xi_g$ or more generally at $\phi = \phi_g$. In the CFD simulations this discontinuity may produce spurious numerical oscillations during the introduction of the gust. In order to resolve this problem, the function $\partial f_\phi / \partial \phi$ should possess a continuous derivative at $\phi = \phi_g$. For this purpose, such a function with a continuous second derivative, denoted ECG-$C^2$, is derived from the classic ECG:

<figure><img src="image_7.png" /><figcaption>Fig. 4. Proposed $C^2$-"1-cosine" gust shape and its derivative for a given scale $L_g^\phi$ and a given central value $\phi_g$.</figcaption></figure>
<table><caption>**Table 1**
Gust configuration for the reference simulation.</caption><thead><tr><th>n<sup>B</sup>|<sub>B<sub>0</sub></sub></th><th>f<sub>i</sub></th><th>f<sub>z</sub></th><th>f<sub>η</sub></th><th>f<sub>ζ</sub></th><th>L<sub>t</sub> u<sub>∞</sub> / L</th><th>L<sub>z</sub> / L</th><th>L<sub>η</sub> / L</th><th>L<sub>ζ</sub> / L</th><th>A<sub>R</sub> / u<sub>∞</sub></th></tr></thead><tbody><tr><td>e<sub>x</sub></td><td>ECG-C<sup>2</sup></td><td>ECG-C<sup>2</sup></td><td>ECG</td><td>ECG</td><td>0.5</td><td>0.5</td><td>0.25</td><td>0.25</td><td>1.0</td></tr></tbody></table>

At first, the derivative of the requested function is defined based on the Extreme Coherent Gust distribution:

The surface and volume integrals are approximated by the midpoint rule. Most flow variables are linearly interpolated to the cell faces leading to a

$$ \frac{\partial f_{\phi}(\phi)}{\partial \phi} = \begin{cases} \frac{2}{L_g} \left( 1 + \cos\left(\frac{4\pi(\phi - \phi_g)}{L_g^{\phi}} + \pi\right) \right) & \text{for } \left(\phi - \left(\phi_g - \frac{L_g^{\phi}}{4}\right)\right) \in \left[-\frac{L_g^{\phi}}{4}, \frac{L_g^{\phi}}{4}\right] \\ -\frac{2}{L_g} \left( 1 + \cos\left(\frac{4\pi(\phi - \phi_g)}{L_g^{\phi}} - \pi\right) \right) & \text{for } \left(\phi - \left(\phi_g + \frac{L_g^{\phi}}{4}\right)\right) \in \left[-\frac{L_g^{\phi}}{4}, \frac{L_g^{\phi}}{4}\right] \\ 0 & \text{else.} \end{cases} \quad (22) $$

After integration the function ECG-C² reads:

* Modified C²-"1-cosine" shape (Extreme Coherent Gust C², ECG-C²):

$$ f_{\phi}(\phi) = \begin{cases} \frac{2}{L_g^{\phi}} \left( \phi - \left( \phi_g - \frac{L_g^{\phi}}{2} \right) + \frac{L_g^{\phi}}{4\pi} \sin \left( \frac{4\pi(\phi - \phi_g)}{L_g^{\phi}} + \pi \right) \right) & \text{for } \left( \phi - \left( \phi_g - \frac{L_g^{\phi}}{4} \right) \right) \in \left[ -\frac{L_g^{\phi}}{4}, \frac{L_g^{\phi}}{4} \right] \\ -\frac{2}{L_g^{\phi}} \left( \phi - \left( \phi_g + \frac{L_g^{\phi}}{2} \right) + \frac{L_g^{\phi}}{4\pi} \sin \left( \frac{4\pi(\phi - \phi_g)}{L_g^{\phi}} - \pi \right) \right) & \text{for } \left( \phi - \left( \phi_g + \frac{L_g^{\phi}}{4} \right) \right) \in \left[ -\frac{L_g^{\phi}}{4}, \frac{L_g^{\phi}}{4} \right] \\ 0 & \text{else.} \end{cases} \quad (23) $$

Fig. 4 shows the ECG-C² distribution and its derivative.

## 4. Flow solver

As mentioned above the Navier-Stokes equations for an incompressible fluid are solved based on a finite-volume method, i.e., an enhanced version of FASTEST-3D (Breuer et al., 2012; Durst and Schäfer, 1996) and corresponding validation studies (De Nayer et al., 2014, 2018a, 2020; De Nayer and Breuer, 2014). The equations are discretized on a curvilinear, block-structured body-fitted grid with a collocated variable arrangement.

second-order accurate central scheme. The convective fluxes are approximated by the technique of flux blending (Ferziger and Perić, 2002; Khosla and Rubin, 1974) to stabilize the simulation. For the current case the flux blending includes 3% of a first-order accurate upwind scheme and 97% of a second-order accurate central scheme. This choice is motivated as follows: Based on the standard grid defined in Section 5 and the reference simulation explained in Section 6.2.1 (Table 1) the local Péclet number Pe<sub>Δx</sub> = ρ u Δx/μ predicted with the grid spacing Δx is much larger than 2. Accordingly, a pure central scheme is particularly susceptible to 2Δx oscillations. The flux blending applied avoids this problem. The subdivision (3%/97%) relies on extensive tests based on the reference case, which on the one hand have shown that below 3% upwind oscillations cannot be avoided. On the other hand, higher portions lead to an increase in numerical diffusion, which can be avoided. Furthermore, the momentum interpolation technique of Rhie and Chow (1983) is applied to couple the pressure and the velocity fields on non-staggered grids.

A semi-implicit predictor-corrector scheme (Breuer et al., 2012) (projection method of second-order accuracy) is used to solve the pressure-velocity coupling problem, which is appropriate for the large-eddy simulation (LES) technique. Here a low-storage multi-stage Runge-Kutta method is applied for time-marching the momentum equations within the predictor-corrector scheme. The subsequent corrector step ensures that mass conservation is achieved in form of a divergence-free velocity field. For this purpose, a Poisson equation for the pressure correction is solved by an incomplete LU decomposition method (Stone, 1968). Overall, that provides second-order accuracy in space and time.

The prediction of turbulent flows relies on LES, where the large scales

<figure><img src="image_8.png" /><figcaption>Fig. 5. Geometry and boundary conditions of the test case.</figcaption></figure>
of the turbulent flow are resolved directly and the small scales are modeled by a subgrid-scale (SGS) model. For this purpose, several models are implemented in FASTEST-3D such as the classical Smagorinsky model (Smagorinsky, 1963) combined with the Van-Driest damping function near rigid walls, the dynamic Smagorinsky model (Germano et al., 1991) and the WALE model (Nicoud and Ducros, 1999). Owing to minor influences of the SGS model considered for the present applications, solely the standard Smagorinsky model with the classical model parameter $C_s = 0.1$ is applied.

For LES an appropriate representation of the inflow boundary conditions is often necessary (Wood et al., 2016; Breuer, 2018). In FASTEST-3D an inflow generator based on the digital filter concept of Klein et al. (2003) mimics incoming turbulent perturbations. The original concept was extended to permit the superposition of these turbulent perturbations as source terms inside the domain (Breuer, 2018; De Nayer et al., 2018b; Schmidt and Breuer, 2017).

## 5. Test cases

In order to test the proposed source-term formulation for injecting gusts, two test cases are chosen. First, a gust is injected into an undisturbed flow and its evolution in space and time is investigated. Based on this test case the first goal of this study can be evaluated in detail, i.e., the injection of a gust within the computational domain of a simulation. Second, a gust is inserted into a turbulent flow to address the second goal, i.e., to simulate gusts and their propagation within turbulent flows.

Note that the objective here is not to predict realistic wind gusts typical for atmospheric conditions as already motivated in the introduction. Thus, neither an atmospheric boundary layer flow is considered, nor specific characteristics of wind gusts under atmospheric conditions such as their coherence in the vertical direction (Branlard, 2009) or vertical profiles of the streamwise velocity across the gust events (Brasseur, 2001) are taken into account.

The geometry of the free flow considered is a cuboid (length × height × depth = $4L \times L \times L$) as sketched in Fig. 5. Without loss of generality the west face is assumed to be the inlet, whereas the east face is the outlet. The other faces are periodic in pairs. The grid is block-structured and contains on the standard level $512 \times 128 \times 128$ cells, which are equidistant in all directions. In order to assure the quality of the results, a grid independence study is carried out in Section 6.1. For parallelization based on domain decomposition the grid is decomposed into 128 blocks and the simulations run on 128 processors.

Unless otherwise indicated, the following gust shape is assumed as the standard case: In η- and ζ-direction the ECG shape function is chosen. Since the derivative of the ξ-function with respect to ξ appears in the final source-term formulation, the ECG is replaced by the ECG-C² shape in ξ-direction as already discussed in Section 3. Assuming that the ξ-direction and the time are linked by the Taylor hypothesis, the ECG-C² shape is also applied for the temporal distribution. Furthermore, the length in ξ-direction $L_ξ$ can be computed based on $L_t$:

$$L_{\xi} = L_t u_{\text{conv}}, \qquad (24)$$

where $u_{\text{conv}}$ is the convective velocity of the flow field.

The direction of the gust is set to $\mathbf{n}^g|_{\partial\Omega_0} = \mathbf{e}_x$ and its center is located at $(\xi_g = x_g = L/2, \eta_g = y_g = 0, \zeta_g = z_g = 0)$ as shown in Fig. 5. In order to carry out parameter studies, the time and length scales of the gust vary in the range 0.25–1. The different values investigated are listed in Tables 2 and 3 in Section 6.2. The background flow without the gust is also directed in x-direction. A constant inlet velocity is set to $\mathbf{u}|_{\text{inlet}} = \mathbf{u}_{\infty} \mathbf{e}_x$ with $u_{\infty} = 1 \text{ m/s}$, which leads to a convective velocity of $u_{\text{conv}} = u_{\infty}$. The time step of the simulation is chosen to $\Delta t = 4 \times 10^{-3}$ s, which guarantees a CFL number less than unity. The density of the fluid is set to $\rho = 1 \text{ kg m}^{-3}$ and the dynamic viscosity to $\mu = 10^{-6} \text{ kg m}^{-1} \text{ s}^{-1}$ leading to a Reynolds number of $Re = 10^6$ based on L.

For the turbulent test case the background flow is superimposed by

<caption>**Table 2**</caption>

Parameter study with varying time scales $L_t$ vs. varying gust amplitudes $A_g$ ($L_η/L = 0.25, L_ζ/L = 0.25$): Time-averaged maximal streamwise velocity.

<table><thead><tr><th>L<sub>t</sub> u<sub>∞</sub>/L</th><th>A<sub>g</sub> u<sub>∞</sub></th><th>u<sub>max</sub><sup>theory</sup> u<sub>∞</sub></th><th>u<sub>max</sub><sup>sim</sup> u<sub>∞</sub></th><th>Δu u<sub>∞</sub></th><th>u<sub>max</sub><sup>sim</sup> u<sub>max</sub> u<sub>∞</sub></th><th>ε = Δu/u<sub>max</sub><sup>theory</sup></th></tr></thead><tbody><tr><td>0.25</td><td>0.5</td><td>1.5</td><td>1.34</td><td>-0.16</td><td></td><td>-11%</td></tr><tr><td>0.25</td><td>1.0</td><td>2.0</td><td>1.77</td><td>-0.23</td><td></td><td>-12%</td></tr><tr><td>0.25</td><td>2.0</td><td>3.0</td><td>2.54</td><td>-0.46</td><td></td><td>-15%</td></tr><tr><td>0.5</td><td>0.5</td><td>1.5</td><td>1.48</td><td>-0.02</td><td></td><td>-1%</td></tr><tr><td>0.5</td><td>1.0</td><td>2.0</td><td>2.07</td><td>0.07</td><td></td><td>+4%</td></tr><tr><td>0.5</td><td>2.0</td><td>3.0</td><td>2.87</td><td>-0.13</td><td></td><td>-4%</td></tr><tr><td>1.0</td><td>0.5</td><td>1.5</td><td>1.56</td><td>0.06</td><td></td><td>+4%</td></tr><tr><td>1.0</td><td>1.0</td><td>2.0</td><td>2.29</td><td>0.29</td><td></td><td>+15%</td></tr><tr><td>1.0</td><td>2.0</td><td>3.0</td><td>3.30</td><td>0.30</td><td></td><td>+10%</td></tr></tbody></table>

homogeneous isotropic turbulence based on the methodology briefly described in Section 4. For this purpose, the digital filter concept of Klein et al. (2003) requires the following input parameters in order to generate turbulent perturbations. Solely the level of the Reynolds stresses and the definition of one integral time scale ($T^{\text{turb}}$) and two integral length scales ($L_y^{\text{turb}}$, $L_z^{\text{turb}}$) are necessary to generate artificial turbulence with proper autocorrelations in time and two-point correlations in space. Owing to the simplifying assumption of isotropy of the approaching flow, which is the usual way for describing the background turbulence found in experimental facilities such as water or wind tunnels, the Reynolds shear stresses are set to zero and all three normal Reynolds stresses are equal $(\overline{u}\overline{u}' = \overline{v}\overline{v}' = \overline{w}\overline{w}'')$.

The constant value $\overline{u}\overline{u}'$ chosen for the entire injection plane located upstream of the gust injection at about $x/L \approx 0.06$ (see Fig. 5) depends on the intended turbulence intensity $TI = \sqrt{\overline{u}\overline{u}'}/u_{\infty}$ chosen. It is worth to note that in principle the synthetically generated turbulence can also be applied as inflow boundary conditions on the west face of the domain. However, here the source-term formulation for the turbulence injection is favored in order to demonstrate that it can be effectively combined with the source-term formulation of the gust. In the present study two different turbulence intensities are investigated, i.e., $TI = 4\%$ and $8\%$. In dimensionless form the integral time scale is assumed to be given by $T^{\text{turb}} u_{\infty}/L = 0.05$. With the aforementioned assumption of isotropy and the Taylor hypothesis the two missing integral length scales are $L_y^{\text{turb}}/L = L_z^{\text{turb}}/L = 0.05$.

## 6. Results

### 6.1. Grid independence study

In order to assure the quality of the results presented in the subsequent sections, a study on the effect of the grid resolution is carried out first. Starting from the standard grid (denoted grid S) with $512 \times 128 \times 128$ cells two coarser and one finer resolutions are taken into account. For this purpose, the number of cells in all directions are halved or

<caption>**Table 3**</caption>

Parameter study with varying length scales $L_η$ vs. varying gust amplitudes $A_g$ ($L_t u_{\infty}/L = 0.5, L_ζ/L = 0.25$): Time-averaged maximal streamwise velocity.

<table><thead><tr><th>L<sub>η</sub> L</th><th>A<sub>g</sub> u<sub>∞</sub></th><th>u<sub>max</sub><sup>theory</sup> u<sub>∞</sub></th><th>u<sub>max</sub><sup>sim</sup> u<sub>∞</sub></th><th>Δu u<sub>∞</sub></th><th>u<sub>max</sub><sup>sim</sup> u<sub>max</sub> u<sub>∞</sub></th><th>ε = Δu/u<sub>max</sub><sup>theory</sup></th></tr></thead><tbody><tr><td>0.25</td><td>0.5</td><td>1.5</td><td>1.48</td><td>-0.02</td><td></td><td>-1%</td></tr><tr><td>0.25</td><td>1.0</td><td>2.0</td><td>2.07</td><td>0.07</td><td></td><td>+4%</td></tr><tr><td>0.25</td><td>2.0</td><td>3.0</td><td>2.87</td><td>-0.13</td><td></td><td>-4%</td></tr><tr><td>0.5</td><td>0.5</td><td>1.5</td><td>1.44</td><td>-0.06</td><td></td><td>-4%</td></tr><tr><td>0.5</td><td>1.0</td><td>2.0</td><td>1.98</td><td>-0.02</td><td></td><td>-1%</td></tr><tr><td>0.5</td><td>2.0</td><td>3.0</td><td>2.74</td><td>-0.26</td><td></td><td>-9%</td></tr><tr><td>1.0</td><td>0.5</td><td>1.5</td><td>1.41</td><td>-0.09</td><td></td><td>-6%</td></tr><tr><td>1.0</td><td>1.0</td><td>2.0</td><td>1.86</td><td>-0.14</td><td></td><td>-7%</td></tr><tr><td>1.0</td><td>2.0</td><td>3.0</td><td>2.50</td><td>-0.50</td><td></td><td>-17%</td></tr></tbody></table>
doubled, respectively. Thus, the finest grid *F* consists of 1024 × 256 × 256 cells, whereas the coarse grid *C* and the very coarse grid *VC* include 256 × 64 × 64 and 128 × 32 × 32 cells, respectively. The time-step size is kept constant.

The evaluation is done for a standard gust configuration, whose pa-
rameters are summarized in Table 1. This set of parameters is considered
as the reference configuration in the whole study. Since an instantaneous
flow field is considered, the evaluation is carried out at two specific in-
stants in time, i.e., shortly after the injection process has finished (t *u*<sub>∞</sub>/
*L* = 0.8) and after the gust has traveled a certain distance through the
domain (t *u*<sub>∞</sub>/L = 1.6). Fig. 6 depicts the distribution of the streamwise
velocity *u*/*u*<sub>∞</sub> and the pressure *p*/*(ρ*u*<sub>∞</sub><sup>2</sup>)* along the streamwise coordinate
at the center of the gust.

As visible in orange in Fig. 6(a) the evolution of the streamwise ve-
locity and the pressure obtained on the very coarse grid VC differs

significantly from the others. This difference increases in form and
amplitude with time (see Fig. 6(b)). The coarse grid C delivers better
results: The streamwise velocity and pressure distributions are similar to
those predicted on the finer meshes. However, they quantitatively differ.
Therefore, the grids VC and C are too coarse and will not be used in the
following.

The results obtained on the standard grid *S* and the finest grid *F* are quasi identical in all areas of the curves just after the injection (see Fig. 6(a)). As time progresses a marginal increase of the deviations appears as shown in the close-ups in Fig. 6(b).

In conclusion, the grid independence study shows that the results
converge with increasing grid resolution and that the deviations
observed between the standard grid *S* and the finest grid *F* are marginal.
Consequently, for the subsequent detailed investigations grid *S* is a good
compromise between efficiency and accuracy.

<figure><img src="image_9.png" /><figcaption>**Fig. 6.** Grid independence study: Dimensionless streamwise velocity and pressure after the injection of a gust into an undisturbed flow at two instants in time (Reference configuration of the gust: ECG-C² shape in ξ and t, standard ECG shape in η- and ζ-direction, A_g/u_∞ = 1, L_t u_∞/L = 0.5, L_η/L = L_ζ/L = 0.25).</figcaption></figure>
## 6.2. Propagation of wind gusts in undisturbed flow

### 6.2.1. Reference simulation

In order to get a deeper insight into the physics of the flow during and after the injection of a gust, a simulation is carried out with the standard parameters summarized in Table 1.

Fig. 7 depicts the evolution of the injected gust in time and space within the undisturbed flow. Since the direction of the gust is equal to the x-axis, the streamwise velocity and pressure along the central axis for y/L = z/L = 0 are investigated at certain instants in time. To complete the

previous graphs, a 3D view of the wind gust is realized in Fig. 8. The values of the streamwise velocity and the pressure in the xy-plane at z/L = 0 are considered and plotted as the third axis.

During the injection ($t u_\infty/L \le 0.5$, which corresponds to the first three instants depicted) the flow is accelerated up to $u_{\text{max}}^{\text{sim}}/u_\infty = 2.24$. In this phase the x-distribution of the streamwise velocity is quasi symmetric with respect to $x_g/L = 0.5$. Moreover, comparing the velocity distribution of the third time instant (see Fig. 7(c)) with the theoretical shape of ECG-C² (see Fig. 4), the simulated shape looks very similar to the theoretical one. This is a particularly noteworthy result because it clearly

<figure><img src="image_10.png" /><figcaption>Fig. 7. Dimensionless streamwise velocity (left) and pressure (right) of a gust injected into an undisturbed flow at different instants in time and for y/L = 0 and z/L = 0 (Reference configuration of the gust: ECG-C² shape in ξ and t, standard ECG shape in η- and ζ-direction, $A_g/u_\infty = 1$, $L_t u_\infty/L = 0.5$, $L_\eta/L = L_\zeta/L = 0.25$).</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>**Fig. 8.** Three-dimensional view of the dimensionless streamwise velocity (left) and pressure (right) of a gust injected into an undisturbed flow at different instants in time in the xy-plane at $z/L = 0$ (Reference configuration of the gust: ECG-C² shape in ξ and t, standard ECG shape in η- and ζ-direction, $A_g/u_∞ = 1$, $L_t u_∞/ L = 0.5$, $L_η/ L = L_ζ/L = 0.25$).</figcaption></figure>

shows that the proposed source-term formulation works as desired. As a reaction to the additional momentum injected into the flow, the pressure takes a negative sinusoidal form. During the injection phase the pressure disturbance $\Delta p/(\rho u_\infty^2)$ reaches its maximum with about 0.18. In the rest of the simulation the dimensionless pressure varies below $\pm 0.1$.

After the injection phase ($t u_\infty/L > 0.5$) the gust travels through the space with a quasi constant propagation velocity equal to approximately $1.4 u_\infty$. This propagation velocity is evaluated based on the position of the velocity maximum as a function of time. It is larger than the velocity of the surrounding fluid, but much smaller than the maximum velocity of the gust. Indeed, the structure of the gust is

composed of fluid elements, which move at highly different velocities: In the middle of the gust the magnitude of their velocities approaches $|\mathbf{u}_{\text{base}} + \max(\mathbf{u}_{\text{deterministic}}^\text{gust})|/u_\infty$, i.e., 2 in the current case. At the periphery this magnitude is nearly equal to $|\mathbf{u}_{\text{base}}|/u_\infty$, i.e., 1 in the present case. The fast central region tries to accelerate the slower areas in the surrounding, which decelerates the faster parts. Hence, the resulting (average) propagation velocity is in between $|\mathbf{u}_{\text{base}}|/u_\infty$ and $|\mathbf{u}_{\text{base}} + \max(\mathbf{u}_{\text{deterministic}}^\text{gust})|/u_\infty$.

During this propagation the shape of the wind gust changes. At first, the velocity distribution becomes asymmetric with respect to the
maximum velocity (see Fig. 7(d)). The fluid elements with the highest velocities in the gust structure flow faster than the slowest ones. A fast front forms and gradually detaches from the tail of the gust (see Fig. 7(f)) as observed in nature (Branlard, 2009).

One important feature which Fig. 7 can not display is clearly visible in Fig. 8. If the gust develops along the x-axis leading to streamwise velocities, which are larger than the base flow, in the region offside the axis the streamwise velocity drops to values below $u_∞$. That is obvious for example in Fig. 8(b), which depicts a trough in the streamwise velocity distribution next to the x-axis. Accordingly, the mass conservation is satisfied.

Similar to the streamwise velocity, the pressure distribution evolves in time: Starting from the sinusoidal form during the injection phase, it shows a kind of reversed Mexican hat shape (see Fig. 2(b)) when the gust travels through the domain.

The different physical observations concerning the propagation of the wind gust for the reference simulation are also valid for the simulations carried out in the following investigations.

### 6.2.2. Effect of the time scale $L_t$ on the effective length

During the generation of the gust three length scales and one time scale are required as shown in Section 2. However, because the Taylor hypothesis is assumed in ξ-direction, $L_ξ$ and $L_t$ are linked by relation (24). Therefore, the number of parameters reduces to three. $L_ξ$, $L_η$ and $L_ζ$ are parameters of the functions $f_ξ$, $f_η$ and $f_ζ$, respectively. Theoretically, for a given dimensionless length scale $L_ξ/L$, the length of the resulting injected gust has to be equal to $L_ξ/L$ and so on for the other two directions. In order to verify if the proposed source-term formulation delivers the correct effective lengths, simulations including a standard gust with an amplitude of $A_g/u_∞ = 1$ but with varying time scales $L_t u_∞/ L = \{0.25, 0.5, 1.0\}$ are carried out (see Table 2).

Fig. 9 compares the effective lengths of the injected gust in x- and y-direction with their associated theoretical values. The results are depicted along the x-axis with $y/L = z/L = 0$ for three different time scales. Furthermore, the distributions along the lateral axis are depicted. Since the results in z-direction are identical with those in y-direction, they are not shown here.

<figure><img src="image_12.png" /><figcaption>Fig. 9. Dimensionless streamwise velocity of a gust injected into an undisturbed flow for varying time scales $L_t u_∞/ L$ at different instants in time: On the left just after the injection ($t = L_t$), on the right at $t u_∞/ L = 2$. The results are taken along the x-axis with $y/L = z/L = 0$. Additionally, the maximum of each streamwise gust velocity profile is selected and its associated streamwise velocity distribution along the lateral y-axis is depicted. (Configuration of the gust: ECG-C² shape in ξ and t, standard ECG shape in η- and ζ-direction, $A_g/u_∞ = 1$, $L_η/L = L_ζ/L = 0.25$).</figcaption></figure>
The following observations can be made:

* Just after the injection of the gust at $t = L_t$ the shape of the gust in x-direction (except for the case $L_t u_\infty / L = 0.25$) is already no longer symmetric, but its effective length fits well to the associated theoretical value $L_\xi / L$. Obviously, the asymmetry increases with increasing values of $L_\xi / L$.

* After the injection the gust propagates through the flow field. A fast front and a slower tail form as mentioned before. The length of the tail is linked to $L_\xi / L$ or $L_t u_\infty / L$, respectively. A small value of $L_\xi / L$ leads to a short tail and vice versa.

* The effective length of the fast front is not strongly influenced by the size of $L_\xi / L$ or $L_t u_\infty / L$, respectively.

* The width of the gust shape in y-direction agrees very well with the theoretical value $L_\eta / L$ directly after the injection and also later. It is not affected by varying the time scale $L_t u_\infty / L$.

### 6.2.3. Effect of the length scale $L_\eta$ (or $L_\zeta$) on the effective length

To investigate the influence of the lateral length scale $L_\eta / L$ on the gust development, this input parameter is now varied within the range {0.25, 0.50, 1.0} (see Table 3). For this purpose, the time scale is fixed to the value $L_t u_\infty / L = 0.5$. Please note that due to the symmetry of the problem the same results can be obtained by varying $L_\zeta$ instead of $L_\eta$.

The results are depicted along the central x-axis ($y/L = z/L = 0$) in Fig. 10. Again, the distributions along the lateral axis are included in the graphs. The following observations can be made:

* Just after the injection of the gust at $t = L_t$ the effective length of the gust shape in x-direction does not depend on the variation of the length scale $L_\eta / L$.

* During the propagation phase the effective length of the gust shape in x-direction is affected by the variation of the initially prescribed length scale $L_\eta / L$. By increasing $L_\eta / L$, the separation between the fast front and the slower tail is attenuated, impacting the length of the gust.

* During the propagation phase the effective width of the gust in y-direction is no longer equal to its originally specified length scale $L_\eta / L$. Indeed, the part with the fast velocity, which is located in the middle of the gust, detaches and moves faster than the rest. The sideward (left and right) edges of the gust are slower. Therefore, the effective width of the gust in y-direction is reduced with increasing time.

### 6.2.4. Effect of the time scale $L_t$ on the strength of the gust

To investigate the influence of the time scale $L_t$ on the strength of the generated gust, the maximum dimensionless gust velocity obtained at each time step is time-averaged over the time period between the end of

the injection phase ($t = L_t$) and the given time $t u_\infty / L = 2$. The resulting value $\bar{u}_{\max}^{\text{sim}} / u_\infty$ is compared with the theoretical value $u_{\text{theory}} / u_\infty$. Table 2 provides an overview of the results obtained for different dimensionless amplitudes of the gust $A_g / u_\infty$ varied in the range {0.5, 1.0, 2.0}. Additionally, the dimensionless time scale $L_t u_\infty / L$ varies between the values {0.25, 0.50, 1.0}. For this analysis the lengths of the gust in $\eta$- and $\zeta$-direction are fixed to $L_\eta / L = L_\zeta / L = 0.25$. According to the results summarized in Table 2 the gust velocities obtained for the shortest time scale $L_t u_\infty / L = 0.25$ are too low (about -11% to -15% error) and the ones obtained for the largest time scale $L_t u_\infty / L = 1.0$ are too large (error between +4% and +15%). However, the error strongly decreases for $L_t u_\infty / L = 0.5$. Based on these results it can be concluded that a dimensionless time scale approximatively two times larger than the dimensionless length scales $L_\eta / L$ and $L_\zeta / L$ is required to achieve the best agreement with the intended maximum gust velocity.

### 6.2.5. Effect of the length scale $L_\eta$ (or $L_\zeta$) on the strength of the gust

For the investigation on the influence of the length scale $L_\eta / L$ in the range {0.25, 0.50, 1.0}, the time scale is fixed to $L_t u_\infty / L = 0.5$. As shown in Table 3 the error made on the time-averaged maximum gust velocity $\bar{u}_{\max}^{\text{sim}} / u_\infty$ as defined in the previous subsection is low and between +4% and -4% for the shortest length scales $L_\eta / L = L_\zeta / L = 0.25$. With an increasing length scale $L_\eta$, the relative error increases. Moreover, the increase of the gust amplitude $A_g / u_\infty$ also increases the relative error obtained.

### 6.2.6. Effect of the different shape functions

The final source-term formulation proposed in Eq. (14) relies on the four shape functions $f_t$, $f_\xi$, $f_\eta$ and $f_\zeta$ and also on the derivative of $f_\xi$ with respect to $\xi$. Since these shape functions determine the spatial and temporal structure of the gust, the different classic shape functions presented in Section 3 are evaluated and compared with the reference case.

In the present method the $\xi$-direction and the time $t$ are assumed to be linked by the Taylor hypothesis. Therefore, it is more consistent to use the same shape functions for $f_\xi$ and $f_t$ as done in the previous investigations. However, to distinguish the effects of each parameter in Eq. (14), the shape functions $f_\xi$ and $f_t$ will nevertheless be varied independently in this subsection.

First, solely the distribution in time $f_t$ is varied. The tested shapes are the ECG, EOG and the Gaussian function. The obtained results are compared in Fig. 11 with the reference simulation based on ECG-C$^2$. The other shapes $f_\xi$, $f_\eta$ and $f_\zeta$ are identical to the reference case ($f_\xi$=ECG-C$^2$, $f_\eta$=ECG and $f_\zeta$=ECG). The deviations between the velocity profiles in x- and y-direction obtained for ECG and ECG-C$^2$ are minimal as expected. Applying the Gaussian distribution the form of the gust in x-direction remains similar, but its strength is slightly decreased. This observation

<figure><img src="image_13.png" /><figcaption>Fig. 10. Dimensionless streamwise velocity of a gust injected into an undisturbed flow for varying length scales $L_\eta / L$ at different instants in time: On the left just after the injection ($t = L_t$), on the right at $t u_\infty / L = 2$. The results are taken along the x-axis with $y/L = z/L = 0$. Additionally, the maximum of each streamwise gust velocity profile is selected and its associated streamwise velocity distribution along the lateral y-axis is depicted. (Configuration of the gust: ECG-C$^2$ shape in $\xi$ and $t$, standard ECG shape in $\eta$- and $\zeta$-direction, $L_t u_\infty / L = 0.5$, $A_g / u_\infty = 1$, $L_\zeta / L = 0.25$).</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>**Fig. 11.** Dimensionless streamwise velocity of a gust injected into an undisturbed flow for varying temporal shape functions $f_t$ at different instants in time: On the left just after the injection ($t = L_t$), on the right at $t u_\infty/L = 2$. The results are taken along the x-axis with $y/L = z/L = 0$. Additionally, the maximum of each streamwise gust velocity profile is selected and its associated streamwise velocity distribution along the lateral y-axis is depicted. (Configuration of the gust: Standard ECG-C$^2$ shape in $\xi$, standard ECG shape in $\eta$- and $\zeta$-direction, $L_t u_\infty/L = 0.5$, $A_g/u_\infty = 1$, $L_\eta/L = L_\zeta/L = 0.25$).</figcaption></figure>

can be explained by the fact that the integral value of the Gaussian function over the injection phase is lower than those of ECG or ECG-C$^2$. The use of the Gaussian function for $f_t$ does not affect the gust structure in y-direction (or z-direction). However, the results predicted with EOG used for $f_t$ differ in form and strength. This can be attributed to the special form of the Mexican hat function (EOG), in which a local minimum occurs before the maximum is reached. Due to the normalization of the Mexican hat function by the difference between the minimum and the maximum value, the maximum itself is also 26% smaller than for the other functions. Since the source term quadratically depends on the gust

velocity, the effect on the resulting amplitude of the gust visible in Fig. 11 is even larger.

In the second series of predictions, $f_t$ is now fixed to ECG-C$^2$ as in the reference case and the distribution in $\xi$-direction $f_\xi$ is varied. The tested shapes are the ECG, ECG-C$^2$ and the Gaussian function. Note that the application of EOG in the main flow direction has been deliberately excluded. The reason for this decision is given by the fact that the special shape of the Mexican hat is lost very quickly during convective transport by the flow, since the fast fluid elements catch up with the slower part very quickly and thereby extinguish it. The shapes in the lateral

<figure><img src="image_15.png" /><figcaption>**Fig. 12.** Dimensionless streamwise velocity of a gust injected into an undisturbed flow for varying spatial shape functions $f_\xi$ at different instants in time: On the left just after the injection ($t = L_t$), on the right at $t u_\infty/L = 2$. The results are taken along the x-axis with $y/L = z/L = 0$. Additionally, the maximum of each streamwise gust velocity profile is selected and its associated streamwise velocity distribution along the lateral y-axis is depicted. (Configuration of the gust: Standard ECG-C$^2$ shape in time, standard ECG shape in $\eta$- and $\zeta$-direction, $L_t u_\infty/L = 0.5$, $A_g/u_\infty = 1$, $L_\eta/L = L_\zeta/L = 0.25$).</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Fig. 13.** Dimensionless streamwise velocity of a gust injected into an undisturbed flow for varying spatial shape functions $f_\eta$ at different instants in time: On the left just after the injection ($t = L_t$), on the right at $t u_\infty/L = 2$. The results are taken along the x-axis with $y/L = z/L = 0$. Additionally, the maximum of each streamwise gust velocity profile is selected and its associated streamwise velocity distribution along the lateral y-axis is depicted. (Configuration of the gust: Standard ECG-C$^2$ shape in $\xi$ and $t$, standard ECG shape in $\zeta$-direction, $L_t u_\infty/L = 0.5$, $A_g/u_\infty = 1$, $L_\eta/L = L_\zeta/L = 0.25$).</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>Fig. 14. Dimensionless streamwise velocity of a gust at different instants in time with different turbulence intensity values (TI = 0, 4 and 8%). Streamwise velocity along the x-axis at y/L = z/L = 0 (Configuration of the gust: ECG-C² shape in ξ and t, standard ECG shape in η- and ζ-direction, L_t u_∞/L = 0.5, A_g/u_∞ = 1, L_η/L = L_ζ/L = 0.25).</figcaption></figure>

directions $f_η$ and $f_ζ$ do not change ($f_η$=ECG and $f_ζ$=ECG). Fig. 12 compares the results with the reference case. By varying $f_ζ$ only marginal changes are visible in the velocity distributions along the central axis. Some minor deviations are observed in the strength of the gust. Furthermore, a change in $f_ζ$ does not perceptibly affect the structure of the gust in y- or z-direction.

Finally, the ECG-C² function is used for $f_t$ and $f_ζ$ and the shape function in η-direction is varied between ECG, EOG and Gaussian. Due to symmetry of the source-term formulation the same results would be obtained by changing $f_ζ$ instead of $f_η$. Therefore, only $f_η$ is investigated.

The gusts predicted are depicted in Fig. 13. Since the ECG and the Gaussian distributions are comparable (see Fig. 2(a) and (c)), the results obtained for these two shapes are very similar along the y-direction. Moreover, this change of the η-shape from ECG to Gaussian does not much affect the spatial distribution of the gust in x-direction and also in time. Interestingly, when $f_t$ changes from ECG to Gauss, it has a noticeable impact on the strength. However, changing $f_η$ from ECG to Gauss only leads to marginal deviations during the injection phase and small difference after propagation. Again EOG delivers different results due to its special shape: In y-direction the central peak is surrounded by two
<figure><img src="image_18.png" /><figcaption>**Fig. 15.** Three-dimensional view of the dimensionless streamwise velocity of a gust injected into turbulent flows (left: $Tl = 4\%$; right: $Tl = 8\%$) at different instants in time in the xy-plane at $z/L = 0$ (Reference configuration of the gust: ECG-$C^2$ shape in $\xi$ and $t$, standard ECG shape in $\eta$- and $\zeta$-direction, $A_g/u_\infty = 1$, $L_t u_\infty / L = 0.5$, $L_\eta / L = L_\zeta / L = 0.25$).</figcaption></figure>

areas with a velocity magnitude below the background flow velocity. These two areas are more pronounced than with ECG or Gauss. According to the theoretical shape, the width of the central peak is smaller than with ECG or Gauss and its strength is also lower. This is also visible in x-direction. It is worth to note that the resulting form of the velocity distribution in x-direction after the injection is apparently not affected by EOG, only its maximum. Since the strength of the EOG gust is lower after the injection phase, its propagation in streamwise direction is slower than for the ECG or Gaussian gusts. Additionally, the two lateral areas of lower velocity obtained with EOG propagate slower than the center of the gust and are located at the end of the gust tail. The effective width of the EOG fast front is comparable with those predicted by ECG or Gauss. Only its strength differs.

### 6.3. Propagation of wind gusts in turbulent flow

A major objective of the present work is to demonstrate the possibility to generate wind gusts within turbulent flows. In the previous sections the Reynolds number is high, but the flow is undisturbed. In order to obtain turbulent flows synthetic turbulent perturbations with different

turbulence intensities are added to the background flow of the following simulations as explained in Section 5. In these simulations the configuration of the wind gust is identical to the reference case in Section 6.2.1.

Before discussing the results the following should be noted about this test case: The presently used simplifications of generating homogeneous isotropic inflow turbulence (see Section 5) is solely an idealized case for the present study and does not claim to be a description of atmospheric turbulence. Of course, the flow in an atmospheric boundary layer requires inhomogeneous and anisotropic turbulent flow fields. However, two issues have to be taken into account here: First, in the present setup a free flow without the effect of any wall and thus boundary layer is considered with the single objective to demonstrate that the formulation of the source term for the gust and the background turbulence can be superimposed. Second, either with the presently used inflow generator or any other inflow generator (artificial or based on precursor simulations) the present restrictions to homogeneous isotropic turbulence can be easily overcome, which is the topic of ongoing work.

Fig. 14 depicts the evolution of the injected gust in time and space for the case with $Tl = 4\%$ and $8\%$ along the x-axis at $y/L = z/L = 0$. To allow a comparison with the undisturbed case the streamwise velocity
and the pressure are also plotted for the previous case of $TI = 0\%$. Similar to the undisturbed flow depicted in Fig. 8 (left) a 3D view of the streamwise velocity is added in Fig. 15.

First, during the injection phase the wind gust shape is only marginally affected by the turbulence perturbations. This observation can be explained by the fact that the amplitude of the perturbations are small compared to the amplitude of the gust. Second, the propagation velocity of the gust is not influenced by the turbulent fluctuations, as visible in Fig. 14(c) and (d). However, the strength of the gust after the injection phase is reduced by the increase of the $TI$ value. The turbulence present in the background augments the destruction of the gust structure. For example, the maximal streamwise velocity of the gust obtained with $TI = 8\%$ at $t u_\infty/L = 2.0$ is only $1.4 u_\infty$ and the one predicted within an undisturbed flow reached $1.93 u_\infty$ (see Fig. 14(d)).

Please note that, if the strength decreases with the increase of $TI$, the effective length in x-direction of the generated gust is not changing significantly. The effective length of the fast front and the tail remains similar. Moreover, no dispersion of the gust in y- and z-direction is noticeable due to the superimposed turbulent fluctuations (see Fig. 15(d)).

## 7. Conclusions

Based on the long-term goal of describing the effect of turbulent wind gusts on flexible structures in the context of coupled fluid-structure interaction simulations, an efficient method for introducing the gusts into the computational domain was developed in a first step. The objectives were to avoid CPU-time intensive techniques and to take the feedback effect of the flow field on the gust itself into account. Thus, a source-term formulation for the momentum conservation equation was derived, which allows to place the injection region of the gust close to the region of interest. Since no source term is added to the continuity equation, the local and global mass conservation is not touched and the divergence-free condition is still satisfied for incompressible fluids. Besides applying the classical shape functions for gusts (ECG, EOG) to define the source term, a modified $C^2$-"1-cosine" shape ECG-$C^2$ was derived which guarantees that the source-term distribution has no kink. After a grid independence study a thorough analysis of the new gust injection technique was carried out in undisturbed and turbulent background flows. The main findings are:

* At the end of the injection phase the predicted shape of the wind gust is very similar to the intended one. Consequently, the proposed source-term formulation works as desired.

* When the source term is no longer active after the injection phase, the conservation equations of mass and momentum are fully satisfied. Consequently, the generated flow field including the wind gust are a valid solution of the Navier-Stokes equations.

* The injected gust travels through the flow field with a nearly constant propagation velocity. In this phase, a fast and steeper front as well as a slower tail of the gust develop as observed in nature.

* The effective length of the fast front does not strongly depend on the size $L_\zeta/L$ (or $L_t u_\infty/L$, respectively), whereas the length of the slower tail is linked to this quantity.

* The prescribed length scale of the gust in the lateral directions significantly influences the separation between the fast front and the slower tail.

* An assessment of the different shape functions yields the finding that the most significant changes in the form and strength of the generated gust is achieved by the application of the EOG shape function, both in time and in $\eta$ or $\zeta$-direction.

* A parameter study on the effect of the time scale $L_t$ has shown that the dimensionless time scale should be about two times larger than the dimensionless length scales $L_\eta/L$ and $L_\zeta/L$.

* The proposed source-term formulation for the wind gusts works in laminar and turbulent flows. It is not restricted to eddy-resolving simulation approaches (DNS or LES) and can also be combined with Reynolds-averaged Navier-Stokes solvers or hybrid LES-URANS techniques.

* The source-term formulation can be joined with other source-term approaches such as the turbulence generator used in the present study.

* The background turbulence does not perceptibly affect the propagation velocity of the gust. However, the strength of the gust is significantly attenuated by the turbulent fluctuations and this trend increases with raising turbulence intensities of the underlying flow field.

Next steps are the application of the gust injection technique in practically relevant flow cases (e.g., in atmospheric boundary layers with inhomogeneous and anisotropic turbulent flow fields) and later on the extension towards FSI.

## Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

## Acknowledgments

The project was initiated with funding from the Interne Forschungsförderung (IFF) of the Helmut-Schmidt-Universität Hamburg. In the meantime, it is financially supported by the Deutsche Forschungsgemeinschaft under the contract number BR 1847/17-1.

## Appendix A. Recipe to apply source-term formulation

The source-term formulation for prescribing gusts within the computational domain is neither restricted to specific (deterministic) gust shapes nor to specific CFD codes. In the following the general procedure how to use the formulation for different gust shapes and within different CFD packages is briefly described:

### 1. Choose a gust shape as a function of space and time:

Besides the typically used deterministic functions described in Section 3 (ECG, EOG, Gaussian) the description of any shape in space and time is possible. Nevertheless, since the current source-term formulation includes a derivative of the function and possesses non-zero values only for $\xi < \xi_g$ (see Eq. (14)) it is advantageous but not mandatory that the second derivative is continuous. In the present study that leads to ECG-$C^2$ (see Eq. (23)).

### 2. Subroutine or user function for the source term:

Depending on the underlying CFD code, the source term for the gust generation according to Eq. (14) can be formulated either in a extra subroutine
or in a user function as available in many commercial codes. Note that in contrast to other methods such as the field velocity method or the split velocity method no further modifications of the convective terms or the time integration are required for the present source-term formulation. Thus, its implementation is straightforward and can be done even in commercial codes.

Furthermore, the source-term formulation described above is written down here for the integral form of the Navier-Stokes equations, which is the basis for the finite-volume scheme applied in this study. If the solution method for the governing equations is based on the differential form of the equations (i.e., for a finite-difference scheme) the formulation in Eq. (14) has to be adjusted by canceling the volume of the control volume. Apart from that the source-term formulation is still valid and directly applicable.

### 3. Consideration of turbulence:

If intended, additional turbulence can either be added at the inflow of the computational domain or as an additional source-term formulation inside the domain (Breuer, 2018; De Nayer et al., 2018b; Schmidt and Breuer, 2017). Again the user is free to either artificially generate the turbulent contributions (e.g., by the digital filter concept, synthetic eddy method, proper orthogonal decompositions, Fourier techniques, ..., see review by Tabor and Baba-Ahmadi (2010)) or to reuse data of auxiliary precursor simulations.

## References

Bierbooms, W., 2005. Investigation of spatial gusts with extreme rise time on the extreme loads of pitch-regulated wind turbines. Wind Energy: An Int. Journal for Progress and Applications in Wind Power Conversion Technology 8 (1), 17–34.

Bierbooms, W., 2006. Constrained stochastic simulation – generation of time series around some specific event in a normal process. Extremes 8 (3), 207–224.

Bierbooms, W., 2007. Specific gust shapes leading to extreme response of pitch-regulated wind turbines. In: Journal of Physics: Conference Series, vol. 75. IOP Publishing, 012058.

Bierbooms, W., Cheng, P.-W., 2002. Stochastic gust model for design calculations of wind turbines. J. Wind Eng. Ind. Aerod. 90 (11), 1237–1251.

Bierbooms, W., Dragt, J.B., Clejne, H., 1999. Verification of the mean shape of extreme gusts. Wind Energy 2, 137–150.

Bos, R., 2017. Extreme Gusts and Their Role in Wind Turbine Design. Ph.D. thesis. Delft University of Technology, The Netherlands.

Branlard, E., 2009. Wind Energy: on the Statistics of Gusts and Their Propagation through a Wind Farm. Master's thesis. Energy Research Center of The Netherlands, Petten, The Netherlands, ECN-Wind-Memo-09-005.

Brasseur, O., 2001. Development and application of a physical approach to estimating wind gusts. Mon. Weather Rev. 129 (1), 5–25.

Breuer, M., 2018. Effect of inflow turbulence on an airfoil flow with laminar separation bubble: an LES study. Flow, Turbul. Combust. 101 (2), 433–456.

Breuer, M., De Nayer, G., Münsch, M., Gallinger, T., Wüchner, R., 2012. Fluid-structure interaction using a partitioned semi-implicit predictor-corrector coupling scheme for the application of large-eddy simulation. J. Fluid Struct. 29, 107–130.

Burton, T., Sharpe, D., Jenkins, N., Bossanyi, E., 2001. Wind Energy Handbook. John Wiley & Sons.

De Nayer, G., Apostolatos, A., Wood, J.N., Bletzinger, K.-U., Wüchner, R., Breuer, M., 2018a. Numerical studies on the instantaneous fluid-structure interaction of an air-inflated flexible membrane in turbulent flow. J. Fluid Struct. 82, 577–609.

De Nayer, G., Breuer, M., 2014. Numerical FSI investigation based on LES: flow past a cylinder with a flexible splitter plate involving large deformations (FSI-PfS-2a). Int. J. Heat Fluid Flow 50, 300–315.

De Nayer, G., Breuer, M., Perali, P., Grollmann, K., 2019. Modeling of wind gusts for large-eddy simulations related to fluid-structure interactions. In: Salvetti, M., Armenio, V., Fröhlich, J., Geurts, B.J., Kuerten, H. (Eds.), ERCOFTAC Series, Direct and Large-Eddy Simulation XI, 11th Int. ERCOFTAC Workshop on Direct and Large-Eddy Simulation: DLES-11, May 29–31, 2017, vol. 25. Springer Nature Switzerland AG 2019, Pisa, Italy, pp. 453–459.

De Nayer, G., Breuer, M., Wood, J.N., 2020. Numerical investigations on the dynamic behavior of a 2-DOF airfoil in the transitional Re number regime based on fully coupled simulations relying on an eddy-resolving technique. Int. J. Heat Fluid Flow 85, 108631.

De Nayer, G., Kalmbach, A., Breuer, M., Sicklinger, S., Wüchner, R., 2014. Flow past a cylinder with a flexible splitter plate: a complementary experimental-numerical investigation and a new FSI test case (FSI-PfS-1a). Comput. Fluid 99, 18–43.

De Nayer, G., Schmidt, S., Wood, J.N., Breuer, M., 2018b. Enhanced injection method for synthetically generated turbulence within the flow domain of eddy-resolving simulations. Comput. Math. Appl. 75 (7), 2338–2355.

Durst, F., Schäfer, M., 1996. A parallel block-structured multigrid method for the prediction of incompressible flows. Int. J. Numer. Methods Fluid. 22 (6), 549–565.

Emeis, S., 2018. Wind Energy Meteorology: Atmospheric Physics for Wind Power Generation. Springer.

Ferziger, J.H., Perić, M., 2002. Computational Methods for Fluid Dynamics, third ed. Springer Berlin.

Feser, F., Barcikowska, M., Krueger, O., Schenk, F., Weisse, R., Xia, L., 2015. Storminess over the north atlantic and northwestern europe – a review. Q. J. R. Meteorol. Soc. 141 (687), 350–382.

Förster, M., Breitsamter, C., 2015. Aeroelastic prediction of discrete gust loads using nonlinear and time-linearized CFD-methods. Journal of Aeroelasticity and Structural Dynamics 3 (3), 19–38.

Frost, W., Long, B.H., Turner, R.E., 1978. Engineering Handbook on the Atmospheric Environmental Guidelines for Use in Wind Turbine Generator Development. Tech.

rep., Tennessee Univ., Tullahoma (USA). Space Inst.; National Aeronautics and Space Administration. George C. Marshall Space Flight Center, Huntsville, AL (USA).

Germano, M., Piomelli, U., Moin, P., Cabot, W.H., 1991. A dynamic subgrid-scale eddy viscosity model. Phys. Fluid. 3, 1760–1765.

Ghoreyshi, M., Greisz, I., Jirasek, A., Satchell, M., 2018a. Simulation and modeling of rigid aircraft aerodynamic responses to arbitrary gust distributions. Aerospace 5 (2), 43.

Ghoreyshi, M., Jirasek, A., Greisz, I., Satchell, M., 2018b. Indicial aerodynamic modeling for arbitrary gust responses. In: 2018 Applied Aerodynamics Conference. AIAA-2018-3965.

Gromke, C., Ruck, B., 2018. On wind forces in the forest-edge region during extreme-gust passages and their implications for damage patterns. Boundary-Layer Meteorol. 168 (2), 269–288.

Gross, J.L., Main, J.A., Phan, L.T., Sadek, F.H., Cauffman, S.A., Jorgensen, D.P., 2013. Final report on the collapse of the Dallas Cowboys indoor practice facility. (Accessed 2 May 2009).

Hawbecker, P., Basu, S., Manuel, L., 2017. Realistic simulations of the july 1, 2011 severe wind event over the buffalo ridge wind farm. Wind Energy 20 (11), 1803–1822.

Heinrich, R., 2014. Simulation of interaction of aircraft and gust using the TAU-code. In: Dillmann, A., Heller, G., Krämer, E., Kreplin, H.P., Nitsche, W., Rist, U. (Eds.), New Results in Numerical and Experimental Fluid Mechanics IX. Notes on Numerical Fluid Mechanics and Multidisciplinary Design, vol. 124. Springer, pp. 503–511.

IEC-Standard, 2002. 61400-21. Measurement and Assessment of Power Quality of Grid Connected Wind Turbines.

Kamimura, K., Gardiner, B., Dupont, S., Finnigan, J., 2019. Agent-based modelling of wind damage processes and patterns in forests. Agric. For. Meteorol. 268, 279–288.

Kasperski, M., 2007. A consistent model for the codification of wind loads. J. Wind Eng. Ind. Aerod. 95 (9), 1114–1124.

Kelleners, P., Heinrich, R., 2016. Simulation of interaction of aircraft with gust and resolved LES-simulated atmospheric turbulence. In: Radespiel, R., Niehuis, R., Kroll, N., Behrends, K. (Eds.), Advances in Simulation of Wing and Nacelle Stall. FOR-1066, 2014. Notes on Numerical Fluid Mechanics and Multidisciplinary Design, vol. 131. Springer, pp. 203–221.

Khosla, P.K., Rubin, S.G., 1974. A diagonally dominant second-order accurate implicit scheme. Comput. Fluid 2 (2), 207–209.

Klein, M., Sadiki, A., Janicka, J., 2003. A digital filter based generation of inflow data for spatially-developing direct numerical or large-eddy simulations. J. Comput. Phys. 186, 652–665.

Knigge, C., Raasch, S., 2016. Improvement and development of one- and two-dimensional discrete gust models using a large-eddy simulation model. J. Wind Eng. Ind. Aerod. 153, 46–59.

Larsen, G.C., Bierbooms, W., Hansen, K.S., 2003. Mean Gust Shapes. Tech. rep., Riso National Laboratory, Riso-R. No. 1133(EN).

Menegozzo, L., Dal Monte, A., Benini, E., Benato, A., 2018. Small wind turbines: a numerical study for aerodynamic performance assessment under gust conditions. Renew. Energy 121, 123–132.

National Institute of Standards and Technology, 2020. Disaster & failure studies. https://www.nist.gov/topics/disaster-failure-studies.

Nicoud, F., Ducros, F., 1999. Subgrid-scale stress modelling based on the square of the velocity gradient tensor. Flow, Turbul. Combust. 62 (3), 183–200.

Nielsen, M., Larsen, G.C., Mann, J., Ott, S., Hansen, K.S., Pedersen, B.J., 2003. Wind Simulation for Extreme and Fatigue Loads. Tech. rep., Riso National Laboratory, Riso-R-1437 (EN).

Norris, S.E., Cater, J.E., Stol, K.A., Unsworth, C.P., 2010. Wind turbine wake modelling using large-eddy simulation. Proceedings of the 17th Australian Fluid Mechanics Conference. Curran Associates, Inc., University of Auckland, Australia.

Onol, A.O., Yesilyurt, S., 2017. Effects of wind gusts on a vertical axis wind turbine with high solidity. J. Wind Eng. Ind. Aerod. 162, 1–11.

Parameswaran, V., Baeder, J.D., 1997. Indicial aerodynamics in compressible flow – direct computational fluid dynamic calculations. J. Aircraft 34 (1), 131–133.

Pasquetti, R., Peres, N., 2015. A penalty model of synthetic micro-jet actuator with application to the control of wake flows. Comput. Fluid 114, 203–217.

Radespiel, R., Francois, D.G., Hoppmann, D., Klein, S., Scholz, P., Wawrzinek, K., Lutz, T., Auerswald, T., Bange, J., Knigge, C., Raasch, S., Kelleners, P., Heinrich, R., Reuß, S.,
Probst, A., Knopp, T., 2013. Simulation of wing stall. In: 43rd AIAA Fluid Dynamics Conference. AIAA-2013-3175.

Rhie, C.M., Chow, W.L., 1983. Numerical study of the turbulent flow past an airfoil with trailing-edge separation. AIAA J. 21 (11), 1525-1532.

Schmidt, S., Breuer, M., 2017. Source term based synthetic turbulence inflow generator for eddy-resolving predictions of an airfoil flow including a laminar separation bubble. Comput. Fluid 146, 1-22.

Singh, R., Baeder, J.D., 1997. Direct calculation of three-dimensional indical lift response using computational fluid dynamics. J. Aircraft 34 (4), 465-471.

Smagorinsky, J., 1963. General circulation experiments with the primitive equations I: the basic experiment. Mon. Weather Rev. 91 (3), 99-165.

Steyn, D.G., Oke, T.R., Hay, J.E., Knox, J.L., 1981. On scales in meteorology and climatology. Climatological Bulletin 39, 1-8.

Stone, H.L., 1968. Iterative solution of implicit approximations of multidimensional partial differential equations. SIAM J. Numer. Anal. 5 (3), 530-558.

Storey, R.C., Norris, S.E., Cater, J.E., 2013. Large eddy simulation of wind events propagating through an array of wind turbines. In: Proceedings of the World Congress on Engineering, vol. 3. London, UK.

Storey, R.C., Norris, S.E., Cater, J.E., June, 2014. Modelling turbine loads during an extreme coherent gust using large eddy simulation. J. Phys. Conf. 524, 012177.

Tabor, G., Baba-Ahmadi, M., 2010. Inlet conditions for large-eddy simulation: a review. Comput. Fluid 39, 553-567.

Tischmacher, M., Ruck, B., 2013. Interaction of gusts and forest edges - an experimental wind-tunnel study. Forestry: Int. J. Financ. Res. 86 (5), 523-532.

Volpe, R., Da Silva, A., Ferrand, V., Le Moyne, L., 2013. Experimental and numerical validation of a wind gust facility. J. Fluid Eng. 135 (1), 011106.

Wales, C., Jones, D., Gaitonde, A., 2014. Prescribed velocity method for simulation of aerofoil gust responses. J. Aircraft 52 (1), 64-76.

Weishäupl, C., Laschka, B., 2001. Euler solutions for airfoils in inhomogeneous atmospheric flows. J. Aircraft 38 (2), 257-265.

Wood, J.N., De Nayer, G., Schmidt, S., Breuer, M., 2016. Experimental investigation and large-eddy simulation of the turbulent flow past a smooth and rigid hemisphere. Flow, Turbul. Combust. 97 (1), 79-119.

Yawar, A., Ebrahem, M., Manzoor, S., Sheikh, N.A., Ali, M., 2019. Transient cross flow and heat transfer over a rotationally oscillating cylinder subjected to gust impulse. Int. J. Heat Mass Tran. 137, 108-123.

Yokel, F.Y., Yancey, C.W., Cattaneo, L.E., Marshall, R.D., 1976. Investigation of wind damage in the metropolitan Washington, DC area, april 3-4, 1975 (NBS TN 909). Tech. rep.