# Assessment of two wind gust injection methods: Field velocity vs. split velocity method

Khaled Boulbrachene, Guillaume de Nayer, Michael Breuer

► To cite this version:

Khaled Boulbrachene, Guillaume de Nayer, Michael Breuer. Assessment of two wind gust injection methods: Field velocity vs. split velocity method. Journal of Wind Engineering and Industrial Aerodynamics, 2021, 218, pp.104790. 10.1016/j.jweia.2021.104790. hal-03374308

HAL Id: hal-03374308

https://hal.science/hal-03374308

Submitted on 13 Oct 2021

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>

# Assessment of two wind gust injection methods: Field velocity vs. split velocity method

Khaled Boulbrachene, Guillaume De Nayer, Michael 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:**

Field velocity method

Split velocity method

Vertical and horizontal wind gust

Feedback effect

ECG

## ABSTRACT

The objective of the present paper is to revisit two well-known wind gust injection methods in a consistent manner and to assess their performance based on different application cases. These are the field velocity method (FVM) and the split velocity method (SVM). For this purpose, both methods are consistently derived pointing out the link to the Arbitrary Lagrangian Eulerian formulation and the geometric conservation law. Furthermore, the differences between FVM and SVM are worked out and the advantages and disadvantages are compared. Based on a well-known test case considering a vertical gust hitting a plate and a newly developed case taking additionally a horizontal gust into account, the methods are evaluated and the deviations resulting from the disregard of the feedback effect in FVM are assessed. The results show that the deviations between the predictions by FVM and SVM are more pronounced for the horizontal gust justifying the introduction of this new test case. The main reason is that the additional source term in SVM responsible for the feedback effect of the surrounding flow on the gust itself nearly vanishes for the vertical gust, whereas it has a significant impact on the flow field and the resulting drag and lift coefficients for the horizontal gust. Furthermore, the correct formulation of the viscous stress tensor relying on the total velocity as done in case of SVM plays an important role, but is found to be negligible for the chosen Reynolds number of the present test cases. The study reveals that SVM is capable of delivering physical results in contradiction to FVM. It paves the way for investigating further complex gust configurations (e.g., inclined gusts) and practical applications towards coupled fluid-structure interaction simulations of engineering structures impacted by wind gusts.

## 1. Introduction

The study of extreme events constitutes an important aspect of many engineering applications nowadays. Investigating extreme events was found to enable the design of safe components that will withstand harsh conditions and circumvent any sudden failure which may occur. In the field of fluid mechanics, wind gusts represent a typical example of such extreme phenomena which might lead to devastating consequences on flexible as well as rigid structures. For instance, aircraft design requires the study of different load cases and the evaluation of their aerodynamic responses (Wu et al., 2019). This ensures that the encounter of aircraft by wake vortices of a preceding plane (Sitaraman and Baeder, 2006; Gordnier and Visbal, 2015), a streamwise wind gust (Granlund et al., 2014) or an impulsive change in the angle of attack (Parameswaran and Baeder, 1997; Singh and Baeder, 1997) is controlled properly. In order to guarantee the safety of trains moving over bridges and to evaluate the risk of derailment, Montenegro et al. (2020) simulated the effect of discrete gusts and turbulent wind models on aerodynamic forces. That allows to determine critical wind speeds,

which is especially of interest for high-speed trains. Beside aeronautics and rail transport, the wind energy sector is constantly dealing with the design of structurally reliable wind turbine tower systems (Kwon et al., 2012) and especially wind turbine blades (for example, see Onol and Yesilyurt (2017)). A turbine blade is, like a helicopter blade, exposed to trailing-edge vortices in the wake of other blades as well as to potential gusts which may strip it away. For example, Hawbecker et al. (2017) reported on extreme structural damages to turbines of the Buffalo Ridge Wind Farm in Minnesota in 2011. Blades from multiple turbines broke away and a tower buckled in the intense winds. They attempted to characterize meteorological conditions over the wind farm area during this event and carried out weather research and forecasting model simulations of the event that considered current and anticipated future operational model setups.

More recently, lightweight thin structures such as membrane-covered and tensile-membrane structures have become widely used in civil engineering applications. Due to the unique flexible characteristics of the fabric membranes such structures allow the engineers

<footnote>* Corresponding author.
E-mail address: breuer@hsu-hh.de (M. Breuer).</footnote>
to create visually exciting and iconic structures while maintaining a faster installation and an overall lower cost in comparison to traditional constructions. However, these structures are vulnerable to wind events and to determine to what extent the aerodynamic loads influence their integrity is essential (for example, see Gross et al. (2013)). Although design standards are prescribed to provide guidelines in order to take such events into account (for example, see Frost et al. (1978), IEC-Standard (2002), Burton et al. (2001), Kasperski (2007)), these are based on simplifying assumptions to facilitate the design process. Therefore, the development of methods capable of providing a rapid and comprehensive investigation of proposed designs would substantially contribute to the design procedure. Highly resolved numerical simulations of the flow field (and the resulting deformations) are considered to be a reliable and powerful method for this objective. That was recently shown in De Nayer et al. (2014), De Nayer and Breuer (2014), De Nayer et al. (2018) and Apostolatos et al. (2019), so far not for gusts but at least for turbulent flows around flexible membranous structures.

In the following, the description is restricted to wind gusts and methodologies to take these short but often violent phenomenon into account. One of the conventional approaches to inject wind gusts into the computational domain is to superimpose them at the inlet boundary and to allow the gusts to freely propagate through the domain. This technique is known in the literature as the far-field boundary condition (FBC) method. To avoid the violation of the global mass balance especially for incompressible fluid flows, this method requires the dynamic adaptation of the inlet velocity away from the gust to ensure a constant total mass flow through the inlet (Norris et al., 2010; De Nayer et al., 2019). Such a requirement is necessary for the convergence of incompressible flow solvers. Owing to the fact that inflow boundaries are often far away from the region of interest, a relatively coarse grid resolution is used in the upstream region to reduce the computational costs. Consequently, numerical dissipation is prevalent in these regions and the induced disturbances are strongly damped while being convected through the flow domain. To circumvent this problem, one would have to use a highly resolved grid throughout the entire computational domain, or a secondary fine moving mesh transporting the gust within the computational domain (e.g., an overset-grid technique such as the resolved gust approach of Heinrich (2014)). Nevertheless, the aforementioned remedies are computationally highly demanding leading to long simulation times.

Several approaches have been proposed to overcome the requirement of using fine grids throughout the computational domain. A category of these methods is referred to as the pure source-term formulation method. Recently, De Nayer and Breuer (2020) proposed an original method for injecting wind gusts in the computational domain based on a source-term formulation. The authors derived a source-term formulation for the momentum conservation equations which allows to place the injection region of the gust close to the region of interest. The gust can then freely propagate through the domain and, consequently, the full coupling (i.e., interaction) between the gust and the surrounding flow is captured.

Contrary to the aforementioned methods which allow gusts to freely travel after the injection phase, another family of methods which defines the position of the imposed field at each instant in time is denoted the prescribed velocity methods. There exist two prescribed velocity methods, namely, the field velocity method (FVM) and the split velocity method (SVM). The field velocity or grid velocity approach (also denoted the disturbance velocity approach by Heinrich and Reimer (2013)) was initially proposed by Parameswaran and Baeder (1997) and Singh and Baeder (1997). FVM works by incorporating gust velocities by suitably modifying the grid velocity and consequently modifying the grid time metrics. However, this modification does not lead to a real grid deformation. Hence, a pseudo grid motion is simulated without an actual grid variation. FVM could be thought of as an extension of the surface transpiration method (Sankar et al., 1986), where instead of applying the velocity correction only to the surface, it is applied to the

entire flow field. Since FVM is similar to moving grid problems, grid time metrics have to be consistently evaluated to satisfy the geometric conservation law (GCL) (Sitaraman and Baeder, 2006).

The history of FVM started in the field of aeronautics, where its main objective was to tackle the problem of evaluating indicial responses (i.e., aerodynamic loads) due to a step change in one of the influencing parameters (e.g., angle of attack or pitch rate) in the context of inviscid three-dimensional compressible flows while having full insight into the flow features (velocity and pressure distributions) (Parameswaran and Baeder, 1997; Singh and Baeder, 1997). The method could successfully solve the challenge of decoupling the step change in the angle of attack from the pitch rate, delivered results that agree with the attainable analytical results, and avoided numerical instabilities associated with this problem. Moreover, in the context of inviscid Euler equations FVM was employed to simulate the response of an airfoil penetrating through a sharp-edged vertical gust (Parameswaran and Baeder, 1997; Sitaraman and Baeder, 2000), and to model the interaction of an isolated vortex with a rotor blade (Singh and Baeder, 1996). Furthermore, FVM was used to aid the prediction of aerodynamic loads on a helicopter blade by incorporating the velocity field caused by trailed tip vortices from all other blades (Sitaraman, 2003; Sitaraman and Baeder, 2006). Nevertheless, FVM is incapable of accounting for the effect of the surrounding fluid flow on the induced field.

In contrast, SVM (originally proposed by Wales et al. (2014)) fully captures the interaction between the gust and the object of interest. SVM is based on decomposing the velocity field to a prescribed gust velocity and the remaining velocity (or background velocity as denoted by Huntley et al. (2017)). The momentum equations are effectively rearranged to be solved for the background velocity. As a result, additional source terms are derived to account for the feedback effect which is missing in FVM. Initial applications of SVM concentrated on applying 1-cosine vertical gusts on several NACA airfoils with different gust lengths and amplitudes (Wales et al., 2014). The simulations, however, were performed in the context of two-dimensional Euler equations. Later, Huntley et al. (2017) extended SVM to investigate the gust response of a flexible aircraft to a series of 1-cosine gusts in the context of three-dimensional viscous flows. Both research groups compared their results with those obtained by FVM and found SVM to be more accurate when short wavelength gusts are imposed. Recently, Biler et al. (2019), Badrya and Baeder (2019) and Badrya et al. (2021) investigated the flow physics resulting from a flat plate-gust encounter. The analysis showed the ability of analytical models (Küssner's solution (Küssner, 1936)) in predicting gust response for a 0° geometric angle of attack while the need for CFD solutions remains in order to formulate empirical formulas predicting responses for high angles of attack.

The steps in this paper are as follows:

1. Provide a comprehensive description of the FVM and SVM methods along with the modifications needed to be made in the governing equations.

2. Deliver best practice guidelines describing the implementation of the methods in a finite-volume code thoroughly.

3. Apply the methods to simulate the encounter of a vertical as well as a horizontal gust (which can be said to be relevant for many engineering applications) on a horizontal flat plate.

4. Highlight the differences among the methods and their limitations.

The main objective is to figure out the most appropriate methodology to simulate the effect of wind gusts impacting on mechanical or civil engineering structures. For the latter especially flexible membranous structures impacted by horizontal gusts are in the focus of medium-term interest, whereby such kind of complex fluid-structure interaction problems are to be simulated by high-fidelity eddy-resolving schemes such as the large-eddy simulation (LES) technique.
## 2. Prescribed gust methods

In the following, two well-established methods for injecting gusts into a CFD simulation by prescribing the position in the computational domain, the shape and the velocity of the gusts at each time step are introduced. In the present work the gust velocity is denoted $u_{g,j}$ ($i \in \{1, 2, 3\}$) in Cartesian coordinates. The flow field surrounding the gust is denoted as the background velocity field $\tilde{u}_i$. The Navier-Stokes equations for an incompressible fluid are solved for this background velocity. The prescribed velocity methods postulate that the divergence of the background velocity is equal to zero (i.e., $\frac{\partial \tilde{u}_i}{\partial x_j} = 0$). To present the total flow field $u_i$ containing the resolved flow field and the prescribed gust, the sum of both contributions (i.e., $u_i = \tilde{u}_i + u_{g,j}$) has to be considered by superposition. Since the divergence of the prescribed gust velocity on a fixed grid (i.e., $\frac{\partial u_{g,j}}{\partial x_i}$) is not necessarily equal to zero, the divergence of the total velocity will not automatically be zero. In order to introduce the prescribed gust velocity into the Navier-Stokes equations while maintaining the postulate of $\frac{\partial \tilde{u}_i}{\partial x_i} = 0$, the Arbitrary Lagrangian Eulerian (ALE) formulation of the Navier-Stokes equations is required as described below.

### 2.1. Field Velocity Method (FVM)

The original idea of this method is the incorporation of unsteady flow conditions via a pseudo-grid movement in the computational domain according to the Arbitrary Lagrangian Eulerian (ALE) formulation. In this way, disturbance velocities (e.g., gust velocities) could be prescribed by changing the grid time metrics while the grid is actually not altered or distorted.

In the context of a temporally varying domain (as in fluid-structure interaction (FSI), see Breuer et al. (2012) for example) the Navier-Stokes equations (mass and momentum) are written in the Arbitrary Lagrangian Eulerian (ALE) form. For an incompressible fluid, the integral form of these equations read (here in a Cartesian coordinate system for the sake of simplicity):

$$ \frac{d}{dt} \int_{V(t)} \rho dV + \int_{S(t)} \rho (\tilde{u}_j - u_j^{\text{grid}}) \cdot n_j \, dS = 0, \quad (1) $$

$$ \frac{d}{dt} \int_{V(t)} \rho \tilde{u}_i \, dV + \int_{S(t)} \rho \tilde{u}_i (\tilde{u}_j - u_j^{\text{grid}}) \cdot n_j \, dS \\ = - \int_{S(t)} \tilde{\tau}_{ij} \cdot n_j \, dS - \int_{S(t)} p \cdot n_i \, dS. \quad (2) $$

Here, the control volumes (CV) have time-dependent volumes $V(t)$ and surfaces $S(t)$. Since in this case the grid is deformable, the grid velocity with which the surface of a CV is moving, is taken into account via $u_j^{\text{grid}}$. $\tilde{u}_i$ describes the velocity of the flow field. A grid deformation in the inner domain will not affect $\tilde{u}_i$. However, the flow velocity $\tilde{u}_i$ changes due to movements or deformations of walls. $\tilde{\tau}_{ij}$ is the stress tensor based on the flow velocity $\tilde{u}_i$.

In the context of a temporally varying domain the ALE formulation has to additionally satisfy the so-called space conservation law (SCL) (Demirdžić and Perić, 1988, 1990) or geometric conservation law (GCL) (Lesoinne and Farhat, 1996):

$$ \frac{d}{dt} \int_{V(t)} dV - \int_{S(t)} u_j^{\text{grid}} \cdot n_j \, dS = 0. \quad (3) $$

This extra conservation law is necessary to assure that no space is lost during the change of CVs. Employing the GCL, the contribution of the grid movement to the mass fluxes will cancel the unsteady term leading to the reduction of the mass conservation equation to:

$$ \int_{S(t)} \rho \tilde{u}_j \cdot n_j \, dS = 0, \quad (4) $$

which is equivalent to the mass conservation equation on a fixed grid. Consequently, the pressure-correction equation in the context of

projection methods does not have to be modified. The discrete GCL is used to evaluate the additional grid fluxes in the momentum equation by evaluating the change of the position or the shape of a CV to compute the unknown grid velocity $u_i^{\text{grid}}$.

The field velocity approach is built per analogy based on the previously written ALE formulation. The grid is now undeformed and thus fixed, i.e., the cell volumes $V$ and surfaces $S$ are no longer time-dependent. Furthermore, the grid velocity $u_i^{\text{grid}}$ is replaced by the negative gust velocity $-u_{g,j}$. Thus, the equations read:

$$ \frac{d}{dt} \int_V \rho dV + \int_S \rho (\tilde{u}_j + u_{g,j}) \cdot n_j \, dS = 0, \quad (5) $$

$$ \frac{d}{dt} \int_V \rho \tilde{u}_i \, dV + \int_S \rho \tilde{u}_i (\tilde{u}_j + u_{g,j}) \cdot n_j \, dS \\ = - \int_S \tilde{\tau}_{ij} \cdot n_j \, dS - \int_S p \cdot n_i \, dS. \quad (6) $$

Even though the grid is fixed, it is still necessary to satisfy a GCL in the context of FVM due to the fact that the prescribed gust velocity component is not necessarily divergence-free. This will again enable the reduction of the mass conservation equation (i.e., Eq. (5) to Eq. (4)) and thus the pressure-correction equation still remains unaffected. Hence, it is clear that a rigorous fulfillment of the GCL will suppress any spurious source terms that might be generated. Despite the undeformed grid ($V^{\text{grid}}$ is the real volume of the cell and taken as constant) a pseudo change in the cell volumes (denoted $\Delta V = V^{*} - V^{\text{grid}}$ and referred to as apparent grid movement in Sitaraman and Baeder (2006)) has to be taken into account in order to ensure the pseudo geometric conservation law:

$$ \frac{d}{dt} \int_{V^{\text{grid}} \to V^{*}} dV + \int_{S^{*}} u_{g,j} \cdot n_j \, dS = 0. \quad (7) $$

Since the pseudo updated surface $S^*$ of each CV is not available, the pseudo geometric conservation law is rewritten applying the Gaussian integral theorem:

$$ \frac{d}{dt} \int_{V^{\text{grid}} \to V^{*}} dV + \int_{V^{*}} \frac{\partial u_{g,j}}{\partial x_j} dV = 0. $$

In a first approximation, when a second-order approximation of the volume integral such as the mid-point rule is employed in the context of a finite-volume scheme, it follows:

$$ \frac{d}{dt} \int_{V^{\text{grid}} \to V^{*}} dV + \gamma \int_{V^{\text{grid}}} \frac{\partial u_{g,j}}{\partial x_j} dV = 0, $$

where the scaling factor $\gamma$ is defined as $\gamma = \frac{V^{*}}{V^{\text{grid}}}$. A back transformation leads to:

$$ \frac{d}{dt} \int_{V^{\text{grid}} \to V^{*}} dV + \gamma \int_{S} u_{g,j} \cdot n_j \, dS = 0. \quad (8) $$

The derived pseudo geometric conservation law now only depends on the pseudo updated volume $V^*$, $S$ is the current real surface of the grid and thus denoted $S^{\text{grid}}$ in the following.

Coming back to the mass and momentum conservation equations, the contributions due to the pseudo grid movement have to be evaluated on the updated volumes $V^*$ and surfaces $S^*$, whereas all other terms related to the background velocity are computed on the real grid:

$$ \frac{d}{dt} \int_{V^{\text{grid}} \to V^{*}} \rho dV + \int_{S^{\text{grid}}} \rho \tilde{u}_j \cdot n_j \, dS + \int_{S^{*}} \rho u_{g,j} \cdot n_j \, dS = 0, \quad (9) $$

$$ \frac{d}{dt} \int_{V^{\text{grid}} \to V^{*}} \rho \tilde{u}_i \, dV + \int_{S^{\text{grid}}} \rho \tilde{u}_i \tilde{u}_j \cdot n_j \, dS + \int_{S^{*}} \rho \tilde{u}_i u_{g,j} \cdot n_j \, dS \\ = - \int_{S^{\text{grid}}} \tilde{\tau}_{ij} \cdot n_j \, dS - \int_{S^{\text{grid}}} p \cdot n_i \, dS. \quad (10) $$

Applying the methodology mentioned above to the mass and momentum conservation leads to:

$$ \frac{d}{dt} \int_{V^{\text{grid}} \to V^{*}} \rho dV + \int_{S^{\text{grid}}} \rho \tilde{u}_j \cdot n_j \, dS + \gamma \int_{S^{\text{grid}}} \rho u_{g,j} \cdot n_j \, dS = 0, \quad (11) $$

$$ \frac{d}{dt} \int_{V^{\text{grid}} \to V^{*}} \rho \tilde{u}_i \, dV + \int_{S^{\text{grid}}} \rho \tilde{u}_i \tilde{u}_j \cdot n_j \, dS + \gamma \int_{S^{\text{grid}}} \rho \tilde{u}_i u_{g,j} \cdot n_j \, dS $$
$$= - \int_{S\text{grid}} \tilde{\tau}_{ij} \cdot n_j \, dS - \int_{S\text{grid}} p \cdot n_i \, dS. \qquad (12)$$

Note that in analogy to the ALE approach, employing the GCL to the mass conservation equation (11) allows to cancel out the contribution of the pseudo-grid movement to the mass fluxes with the unsteady term such that the mass conservation equation again reduces to Eq. (4). Since the gust velocity is a prescribed velocity, the additional mass fluxes in the momentum equation can be directly evaluated.

For a flow solver already able to solve ALE problems, the integration of the FVM is straightforward. Further advantages and disadvantages are as follows. In FVM the gust is prescribed and the remaining field is solved, hence no numerical dissipation of the gust is appearing. Moreover, the gust shape and position can be prescribed at any time in the computational domain. It has the advantage that even when the amplitude of the gust is significantly larger than the background velocity, the simulation will converge. Despite the reasonable results delivered by FVM, it is not supported by a clear description of simplifications used in its derivation. Even though FVM accounts for the influence of the gust on the surrounding flow field, its major drawback lies on its failure to capture the effect of the surrounding on the gust itself (often denoted as the feedback effect). In addition, FVM was found to be accurate when solving the Euler equations. The reason is that the gust velocity is taken into account when the convective fluxes are calculated across cell faces. However, when solving the Navier-Stokes equations with a non-uniform gust velocity distribution, viscous flux correction is required. Gu et al. (2012) demonstrated this issue for accurately simulating practical flow phenomena (such as pitching and plunging motions of an airfoil). Their contribution relies on the principle of relative motion to show that the velocity which needs to be used for the calculation of the viscous fluxes must be a summation of the background and the gust velocity. They show that this correction is indispensable if the gradient of the gust velocity is different from zero. Consequently, the magnitude of the applied correction is depending on how steep the gust velocity profile is. The drawbacks of the original FVM are not existent when using the split velocity method (SVM) explained next.

## 2.2. Split Velocity Method (SVM)

The split velocity method was proposed by Wales et al. (2014). As before, the total velocity $u_i$ is decomposed into two components, namely, a background $\tilde{u}_i$ and a prescribed gust velocity $u_{g,i}$. However, now the governing equations are written on a fixed grid based on the total velocity $u_i$ and rearranged on that grid without the need for any simplifying assumptions. As explained below, in order to fulfill the postulate of $\frac{\partial \tilde{u}_i}{\partial x_i} = 0$, the ALE formulation of the Navier-Stokes equations is used once more including the necessary application of the pseudo GCL. The SVM has been applied to 2-D Euler equations (Wales et al., 2014) and to simulate 2-D as well as 3-D viscous flows (Huntley et al., 2016).

The formulation starts by writing the integral form of the Navier-Stokes equations on a fixed grid:

$$\int_{V\text{grid}} \frac{\partial \rho}{\partial t} dV + \int_{S\text{grid}} \rho u_i \cdot n_i \, dS = 0, \qquad (13)$$

$$\int_{V\text{grid}} \frac{\partial (\rho u_i)}{\partial t} dV + \int_{S\text{grid}} \rho u_i u_j \cdot n_j \, dS = - \int_{S\text{grid}} \tau_{ij} \cdot n_j \, dS - \int_{S\text{grid}} p \cdot n_i \, dS. \qquad (14)$$

$\tau_{ij}$ is the stress tensor based on the total velocity $u_i$. Substituting the decomposed total velocity $u_i$, the conservation equations read:

$$\int_{V\text{grid}} \frac{\partial \rho}{\partial t} dV + \int_{S\text{grid}} \rho (\tilde{u}_j + u_{g,j}) \cdot n_j \, dS = 0, \qquad (15)$$

$$\int_{V\text{grid}} \frac{\partial (\rho (\tilde{u}_j + u_{g,j}))}{\partial t} dV + \int_{S\text{grid}} \rho (\tilde{u}_i + u_{g,i}) (\tilde{u}_j + u_{g,j}) \cdot n_j \, dS = 0. \qquad (16)$$

$$= - \int_{S\text{grid}} \tau_{ij} \cdot n_j \, dS - \int_{S\text{grid}} p \cdot n_i \, dS. \qquad (16)$$

The terms of the momentum equation are rewritten as follows:

1. Transient term (note that in the context of fixed grids, the temporal derivative inside the volume integral can be directly replaced by a temporal derivative outside the integral):

$$\int_{V\text{grid}} \frac{\partial (\rho (\tilde{u}_i + u_{g,i}))}{\partial t} dV = \int_{V\text{grid}} \frac{\partial (\rho \tilde{u}_i)}{\partial t} dV + \int_{V\text{grid}} \frac{\partial (\rho u_{g,i})}{\partial t} dV = \frac{\partial}{\partial t} \int_{V\text{grid}} \rho \tilde{u}_i \, dV + \int_{V\text{grid}} \rho \frac{\partial u_{g,i}}{\partial t} \, dV + \underbrace{\int_{V\text{grid}} u_{g,i} \frac{\partial \rho}{\partial t} \, dV}_{\bullet}. \qquad (17)$$

2. Convection term:

$$\begin{aligned} & \int_{S\text{grid}} \rho (\tilde{u}_i + u_{g,i}) (\tilde{u}_j + u_{g,j}) \cdot n_j \, dS \\ &= \int_{V\text{grid}} \frac{\partial}{\partial x_j} (\rho (\tilde{u}_i + u_{g,i}) (\tilde{u}_j + u_{g,j})) \, dV \\ &= \int_{V\text{grid}} \frac{\partial}{\partial x_j} (\rho \tilde{u}_i (\tilde{u}_j + u_{g,j})) \, dV + \int_{V\text{grid}} u_{g,i} \frac{\partial}{\partial x_j} (\rho (\tilde{u}_j + u_{g,j})) \, dV \\ & \quad + \int_{V\text{grid}} \rho (\tilde{u}_j + u_{g,j}) \frac{\partial u_{g,i}}{\partial x_j} \, dV \\ &= \int_{S\text{grid}} \rho \tilde{u}_i (\tilde{u}_j + u_{g,j}) \cdot n_j \, dS + \underbrace{\int_{V\text{grid}} u_{g,i} \frac{\partial}{\partial x_j} (\rho (\tilde{u}_j + u_{g,j})) \, dV}_{\bullet} \\ & \quad + \int_{V\text{grid}} \rho (\tilde{u}_j + u_{g,j}) \frac{\partial u_{g,i}}{\partial x_j} \, dV. \end{aligned} \qquad (18)$$

The summation of the terms marked by • represents the product of the gust velocity $u_{g,i}$ with the mass conservation equation (13) expressed by $u_i = \tilde{u}_i + u_{g,i}$:

$$\int_{V\text{grid}} u_{g,i} \frac{\partial \rho}{\partial t} dV + \int_{V\text{grid}} u_{g,i} \frac{\partial}{\partial x_j} (\rho u_i) \, dV = \int_{V\text{grid}} u_{g,i} \left( \frac{\partial \rho}{\partial t} + \frac{\partial (\rho u_j)}{\partial x_j} \right) dV = 0. \qquad (19)$$

Hence, the momentum equation can be rewritten in its simplified form as:

$$\begin{aligned} & \frac{\partial}{\partial t} \int_{V\text{grid}} \rho \tilde{u}_i \, dV + \int_{S\text{grid}} \rho \tilde{u}_i (\tilde{u}_j + u_{g,j}) \cdot n_j \, dS \\ &= - \int_{V\text{grid}} \rho \frac{\partial u_{g,i}}{\partial t} \, dV - \int_{V\text{grid}} \rho (\tilde{u}_j + u_{g,j}) \frac{\partial u_{g,i}}{\partial x_j} \, dV \\ & \quad - \underbrace{- S_{\text{momentum},i}} \\ & \quad - \int_{S\text{grid}} \tau_{ij} \cdot n_j \, dS - \int_{S\text{grid}} p \cdot n_i \, dS. \end{aligned} \qquad (20)$$

In this formulation the convective fluxes are described by mass fluxes, which are predicted by the total velocity (i.e., the sum of the background and the gust velocity). Furthermore, the decomposition of the velocity leads to the formation of a source term. This term models the convective effect of the object being hit by the gust on the gust itself. The source term can be written in its general form as:

$$S_{\text{momentum},i} = \int_V \rho \left\{ \frac{\partial u_{g,i}}{\partial t} + (\tilde{u}_j + u_{g,j}) \frac{\partial u_{g,i}}{\partial x_j} \right\} dV. \qquad (21)$$

In the framework of fixed grids, the partial derivative in time inside the volume integral can be replaced by a total derivative outside the volume integral. Hence, the mass conservation equation (15) reads:

$$\frac{d}{dt} \int_{V\text{grid}} \rho \, dV + \int_{S\text{grid}} \rho (\tilde{u}_j + u_{g,j}) \cdot n_j \, dS = 0. \qquad (22)$$
The transient term in Eq. (22) vanishes for an incompressible fluid implying that the divergence of the total velocity field is equal to zero. Since the prescribed velocity is not necessarily divergence-free, the postulate $\frac{\partial \tilde{u}_i}{\partial x_i} = 0$ might not be fulfilled. In order to comply with this condition, the pseudo GCL (i.e., Eq. (7)) is employed on the transient term of Eq. (22):

$$ \frac{d}{dt} \int_{V\text{grid} \to V^*} \rho dV = - \int_{S^*} \rho u_{g,j} \cdot n_j dS. \qquad (23) $$

Substituting Eq. (23) into Eq. (22), Eq. (4) is retrieved. In a similar manner, the partial temporal derivative in time in the momentum equation (20) on a fixed grid is substituted by a total derivative and the ALE formulation is employed consistently:

$$ \frac{d}{dt} \int_{V\text{grid} \to V^*} \rho \tilde{u}_i dV + \int_{S\text{grid}} \rho \tilde{u}_i \tilde{u}_j \cdot n_j dS + \int_{S^*} \rho \tilde{u}_i u_{g,j} \cdot n_j dS = -S_{\text{momentum},i} - \int_{S\text{grid}} \tau_{ij} \cdot n_j dS - \int_{S\text{grid}} p \cdot n_i dS. \qquad (24) $$

Employing equation (8) deduced in the FVM section, the final governing equations solved for SVM containing the scaling factor $\gamma$ are:

$$ \frac{d}{dt} \int_{V\text{grid} \to V^*} \rho dV + \int_{S\text{grid}} \rho \tilde{u}_j \cdot n_j dS + \gamma \int_{S\text{grid}} \rho u_{g,j} \cdot n_j dS = 0, \qquad (25) $$

$$ \frac{d}{dt} \int_{V\text{grid} \to V^*} \rho \tilde{u}_i dV + \int_{S\text{grid}} \rho \tilde{u}_i \tilde{u}_j \cdot n_j dS + \gamma \int_{S\text{grid}} \rho \tilde{u}_i u_{g,j} \cdot n_j dS = -S_{\text{momentum},i} - \int_{S\text{grid}} \tau_{ij} \cdot n_j dS - \int_{S\text{grid}} p \cdot n_i dS. \qquad (26) $$

Comparing the mass conservation equation (11) for FVM with the corresponding equation (25) for SVM, both formulations are identical and both allow a reduction to the original mass conservation equation (4) by the GCL. Obviously, the differences between both approaches arise from the momentum equations. In comparison to the field velocity approach, the evaluation of the additional source terms in Eq. (21) constitutes an important aspect of the split velocity method, as these terms enable to capture the effect of the body on the wind gust. Furthermore, in contrast to FVM, SVM uses the total velocity field for the computation of the viscous fluxes. Hence, it is suitable for the simulation of viscous flow problems. So the major drawbacks of FVM are not present in SVM. However, the implementation of SVM needs more changes in the code than FVM.

Note that even in case of LES predictions the viscous terms cannot always be neglected in comparison with the resolved stresses. For flow problems with separation for example at curved walls or with transition to turbulence, the viscous stresses are decisive for the correct prediction of the location of these phenomena and thus the entire flow development.

## 2.3. Boundary conditions

At the boundary patches of a computational domain, velocities have to be assigned and convective and diffusive fluxes have to be computed depending on the boundary condition type (e.g., inlet, outlet, symmetry, periodic and wall). The distinction made between the background and gust velocities in the prescribed velocity methods may act as a potential source of misconception. Therefore, special care has to be taken on the application of the boundary conditions in the context of FVM and SVM.

1. **Velocities:** In both methods, the variable for which the momentum equations are solved is the background velocity $\tilde{u}_i$. Therefore, velocities at the faces belonging to a boundary patch are assigned based on this velocity. For example, for a stationary

<figure><img src="image_4.png" /><figcaption>Fig. 1. Local basis definition $B_1 = (g_1, g_2, g_3)$ of the gust.</figcaption></figure>

wall, the total velocity is equal to zero fulfilling the no-slip and impermeability conditions. Hence, the background flow velocity at the wall is set to $-u_{g,i}$. Similarly, the velocities at the inlet, outlet and symmetry faces are assigned based on the primary variable $\tilde{u}_i$.

2. **Convective fluxes:** In both methods, the convective fluxes should eventually include the contribution of both velocities (i.e., the background as well as the gust velocity). For example, at a no-slip wall, mass fluxes are first computed based on the assigned background velocity and then corrected using the gust velocity. Hence, this yields a mass flux fulfilling the impermeability condition of a no-slip wall. The same treatment is applied to all boundary conditions.

3. **Diffusive fluxes:** The only difference between FVM and SVM in the application of boundary conditions lies on the computation of the diffusive fluxes. While the diffusive fluxes on boundary faces (e.g., wall shear stress of a no-slip wall or the wall normal stress of a symmetry patch) are solely computed based on the background velocity in FVM, these are evaluated using the contribution of both velocities in SVM.

## 3. Definition of the gust shape

The gust shape is defined in the local orthonormal basis $B_1 = (g_1, g_2, g_3)$ depicted in Fig. 1 assuming that the only non-zero gust velocity component is that of the first basis vector of the local coordinate system (i.e., of $g_1$).

Letting $u_{g,j}$ and $\theta_j$ denote the components of the gust velocity vector in the Cartesian basis $B_0 = (e_1, e_2, e_3)$ and in the local basis $B_1 = (g_1, g_2, g_3)$, respectively, the gust velocity reads:

$$ \mathbf{u}_g = u_{g,1} \mathbf{e}_1 + u_{g,2} \mathbf{e}_2 + u_{g,3} \mathbf{e}_3 = \theta_1 \mathbf{g}_1 + \theta_2 \mathbf{g}_2 + \theta_3 \mathbf{g}_3 = \theta_1 \mathbf{g}_1. \qquad (27) $$

The first basis vector of the local coordinate system, i.e., $g_1$ is imposed by the user. $g_2$ and $g_3$ can be imposed or computed so that $(g_1, g_2, g_3)$ forms an orthonormal basis. The transformation matrix $T_{B_1 \to B_0}$ from the local basis $B_1$ to the Cartesian basis $B_0$ is required. Using the abbreviation $\beta_{ij} = \frac{\partial u_{g,j}}{\partial \theta_i}$, one could write the transformation in a matrix notation as:

$$ \begin{gathered} [g_i] = [\beta_{ij}] [e_j], \\ \underbrace{\begin{bmatrix} g_1 \\ g_2 \\ g_3 \end{bmatrix}}_{\text{old basis}} = \underbrace{\begin{bmatrix} \frac{\partial u_{g,1}}{\partial \theta_1} & \frac{\partial u_{g,2}}{\partial \theta_1} & \frac{\partial u_{g,3}}{\partial \theta_1} \\ \frac{\partial u_{g,1}}{\partial \theta_2} & \frac{\partial u_{g,2}}{\partial \theta_2} & \frac{\partial u_{g,3}}{\partial \theta_2} \\ \frac{\partial u_{g,1}}{\partial \theta_3} & \frac{\partial u_{g,2}}{\partial \theta_3} & \frac{\partial u_{g,3}}{\partial \theta_3} \end{bmatrix}}_{\text{new basis}} \underbrace{\begin{bmatrix} e_1 \\ e_2 \\ e_3 \end{bmatrix}}_{\text{old basis}} = \begin{bmatrix} \beta_{11} & \beta_{12} & \beta_{13} \\ \beta_{21} & \beta_{22} & \beta_{23} \\ \beta_{31} & \beta_{32} & \beta_{33} \end{bmatrix} \begin{bmatrix} e_1 \\ e_2 \\ e_3 \end{bmatrix}. \end{gathered} \qquad (28) $$
Using the relation deduced in Eq. (28), the definition of the gust velocity with respect to the physical Cartesian basis reads:

$$ \mathbf{u}_g = \theta_1 \mathbf{g}_1 = \theta_1 (\beta_{11} \mathbf{e}_1 + \beta_{12} \mathbf{e}_2 + \beta_{13} \mathbf{e}_3). \qquad (29) $$

In the context of deterministic gust models, the gust velocity component $\theta_1$ describes the spatial and temporal distributions of the gust using analytic functions:

$$ \theta_1(t, \xi, \eta, \zeta) = A_g f_t(t) f_1(\xi) f_2(\eta) f_3(\zeta). \qquad (30) $$

$A_g$ is the user-defined amplitude of the gust. $f_t$, $f_1$, $f_2$ and $f_3$ are user-defined analytic functions representing the shape of the gust in time and space (in $\mathbf{g}_1$, $\mathbf{g}_2$- and $\mathbf{g}_3$-direction), respectively. $\xi, \eta, \zeta$ are the spatial coordinates in the local orthonormal basis ($\mathbf{g}_1$, $\mathbf{g}_2$, $\mathbf{g}_3$). The local coordinates $\xi, \eta$ and $\zeta$ can be expressed depending on the Cartesian coordinates (see De Nayer and Breuer (2020) for more details).

Typically used deterministic functions required by Eq. (30) for the spatial and temporal distributions related to a gust are the Extreme Coherent Gust (ECG) and the Extreme Operating Gust (EOG) defined in the IEC-Standard (2002). Other options are a Gaussian distribution or the model recently proposed by Knigge and Raasch (2016). The latter was derived from turbulent flow fields computed by LES and is supposed to describe gusts more realistically. Here, solely the ECG is applied which relies on the “1-cosine” shape. The original “1-cosine” shape found in the literature (IEC-Standard, 2002) is expressed as follows:

$$ f_i(\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} \qquad (31) $$

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

$$ f_i(\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} \qquad (32) $$

The subscript $i$ is equal to 1, 2, 3 or t. The variable $\phi$ corresponds to the coordinate $\xi, \eta, \zeta$ or to the time $t$ for $i = 1, 2, 3$ or t, respectively. The constant $\phi_g$ is the user-defined central value of the gust distribution for the corresponding $\phi$ and $L_g^\phi$ is its user-defined length or time scale. Note that the present methodology allows to prescribe the wind gust as a three-dimensional instantaneous phenomenon.

# 4. Implementation into a finite-volume NS solver

## 4.1. Flow solver

The prescribed gust injection techniques mentioned above are incorporated into an incompressible Navier-Stokes solver, i.e., an enhanced and well validated (De Nayer et al., 2014; De Nayer and Breuer, 2014) version of FASTEST-3D (Durst and Schäfer, 1996; Breuer et al., 2012). The Navier-Stokes equations are discretized based on the finite-volume technique on a curvilinear, block-structured body-fitted grid with a collocated variable arrangement. The surface and volume integrals are approximated by the midpoint rule. Most flow variables are linearly interpolated to the cell faces leading to a second-order accurate central scheme. The convective fluxes are approximated by the technique of flux blending (Khosla and Rubin, 1974; Ferziger and Perić, 2002) 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. In order to avoid unwanted oscillations, the momentum interpolation technique of Rhie and Chow (1983) for non-staggered grids is applied to couple the pressure and the velocity fields.

A semi-implicit predictor-corrector scheme (Breuer et al., 2012) of second-order accuracy is used to solve the pressure-velocity coupling problem. First, the momentum equations are time-marched by a low-storage multi-stage Runge-Kutta method to obtain an intermediate velocity. Then, the 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). The whole procedure provides second-order accuracy in space and time. The solver can be applied for direct numerical simulations and large-eddy simulations (Breuer et al., 2012; De Nayer et al., 2014; De Nayer and Breuer, 2014, 2020).

## 4.2. Additional implementations for FVM and SVM

In order to implement the field and the split velocity method into a CFD code, the gust velocity $u_{g,j}$ is prescribed in the domain at the beginning of each time step. The pseudo updated volume $V^*$ is then computed and applied to the time discretization during the solution of the momentum equation. The momentum equation includes convective fluxes relying on $u_{g,j}$ and $\bar{u}_i$. For a three-dimensional finite-volume CFD solver working on curvilinear grids (see Section 4.1) these steps are detailed below:

* The gust shape is prescribed as described in Section 3.

* Based on the definition of the gust velocity, the volumetric flow rate $Q$ through the faces of a finite volume due to the gust velocity reads:

$$ Q = \theta_1 (\beta_{11} A_x + \beta_{12} A_y + \beta_{13} A_z). \qquad (33) $$

$A_x, A_y$ and $A_z$ are the face area projected on $\mathbf{e}_1, \mathbf{e}_2$ and $\mathbf{e}_3$, respectively.

* In the context of the prescribed gust methods, a volume change due to a pseudo grid movement takes place if a gust velocity component varies along its direction. The pseudo cell volume $V^*$ is computed based on the definition of the pseudo GCL according to Eq. (8):

$$ \frac{d}{dt} \int_{V_{\text{grid}} \to V^*} dV = \gamma \sum_a \int_{S_a} -u_{g,j} \cdot n_j \, dS, \\ \frac{V^* - V_{\text{grid}}}{\Delta t} = - \frac{V^*}{V_{\text{grid}}} (Q_e + Q_w + Q_n + Q_s + Q_t + Q_b), $$

which leads to:

$$ V^* = \frac{V_{\text{grid}}}{1 + \frac{\Delta t}{V_{\text{grid}}} (Q_e + Q_w + Q_n + Q_s + Q_t + Q_b)}. \qquad (34) $$

`{e,w,n,s,t,b}` are the indices of the east, west, north, south, top and bottom faces of a hexahedral control volume, respectively.

* The additional mass fluxes due to the prescribed gust methods can now be evaluated based on the previously mentioned volumetric flow rate and the pseudo updated volume. For example, for the east face the additional mass flux $\dot{m}_e$ reads:

$$ \dot{m}_e = \int_{S_e^*} \rho u_{g,j} \cdot n_j \, dS \approx \rho \underbrace{\frac{V^*}{V_{\text{grid}}}}_{\gamma} \underbrace{\theta_1 (\beta_{11} A_{x,e} + \beta_{12} A_{y,e} + \beta_{13} A_{z,e})}_{Q_e \text{ from Eq. (33)}} = \rho \gamma Q_e. \qquad (35) $$

* Similar to moving grids, the pseudo change in volume has to be also included in the integration scheme in time. In FASTEST-3D, a three-substeps second-order explicit Runge-Kutta time integration scheme is employed to march the solution in time. For a conserved quantity the momentum equation is integrated with
respect to time as follows:

$$ 
\begin{align} 
& \int_{t^n}^{t^{n+1}} \frac{d}{dt} (\rho \phi V) dt = \int_{t^n}^{t^{n+1}} \text{RHS} dt, \nonumber \\ 
& [\rho \phi V]^{n+1} - [\rho \phi V]^n = \Delta t \text{RHS}(\phi^n), \tag{36} \\ 
& [\rho \phi]^{n+1} V^* - [\rho \phi]^n V^{\text{grid}} = \Delta t \text{RHS}(\phi^n), \nonumber \\ 
& \phi^{n+1} = \left( V^{\text{grid}} \phi^n + \frac{\Delta t}{\rho} \text{RHS}(\phi^n) \right) \frac{1}{V^*}. \nonumber 
\end{align} 
$$

Incorporating Runge-Kutta coefficients ($\alpha_k, k = 1, 2, 3$), the conserved quantity is advanced from the time step $t^n$ to $t^{n+1}$ in substeps by:

$$ \phi^k = \left( V^{\text{grid}} \phi^n + \alpha_k \frac{\Delta t}{\rho} \text{RHS}(\phi^{k-1}) \right) \frac{1}{V^*}. $$

Following the work of Münsch (2015), a linear variation of the volume of the cell is assumed within the Runge-Kutta substeps:

$$ \phi^k = \left( V^{\text{grid}} \phi^n + \alpha_k \frac{\Delta t}{\rho} \text{RHS}(\phi^{k-1}) \right) \frac{1}{V^{\text{grid}} + \alpha_k (V^* - V^{\text{grid}})}. $$

The momentum equation is solved for the background velocity $\tilde{u}_i$. Replacing the quantity $\phi$ by $\tilde{u}_i$ delivers the final time discretization scheme:

$$ \tilde{u}_i^k = \left( V^{\text{grid}} \tilde{u}_i^n + \alpha_k \frac{\Delta t}{\rho} \text{RHS}(\tilde{u}_i^{k-1}) \right) \frac{1}{V^{\text{grid}} + \alpha_k (V^* - V^{\text{grid}})}. \quad (37) $$

* The prescribed gust methods can be coupled with a temporally varying domain in case of FSI. In the case of moving grids, the temporally discretized momentum equation is (Breuer et al., 2012):

$$ \tilde{u}_i^k = \left( V^n \tilde{u}_i^n + \alpha_k \frac{\Delta t}{\rho} \text{RHS}(\tilde{u}_i^{k-1}) \right) \frac{1}{V^n + \alpha_k (V^{n+1} - V^n)}. $$

where $V^n$ and $V^{n+1}$ are the geometric cell volumes corresponding to the deformed grid at $t^n$ and $t^{n+1}$, respectively. The insertion of the prescribed gust methods relies on the addition of the pseudo volume change to the geometric cell volumes of the current deformed grid, i.e., $V^{n+1}$. Hence, the solution is now marched in time using the updated cell volume (i.e., $(V^{n+1})^*$ evaluated by Eq. (34)) as:

$$ \tilde{u}_i^k = \left( V^n \tilde{u}_i^n + \alpha_k \frac{\Delta t}{\rho} \text{RHS}(\tilde{u}_i^{k-1}) \right) \frac{1}{V^n + \alpha_k ((V^{n+1})^* - V^n)}. \quad (38) $$

The implementation of the additional source terms required for SVM is straightforward and thus needs no further explanation.

## 5. Test case: Flat plate with rounded edges

The aim of the present test cases originally defined in Biler et al. (2019), Badrya and Baeder (2019) and Badrya et al. (2021) is to compare the impact of gusts prescribed by FVM or SVM on the aerodynamic forces of a fixed rigid body. To allow a direct comparison with some of the results of these previous studies, the geometry of the flat plate, the Reynolds number and the properties of the vertical gusts are taken over. However, since different numerical methods (structured overset grid vs. curvilinear, block-structured body-fitted grid) are applied, the computational domain and the numerical grid have to be adjusted accordingly.

### 5.1. Setup and initial conditions

The geometry of the body is taken as simple as possible: A thin flat plate (chord $c$) with a small thickness ($h = 3 \times 10^{-2} c \ll c$) and rounded front and rear edges (half-cylinder with a radius of $h/2$) is oriented along the horizontal axis $e_1 = e_x$ as depicted in Fig. 2. Rounded edges have the advantage that in comparison to a sharp edge the location of separation is not fixed, which renders the test case more challenging

<figure><img src="image_5.png" /><figcaption>Fig. 2. Geometry of the test case and definition of the boundary conditions.</figcaption></figure>

<table><caption>Table 1
Geometrical setup of the flat plate case.</caption><tr><td>Plate thickness</td><td>$h/c = 3 \times 10^{-2}$</td></tr><tr><td>Plate width</td><td>$B/c = 4$</td></tr><tr><td>Domain length</td><td>$D/c = 14$</td></tr><tr><td>Domain height</td><td>$H/c = 13.5$</td></tr><tr><td>Domain width</td><td>$W/c = 14$</td></tr></table>

for the prediction of viscous flows. The separation point depends on the Reynolds number chosen, which is not the case for the sharp counterpart. The plate has a finite width of $B = 4c$. The side edges of the plate are sharp. The x-y-cross-section of the computational domain is ovally shaped with its center at the center of the plate. Its large radius oriented in streamwise direction is equal to $R = 7c$ and its small radius oriented in wall-normal direction is equal to $R = 6.75c$. Thus, the boundaries of the domain are far away from the plate to avoid spurious numerical effects in the far-field when gusts encounter the body. Note that the distances are above the recommendations provided by classical best practice guidelines (Franke and Baklanov, 2007; Franke et al., 2011). In spanwise direction the domain has a width of $14c$. The plate is in the middle and two spanwise side domains with a width of $5c$ are present. Thus, the recommendations provided by classical best practice guidelines are also satisfied in the lateral direction.

In the far-field the flow has a uniform inlet velocity ($u_1^{\text{inlet}} = u_\infty, u_2^{\text{inlet}} = 0, u_3^{\text{inlet}} = 0$). No-slip and impermeability boundary conditions are applied on the surface of the plate. Symmetry boundary conditions are assumed on the lateral sides of the domain in the spanwise direction. In order to avoid reflections and disturbances at the outlet, a convective boundary condition (Breuer, 2002) is defined on the outlet patch with a convective velocity ($u_1^{\text{conv}} = u_\infty, u_2^{\text{conv}} = 0, u_3^{\text{conv}} = 0$).

The geometrical and physical characteristics of the test case are summarized in Table 1. Based on the chord length $c$ of the plate, the far-field velocity $u_\infty$, the density $\rho$ and the dynamic viscosity $\mu$ of the fluid, the Reynolds number $Re_c = \rho c u_\infty / \mu$ is set to 20,000.

One of the important aspects that has to be taken into account when generating a grid for such a test case is to guarantee a reasonably fine grid near the plate to capture the flow details near the boundary. These viscous effects are significant and lead to the formation of thin boundary layers at the wall. In order to quantify the maximum cell size in the vicinity of the wall, the flat-plate boundary layer theory is employed. For a Reynolds number of $Re_c = 20,000$, the flow is assumed to be laminar. Hence, the skin friction coefficient at the trailing edge
<figure><img src="image_6.png" /><figcaption>Fig. 3. x-y plane of the computational O-type grid. Only every second grid line in each direction is shown.</figcaption></figure>

of the plate can be evaluated to $c_f = 0.664/\sqrt{Re_c} = 4.7 \times 10^{-3}$. The wall shear stress at this position is equal to $\tau_w = c_f \rho u_\infty^2/2$ and the friction velocity normalized by the free-stream velocity is given by $u_\tau/u_\infty = \sqrt{\tau_w/(\rho u_\infty^2)} = \sqrt{c_f/2} = 4.85 \times 10^{-2}$. Hence to resolve the near-wall region by a few grid points one could approximate the maximum height of the first grid cell off the wall ($\Delta y/c = 2 \times 0.2/(Re_c \sqrt{c_f/2}) \approx 4 \times 10^{-4}$). Therefore, a dimensionless spacing of $4 \times 10^{-4}$ is used in the wall-normal direction and the grid is stretched according to a geometric series with a mild stretching factor of 1.05. A zoom of the O-type grid in the x-y plane in the vicinity of the flat plate is shown in Fig. 3. Note that at the leading edge the skin friction theoretically rises to infinity. An estimation at a distance of 5% chord length from the leading edge shows that even in this region the current wall-normal resolution satisfies the condition that the first cell center is below $y_{cc}^+ = 1$.

An equidistant grid spacing of $2 \times 10^{-2}c$ is assigned in the spanwise direction to the cells of the middle block (4c) above and below the plate to capture the three-dimensional flow structures during the gust encounter. To guarantee a smooth grid transition on the side edges in the spanwise direction of the other blocks, a first cell size of $2 \times 10^{-2}c$ and a stretching factor of 1.05 are set. A summary of the fine grid on which the simulations are performed is provided in Table 2. In addition to this grid, two coarser grids are applied for a grid convergence study (see Appendix B).

To ensure the stability of the explicit time integration scheme used in the following simulations, a dimensionless time step size of $\Delta t u_\infty/c = 2 \times 10^{-4}$ is applied. That leads to a CFL number of 0.35, where in case of an incompressible fluid CFL is based on the characteristic fluid velocity $u_\infty$ instead of the speed of sound.

The initial flow conditions are taken from the quasi steady-state flow solution in order to suppress the effects of the flow development stage. The drag and lift coefficients are chosen as the parameters indicating whether the simulation reaches its quasi steady-state solution. The dimensionless force coefficients are defined as the force components normalized by the product of the dynamic pressure (i.e., $\frac{1}{2} \rho u_\infty^2$) and the area of the plate (i.e., $A = c B$):

$$C_D = \frac{F_D}{\frac{1}{2} \rho u_\infty^2 A}, \quad C_L = \frac{F_L}{\frac{1}{2} \rho u_\infty^2 A} \quad (39)$$

The simulation was carried out for a total dimensionless time interval of $t^* = t u_\infty / c = 10$ starting from initial conditions given by a parallel flow with $u = u_\infty$. The force coefficients recorded are shown in Fig. 4.

While the drag coefficient has converged to a value of about 0.025, the lift coefficient is fluctuating around a mean value of zero. This was expected as a result of the flow separation and the resulting vortex shedding taking place at the trailing edge of the plate depicted by the flow fields shown in Fig. 5. Note that the amplitudes of the lift variation are small compared to those expected by the vertical gusts.

## 5.2. Vertical gusts

The first configuration is devoted to the modeling of gusts having their direction perpendicular to the streamwise direction of the main flow. Hence, the name vertical gust is used in the present work. Note that wind gusts are in general three-dimensional structures. As mentioned above, the methodology described in Section 3 allows to prescribe the gust in this manner. However, for fundamental investigations canonical one-dimensional wind gusts are more suitable and are thus considered here.

A generic vertical gust configuration allows the evaluation of the rapid increase in lift forces when a structure is subjected to a vertical wind gust. These gusts are widely used in the literature (Sitaraman, 2003; Zaide and Raveh, 2006; Heinrich and Reimer, 2013; Wales et al., 2014; Ghoreyshi et al., 2018; Biler et al., 2019; Badrya et al., 2021) to investigate the response of an airfoil penetrating through a gust and consequently facilitating the design optimization process. In these investigations, a variety of gust shapes is employed. However, the current work is restricted to the application of gusts with an ECG shape not varying in time (see Heinrich and Reimer (2013), Biler et al. (2019) and Badrya et al. (2021)).

Based on the definition of the gust introduced in Section 3, the basis vector $g_1$ (defining the gust direction) has to be oriented towards the basis vector $e_2$. The gust velocity component $\theta_1$ given by Eq. (30) is defined as follows. The amplitude $A_g$ is set to 0.8 as in Badrya and Baeder (2019) and Badrya et al. (2021). Since the gust shape is independent of time, because it is purely convected downstream, the temporal shape function $f_t(t)$ is a uniform value (i.e., $f_t(t) = 1$). Along the $g_1$ and $g_3$ directions (i.e., the vertical and spanwise directions, respectively), the gust is uniform. Hence, a uniform value is assigned to $f_1(\xi)$ and $f_3(\xi)$. Along the $g_2$ direction (i.e., the streamwise direction in the case of a vertical gust), the ECG shape is set ($f_2(\eta) = ECG$). The center $\eta_g$ of the ECG shape is convected downstream in time. In this vertical configuration, the displacement of the gust in the streamwise direction is prescribed by the convective velocity $u_{conv} = u_\infty$ leading to $\eta_g = \eta_g^0 + (t - t_0) u_{conv}$. At $t = t_0$, the gust is prescribed at the center $\eta_g^0$. Different gust wavelengths $L_g^\eta$ in $g_2$-direction are evaluated. The value of $\eta_g^0$ is chosen so that the front of the widest gust, i.e., $L_g^\eta = 4c$, starts just before the leading edge of the plate. Thus, if the origin of the frame is taken at the leading edge, $\eta_g^0 = -2.1c$. The gust characteristics of the simulations performed are summarized in Table 3.

The length scale of the gust is chosen between 0.5 and 4 times the chord length $c$ of the plate. As mentioned in the introduction, in the intermediate run investigation of the effect of wind gusts on flexible (mechanical and civil engineering) structures are intended. Thus, this is the range of length scales which might be most important for the structural integrity and safety of the design. If the scales of the wind gusts are much smaller, the length scales can be considered as a part of the spectrum of the turbulent flow. If the scales are much larger, the corresponding flow fields are seen as long-term variations of the flow which are less critical than short-term highly instantaneous loads. The latter may lead to high tensile stresses, which might exceed the yield stress of the material and thus are critical for the safety. That explains the range of the investigated length scales.

## 5.3. Horizontal gusts

The second configuration is devoted to the modeling of gusts (also with an ECG shape) having their direction parallel to the streamwise direction of the main flow. Such a gust configuration is often experienced in engineering applications. For example, civil engineers are concerned with the effect of horizontal gusts on the integrity of flexible structures as well as on the stability of bridges and high-rise buildings. Moreover, energy harvesting wind turbines are also vulnerable to such gusts and this has consequently drawn the attention of the wind energy sector. A practical application in the field of aeronautics, where such a test case
<table><caption>**Table 2**
Grid information for the fine grid.</caption><thead><tr><th>Grid</th><th>Total no. of CVs</th><th>CVs along<br>streamwise<br>direction</th><th>CVs in<br>wall-normal<br>direction</th><th>CVs in<br>spanwise<br>direction</th></tr></thead><tbody><tr><td>Middle block</td><td>27,014,400</td><td>1008</td><td>134</td><td>200</td></tr><tr><td>Side blocks</td><td>2 × 7,934,400</td><td>-</td><td>-</td><td>50</td></tr></tbody></table>

<figure><img src="image_7.png" /><figcaption>**Fig. 4.** Time history of the drag and lift coefficients.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Fig. 5.** Snapshot of the quasi-steady flow fields.</figcaption></figure>

<table><caption>**Table 3**
Vertical gust characteristics.</caption><thead><tr><th rowspan="2">g<sub>1</sub></th><th rowspan="2">A<sub>g</sub></th><th rowspan="2">f<sub>t</sub>(t) = 1</th><th rowspan="2">f<sub>1</sub>(ξ) = 1</th><th colspan="2">f<sub>2</sub>(η)=ECG</th><th rowspan="2">f<sub>3</sub>(ζ) = 1</th></tr><tr><th>η<sup>0</sup><sub>g</sub></th><th>L<sup>η</sup><sub>g</sub></th></tr></thead><tbody><tr><td>e<sub>2</sub></td><td>0.8</td><td>-</td><td>-</td><td>-2.1 c</td><td>0.5 c</td><td>1 c</td><td>1.5 c</td><td>2 c</td><td>4 c</td><td>-</td></tr></tbody></table>

is critical, can be found in high-altitude and long-endurance aircraft employing high aspect-ratio wings. Such wings are susceptible to flow-induced vibrations (i.e., flutter) and associated instabilities. Subjecting such a geometry to a horizontal gust may amplify the streamwise velocity to levels exceeding the critical flutter speed.

Here, the vector $g_1$ is oriented towards the streamwise direction $e_1$. The horizontal gust was chosen to vary solely in the streamwise direction and the computation of the gust velocity component $\theta_1$ was as follows. The same amplitude $A_g$ as for the previous configuration is used (i.e., $A_g = 0.8$). Again, no gust shape variation in time is allowed. Therefore, the temporal shape function $f_t(t)$ is set to a uniform value (i.e., $f_t(t) = 1$). Along the $g_1$ direction (i.e., the streamwise direction), the gust varies with an ECG shape function ($f_1(\xi) = ECG$). The center coordinate $\xi_g$ of the ECG shape is convected downstream in time. In this configuration, the displacement of the gust in the streamwise direction is prescribed by the convective velocity $u_{conv} = u_{\infty}$ following $\xi_g = \xi_g^0 + (t - t_0) u_{conv}$. At $t = t_0$, the gust is introduced at the gust center coordinate $\xi_g^0$. Different gust wavelengths $L_g^\xi$ in $g_1$ direction are investigated (see Table 4). As for the vertical gust, $\xi_g^0$ is set so that the front of the widest gust, i.e., $L_g^\xi = 4c$, starts just before the leading edge of the plate. Along the $g_2$ and $g_3$ directions (i.e., the vertical and spanwise directions, respectively), uniform values $f_2(\eta) = 1$ and $f_3(\zeta) = 1$ are used. As motivated in Section 5.2, the length scale of the gust is again chosen between 0.5 and 4 times the chord length $c$ of the plate.

## 6. Results: FVM vs. SVM

In the following, the results of both methods are evaluated based on the vertical as well as the horizontal gust configurations. In a first step, the flow field of one specific setup with an intermediate gust length of $2c$ is visualized at different characteristic instants in time in order to illustrate the flow development when the gust encounters the plate. This also allows to work out deviations observed between the predictions based on FVM and SVM and to illuminate the resulting time histories of the lift and drag coefficients. In a second step, the full range of gust lengths is exploited and the resulting time histories of the force coefficients are evaluated.

### 6.1. Vertical gusts

#### 6.1.1. Flow development at $L_g^\eta = 2c$

Fig. 7 depicts the flow field based on the vertical background velocity $\tilde{u}_2/u_\infty$, the total vertical velocity $u_2/u_\infty$ and the pressure normalized by $\rho u_\infty^2$. Six instants in time are chosen to present the development of the flow. At the first snapshot (see Fig. 6, time 1) the vertical gust is still in front of the plate and the ECG shape just reaches the leading edge of the plate. As expected in both predictions the velocity field is the same as shown in Fig. 5 for the quasi-steady initial conditions. At the second snapshot (see Fig. 6, time 2) the center of the gust is at the leading edge of the plate as visible in the contours of the total vertical velocity. The gust hits the plate from below leading to a flow separation
<table><caption>**Table 4**
Horizontal gust characteristics.</caption><thead><tr><th rowspan="2">g<sub>1</sub></th><th rowspan="2">A<sub>g</sub></th><th rowspan="2">f<sub>1</sub>(t) = 1</th><th colspan="2">f<sub>1</sub>(ξ)=ECG</th><th rowspan="2">f<sub>2</sub>(η)=1</th><th rowspan="2">f<sub>3</sub>(ξ) = 1</th></tr><tr><th>ξ<sup>0</sup><sub>g</sub></th><th>L<sup>ξ</sup><sub>g</sub></th></tr></thead><tbody><tr><td>e<sub>1</sub></td><td>0.8</td><td>-</td><td>-2.1 c</td><td>0.5 c</td><td>1 c</td><td>1.5 c</td><td>2 c</td><td>4 c</td><td>-</td><td>-</td></tr></tbody></table>

on the top and the development of a strong vortical clockwise rotating structure as visible in the pressure minimum in the vicinity of the leading edge. At this stage already some first small deviations between the results achieved by FVM and SVM are obvious. The third instant in time (see Fig. 7, time 3) is the one at which the plate experiences the largest lift force. The center of the ECG shape is at a distance of about 0.2 c from the leading edge. A second pressure minimum indicating a further strong clockwise rotating vortical structure can be seen close to the leading edge. Again small deviations between the results of FVM and SVM are apparent, but the differences have not grown visibly. At the fourth snapshot (see Fig. 7, time 4) the center of the gust is located at the center of the plate. The two leading-edge vortices (LEV) are convected downstream in comparison to time 3. At this phase the lift force is already decreasing again. At the fifth instant in time (see Fig. 8, time 5) the center of the gust reaches the trailing edge of the plate. Now the first LEV is located further downstream and its strength has significantly decreased, whereas the second vortical structure is hardly visible in the FVM prediction but still apparent in the SVM results. At the last snapshot (see Fig. 8, time 6) the gust already detaches from the plate completely. Now the deviations between the FVM and SVM predictions are the largest. Nevertheless, the deviations observed between FVM and SVM are rather small. This is especially true for the pressure field. The reason for this behavior is the nearly vanishing source term in the SVM prediction in case of the vertical gust (see Appendix A.1). As will be shown below, the source term is not disappearing in the horizontal gust case and thus significant differences between FVM and SVM results occur. Furthermore, in case of the vertical gust no unphysical flow phenomena are observed either for FVM nor for SVM.

### 6.1.2. Force coefficients and effect of gust length

Beside the focus of this configuration to illustrate the effect of the gust width on the forces of the flat plate, it was also used to validate the implementation of SVM in FASTEST-3D based on the results in the literature (Badrya and Baeder, 2019; Badrya et al., 2021). Fig. 9 shows the agreement achieved between the aerodynamic responses predicted by SVM using FASTEST-3D and those from the literature. Here, the aerodynamic responses of the plate are plotted as a function of the dimensionless time $t^*$.<fnref n="1" /> Small deviations (see Table 5) are observed among the responses and this could be related to the use of different discretization schemes of the underlying equations. For instance, Badrya and Baeder (2019) and Badrya et al. (2021) employed the Monotonic Upstream-centered Scheme for Conservation Laws (MUSCL) for the computation of the inviscid terms. The use of flux limiters makes the solution total variation diminishing and, therefore, spurious oscillation free. On the contrary, FASTEST-3D relies on a second-order accurate central scheme with 3% flux blending of a first-order accurate upwind scheme. Such a low-dissipative scheme is favored since it does not damp out small-scale structures typically occurring in large-eddy simulations which are the final objective of this study. Obviously, such small-scale structures develop in the present simulations leading to local maxima for example in the time history of the lift coefficient. Despite the small deviations between the present results and the literature data (see Table 5), a good arrangement is achieved for the purpose of SVM validation.

At $t^* = 0$, the plate-gust interaction commence at the leading edge of the plate yielding an increase in the lift and a drop in the drag coefficients. While the rise in the lift coefficients is a consequence of an increase in the effective angle of attack of the plate, the drop in the drag coefficients is due to the formation of a clockwise rotating leading-edge vortex (LEV). The elevation in the lift coefficients continues until a peak value is reached. Then, these coefficients start decaying to eventually recover the steady-state values when the gusts have completely left the plate. On the other hand, drag coefficients keep on decreasing until the leading edge vortex has reached its largest size. Finally, the steady-state values of the drag coefficient are retrieved following the convection of the leading-edge vortex downstream along the chord length of the plate.

In Fig. 10 the aerodynamic responses of the plate due to the gust encounter predicted by FVM are plotted against those delivered by SVM. One could notice that only some minor deviations are present between the results of the two methods (see Table 5). These deviations can be detected in the lift coefficients of gusts with short wavelengths (e.g., 0.5 c, 1 c and 1.5 c), where FVM predicts a second local peak when the gusts have entirely left the plate.

Based on the results found, one could deduce three main effects of the gust width on the aerodynamic responses appearing in FVM and SVM:

1. The rate at which the forces change: It could be observed that the shorter the gust is, the higher is the rate at which the forces increase. This could be related to the fact that a gust with a shorter wavelength has larger gradients along its profile, and therefore, leads to higher force rates.

2. The peak forces recorded: The instantaneous lift coefficient steadily increases to reach a maximum value which is depending on the gust width. In fact, the maximum lift coefficient non-linearly increases with the gust width as found in (Biler et al., 2019) and shown in Fig. 11(a), where the slope decreases for larger gust lengths. This behavior is expected since for a wider gust, the plate experiences a higher effective angle of attack along the chord for a longer period of time. Fig. 11(a) also shows that the current results of both SVM and FVM are in close agreement with the predicted data by Badrya et al. (2021) (see Table 5).

3. The instant in time at which the peaks are recorded: It was found that the peak lift occurs when the center of the gust is located between the leading edge and the mid-chord of the plate. Since wider gusts require more time to reach this region, their resulting maximum peak values occur at later times. It was also found in Biler et al. (2019) that the convective time associated with the maximum lift linearly increases with the gust width as shown in Fig. 11(b). Again the present results are found to be in good agreement with the predictions in Badrya et al. (2021) (see Table 5).

### 6.2. Horizontal gusts

#### 6.2.1. Flow development at $L_g^\xi = 2c$

The flow fields induced by a prescribed horizontal gust of $L_g^\xi = 2c$ are depicted in Figs. 12 and 13 at six instants in time using the streamwise background velocity $\bar{u}_1/u_\infty$, the streamwise total velocity $u_1/u_\infty$ and the pressure normalized by $\rho u_\infty^2$. Fig. 14 highlights the differences in the results between both methods at the same time

<footnote n="1">¹ $t^* = u_\infty t/c$ is the convective time defined in terms of the number of chords traveled by the gust.</footnote>
<figure><img src="image_9.png" /><figcaption>Fig. 6. Vertical gust passing around the rounded flat plate: Predictions based on FVM (left) and SVM (right) at time 1 and 2.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>Fig. 7. Vertical gust passing around the rounded flat plate: Predictions based on FVM (left) and SVM (right) at time 3 and 4.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>Fig. 8. Vertical gust passing around the rounded flat plate: Predictions based on FVM (left) and SVM (right) at time 5 and 6.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>**Fig. 9.** Vertical gust passing around the rounded flat plate: Instantaneous drag and lift coefficients. The gust shape is ECG with a gust amplitude of $A_g = 0.8$: FASTEST-3D vs. Badrya and Baeder (2019) SVM results.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Fig. 10.** Vertical gust passing around the rounded flat plate: Instantaneous drag and lift coefficients. The gust shape is ECG with a gust amplitude of $A_g = 0.8$: FVM vs. SVM.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Fig. 11.** Vertical gust passing around the rounded flat plate: Lift coefficient analysis: (a) Maximum lift coefficient as a function of the gust wavelength and (b) Convective time corresponding to the maximum lift coefficient as a function of the gust wavelength.<fnref n="2" /></figcaption></figure>

ininstants. Please note that the current six time instants are partially different from those used for the vertical gust configuration.

At the first snapshot (time 1) the horizontal gust has reached the front of the plate and the first quarter of the gust has already penetrated the leading edge of the plate. The predictions (FVM and SVM) show that the streamwise background velocity field and the streamwise total velocity field (see Fig. 12, time 1) are very similar to each other

<footnote n="2">² For FVM and a gust wavelength of 1.5c the second local maximum is used.</footnote>
<table><caption>**Table 5**

Normalized root-mean-square error of the lift and drag coefficients found in Fig. 9 (present simulation vs. Badrya and Baeder (2019)) and Fig. 10 (FVM vs. SVM). The relative errors are normalized by the corresponding ranges (maximum–minimum).</caption><thead><tr><th>Case</th><th>Normalized RMSE [%]</th><th>0.5 c</th><th>1.0 c</th><th>1.5 c</th><th>2.0 c</th><th>4.0 c</th></tr></thead><tbody><tr><td rowspan="2">Present sim. vs. B&B</td><td>ΔC<sub>D</sub></td><td>5.0</td><td>4.9</td><td>6.3</td><td>4.9</td><td>4.7</td></tr><tr><td>ΔC<sub>L</sub></td><td>5.0</td><td>4.6</td><td>4.6</td><td>2.8</td><td>1.3</td></tr><tr><td rowspan="2">FVM vs. SVM</td><td>ΔC<sub>D</sub></td><td>3.5</td><td>4.6</td><td>4.4</td><td>4.3</td><td>2.4</td></tr><tr><td>ΔC<sub>L</sub></td><td>3.5</td><td>3.3</td><td>2.2</td><td>3.3</td><td>3.5</td></tr></tbody></table>

and also very similar to the quasi-steady flow shown in Fig. 5. However, strong differences are visible in the pressure distribution (see Fig. 14(a)). In case of SVM an area of negative pressure forms at the prescribed gust location<fnref n="3" />. Obviously, the observed pressure distribution is due to the contribution of the SVM source term (see Appendix A.2 for further details), since a similar observation cannot be made for the FVM prediction, which does not use such a source term. A region with a strong pressure drop at the location of the gust is reasonable, since an increase of the total streamwise velocity occurs. Therefore, the pressure sinks. At the second snapshot (Fig. 12, time 2) the center of the gust is exactly at the leading edge of the plate. A minor augmentation of the background velocity is observed at the leading edge. This acceleration of the flow is more pronounced in the SVM results. The previously mentioned area of the pressure drop follows the gust prescribed by SVM. For FVM a locally occurring pressure reduction is observed near the leading edge, where the flow is locally accelerated during the deflection around the leading edge. That is also the case in the SVM prediction. However, it is not visible due to the presence of the large region of negative pressure. Contrary to the vertical gust configuration the gust hits a much smaller area of the plate. Therefore, no strong vortical structure is shed in the vicinity of the front of the plate. At the third instant in time (Fig. 12, time 3) the center of the gust is located at the center of the plate. The largest deviations are again found in the pressure. Assuming that the whole area of negative pressure following the gust observed in SVM is physical (velocity increases, thus the pressure decreases), the pressure field connected with FVM is more doubtful. The gust velocity prescribed by FVM does not directly influence the pressure. Only two negative pressure areas are observed at both plate ends, where the flow is accelerated due to the deflection induced by the plate. In the middle of the structure, where the gust is the strongest, only a mild decrease of the pressure is visible. Contrary to time 1 and 2, the streamwise background velocity field predicted by FVM and SVM starts to differ in the vicinity of the plate and in the wake (see Fig. 14(c)). The additional source term and the evaluation of the shear stress tensor based on the total velocity are the two decisive differences characterizing SVM.

In order to track which of these two issues is responsible for the difference observed in the background velocity, an altered SVM simulation, where the shear stress tensor is evaluated based on the background velocity (as for FVM) instead of the total velocity, is carried out. Since the streamwise background velocity field predicted by SVM and the altered SVM are in close agreement, it can be concluded that the SVM source term is responsible for the change in the streamwise background velocity field near the plate. The SVM source term generates the previously mentioned pressure drop area, which implies an acceleration of the flow in the tail of the gust and a deceleration in front of the gust (particularly visible near the plate in Fig. 14(c)).

The gust reaches the trailing edge of the plate in the fourth snapshot (Fig. 13, time 4). Looking at the pressure field, where the changes are most obvious, it can be noticed that the vortex shedding at the trailing edge becomes stronger in both simulations. This can be further

observed at the fifth instant in time (Fig. 13, time 5), when the plate is exposed to the last quarter of the gust. In this snapshot the background velocity predicted by FVM strongly differs from the SVM results (see Fig. 14(e)). Several short local areas of increased streamwise velocity connected with negative pressure areas are present in the vicinity of the plate in case of FVM (see Fig. 13, time 5). These small vortices are located near the plate in the tail of the gust. These strong unphysical vortical structures are neither present in SVM nor in the altered SVM. The physical explanation is that the acceleration of the flow field generated by the SVM pressure drop area induced by the source term in the tail of the gust avoids the formation of these vortices. At the last snapshot (Fig. 13, time 6) the gust detaches from the plate completely. In case of SVM the flow looks very similar to the quasi-steady case (see Fig. 5). However, the unphysical vortices are further present in the FVM simulation (see Fig. 13, time 6 and Fig. 14(f)).

In summary, the test case of the horizontal gust reveals two findings, which are not so clearly visible in the vertical gust case. These are directly related to the differences in the formulation of the FVM and SVM methods. First, the additional source term in SVM is responsible for a meaningful forecast of the pressure field and consequently of the background velocity field near the structure (feedback effect). Second, the difference in the formulation of the viscous shear stresses between FVM and SVM does not have a major effect on the predictions, at least for the present configuration. This could be related to the fact that the diffusive contribution of the gust velocity (i.e., $\tau_{ij}(u_{g,j})$) is depending on two main parameters, namely, the gust profile and the reciprocal of the Reynolds number. In the present case, the product of the spatial derivative of the gust profile with the reciprocal of the Reynolds number has an order of magnitude which is significantly smaller than that of the SVM source term.

### 6.2.2. Force coefficients and effect of gust length

Similar to the vertical gust case Fig. 15 depicts the time histories of the drag and lift coefficients for the case of the horizontal gust. To ensure a better visibility of the results, separate subfigures are provided for each gust length. At the dimensionless time instant $t^* = 0$ the branch of the gust reaches the leading edge of the plate. As visible from the drag coefficient depicted in Fig. 15(a), the duration of the gust-plate interaction lasts slightly longer than about $1 + L_g^c/c$ convective time units, which is the time the gusts need to travel completely over the plate. Obviously, the deviations observed between the drag coefficients predicted by FVM and SVM are significantly larger than for the vertical gust case. Furthermore, the lift coefficient has obviously less importance. However, it is directly connected to the van Kármán vortex shedding occurring at the trailing edge. The frequency of the vortex shedding rises during the passage of the gust over the plate. Assuming a constant Strouhal number of the shedding process the rise of the frequency can be explained by the increase of the velocity due to the horizontal gust.

In order to better understand the time history of the drag coefficient, the pressure part of the drag coefficient $C_D^p$ and the shear part of the drag coefficient $C_D^s$ are separately plotted in Fig. 16 for FVM and SVM for the gust length of $L_g^c = 2c$. It is obvious that the temporal evolution of the drag coefficient $C_D$ predicted by FVM mainly corresponds to the temporal evolution of its shear contribution. On the contrary, the curve of $C_D$ obtained by SVM is directly connected to the pressure part. This

<footnote n="3">³ Note that the reference point of the pressure is located just after the inlet patch at the height of the plate. Thus, a negative pressure has to be interpreted as the difference to the pressure at the reference point.</footnote>
<figure><img src="image_15.png" /><figcaption>Fig. 12. Horizontal gust passing around the rounded flat plate: Predictions based on FVM (left) and SVM (right) at time 1, 2 and 3.</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>Fig. 13. Horizontal gust passing around the rounded flat plate: Predictions based on FVM (left) and SVM (right) at time 4, 5 and 6.</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>Fig. 14. Horizontal gust passing around the rounded flat plate: Flow differences between FVM and SVM at six instants in time.</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>Fig. 15. *Horizontal gust passing around the rounded flat plate:* Instantaneous drag and lift coefficients for different gust lengths.</figcaption></figure>

fits to the observations made on the flow fields in the previous section. The SVM source term significantly influences the pressure distribution at the gust location and this has a strong impact on the temporal evolution of the pressure part of the drag coefficient.

At $t^* = 1$, the center of the gust reaches the leading edge. Consequently, the maximum pressure drop is also located at the leading edge. Therefore, a decrease in the drag coefficient is observed at this time instant. At $t^* = 2$, the center of the gust reaches the trailing edge and the pressure drop is located at this end of the plate. Thus, a maximum of the drag coefficient occurs at this instant in time. Of further interest is also the temporal evolution of the shear stresses. In case of SVM the shear stress part of the drag coefficient looks like the shape of the gust. On the contrary, the time history of $C_D^r$ obtained by FVM reveals that the shear part of the drag coefficient temporarily drops below the quasi-steady value. This observation is made for the second half of the period of the gust-plate interaction and can hardly be explained by the physics of the event. The additional source term and the evaluation of the shear stress tensor based on the total velocity are the two main differences between FVM and SVM.

As in the previous section, the altered SVM simulation, where the shear stress tensor is evaluated based on the background velocity (as for FVM) instead of the total velocity is applied to clarify the source of the changes. The temporal evolution of the force coefficients predicted by SVM and the altered SVM are quasi identical. This means that the SVM source term is responsible for the change in the $C_D^r$ shape. As explained in the previous section, this outcome is expected as the SVM source term contribution is more pronounced.

Please note that the evaluation of the shear stress tensor based on the total velocity does not have a major impact in the present case owing to the currently selected gust configuration. Indeed, the gust shape only changes in g₁-direction. If the gust would also vary in g₂- and/or g₃-direction, a shear stress tensor based on the total velocity would have a stronger effect on the results.

## 7. Conclusions

The field velocity method and the split velocity method are wind gust injection methods which are often used in the literature. Nevertheless, a comprehensive and comprehensible derivation in the context of a finite-volume scheme applying the Arbitrary Lagrangian Eulerian

formulation and the geometric conservation law was missing. Such a consistent derivation is provided here. This works out the similarities between the two approaches but also their differences. Their effects on the results are assessed based on vertical and horizontal gusts hitting a plate. The following results can be recorded:

* The derivation has highlighted that both methods are strongly linked to the Arbitrary Lagrangian Eulerian formulation and the geometric conservation law. The difference between both approaches are on the one hand the additional source terms in SVM describing the feedback effect of the surrounding flow field on the gust itself and the formulation of the viscous stress tensor by the background and the total velocity, respectively.

* From the numerical point of view and the required computational effort both methods are more or less equivalent. Therefore, using FVM does not really make sense anymore. This statement is supported by the computational results achieved by both methods.

* For the classical case of a vertical gust often used in the context of aircraft, the deviations in the results (flow field and force coefficients) are rather small, i.e., the normalized root-mean-square errors are in the range between 2.4% and 4.6%. An analytical derivation revealed that the source term vanishes in the regions where the dimensionless background velocity $\bar{u}_1$ is equal to unity (e.g., the far-field region), which explains its minor impact on the predictions in the vertical gust case.

* For the newly introduced horizontal gust hitting a plate the observed differences both in the flow field and in the resulting force coefficients are significantly more pronounced. The non-vanishing source term as well as the modified formulation of the viscous shear stresses can be made responsible for this observation. However, an additional study clearly showed that for the chosen Reynolds number the consideration of the feedback effect due to the source term is the decisive difference between FVM and SVM, which prevents unphysical flow structures from being observed on the plate.

The significance of the outcome of this study lies on the assessment of the capability of the prescribed velocity methods to introduce wind gusts of different configurations. Especially the horizontal gust case reveals the shortcomings of the underlying assumptions used in FVM and
<figure><img src="image_19.png" /><figcaption>Fig. 16. Horizontal gust passing around the rounded flat plate: Time history of the drag coefficient $C_D$, of the pressure part of the drag coefficient $C_D^p$ and of the shear part of the drag coefficient $C_D^r$ for FVM and SVM. The gust has the length of $L_g^i = 2c$.</figcaption></figure>

suggests the use of SVM for the injection of gusts having a streamwise velocity component. That will allow to tackle practically relevant flow problems in the near future, e.g., to predict the impact of strong wind gusts on flexible membranous structures and to answer the question of their structural integrity.

### CRediT authorship contribution statement

**Khaled Boulbrachene:** Conceptualization, Methodology, Software, Validation, Investigation, Visualization, Writing – original draft, Writing – review & editing. **Guillaume De Nayer:** Conceptualization, Methodology, Software, Validation, Investigation, Visualization, Writing – original draft. **Michael Breuer:** Conceptualization, Methodology, Resources, Writing – original draft, Writing – review & editing, Supervision, Project administration, Funding acquisition.

### 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 is financially supported by the *German Research Foundation (Deutsche Forschungsgemeinschaft)* under the contract number BR 1847/17-1.

## Appendix A. Evaluation and effect of the svm source term

### A.1. Vertical gusts

For a vertical gust, the gust velocity points only in the y-direction and is given by $\mathbf{u}_g = [0, u_{g,2}(x, t), 0]$. Hence, the source terms of the momentum equations simplify to $\mathbf{S} = [0, S_2(u_{g,2}(x, t)), 0]$. Employing Eq. (21), the simplified source term of the y-momentum equation reads:

$$S_2 = \rho \left\{ \frac{\partial u_{g,2}}{\partial t} + \tilde{u}_1 \frac{\partial u_{g,2}}{\partial x_1} \right\}. \quad (A.1)$$

The vertical gust component can be written in the Cartesian basis as:

$$u_{g,2}(x, t) = \begin{cases} \frac{v_{g,max}}{2} \left( 1 + \cos \left( \frac{2\pi (x_1 - x_{g,1}(t))}{L_g^{x_1}} \right) \right) & \text{for } (x_1 - x_{g,1}) \in \left[ -\frac{L_g^{x_1}}{2}, \frac{L_g^{x_1}}{2} \right] \\ 0 & \text{else,} \end{cases} \quad (A.2)$$

with the convected gust center defined as:

$$x_g(t) = x_g^0 + (t - t_0) u_\infty. \quad (A.3)$$

The derivatives in Eq. (A.2) are then:

$$\begin{aligned} \frac{\partial u_{g,2}}{\partial t} &= + \frac{\pi v_{g,max} u_\infty}{L_g^{x_1}} \sin \left( \frac{2\pi (x_1 - x_{g,1}(t))}{L_g^x} \right), \\ \frac{\partial u_{g,2}}{\partial x_1} &= - \frac{\pi v_{g,max}}{L_g^{x_1}} \sin \left( \frac{2\pi (x_1 - x_{g,1}(t))}{L_g^{x_1}} \right). \end{aligned} \quad (A.4)$$

The derivatives are therefore nearly identical but of opposite signs. Consequently, the source term (Eq. (A.1)) would vanish if the background velocity $\tilde{u}_1$ is equal to $u_\infty$. In fact, this is the case in the far-field region (far away from the plate), where the background velocity $\tilde{u}_1$ is equal to the free-stream velocity $u_\infty$. In that region, the source term plays no role and has no effect on the pressure field in contrast to the effects remarked in the horizontal gust configuration. As a result, no major deviations between FVM and SVM are realized for the vertical gust configuration (see Figs. 6 to 8).

### A.2. Horizontal gusts

For a horizontal gust, the gust velocity points only in the x-direction and is given by $\mathbf{u}_g = [u_{g,1}(x, t), 0, 0]$. Hence, the source terms of the momentum equations simplify to $\mathbf{S} =$

$$[S_1(u_{g,1}(x, t)), 0, 0]$. Employing Eq. (21), the simplified source term of the x-momentum equation reads:

$$S_1 = \rho \left\{ \frac{\partial u_{g,1}}{\partial t} + (\tilde{u}_1 + u_{g,1}) \frac{\partial u_{g,1}}{\partial x_1} \right\}. \quad (A.5)$$

The horizontal gust component can be written in the Cartesian basis as:

$$u_{g,1}(x, t) = \begin{cases} \frac{u_{g,max}}{2} \left( 1 + \cos \left( \frac{2\pi (x_1 - x_{g,1}(t))}{L_g^{x_1}} \right) \right) & \text{for } (x_1 - x_{g,1}) \in \left[ -\frac{L_g^{x_1}}{2}, \frac{L_g^{x_1}}{2} \right] \\ 0 & \text{else,} \end{cases} \quad (A.6)$$

with the convected gust center defined as:

$$x_g(t) = x_g^0 + (t - t_0) u_\infty. \quad (A.7)$$

The derivatives in Eq. (A.2) are then:

$$\begin{aligned} \frac{\partial u_{g,1}}{\partial t} &= + \frac{\pi u_{g,max} u_\infty}{L_g^x} \sin \left( \frac{2\pi (x_1 - x_{g,1}(t))}{L_g^{x_1}} \right), \\ \frac{\partial u_{g,1}}{\partial x_1} &= - \frac{\pi u_{g,max}}{L_g^x} \sin \left( \frac{2\pi (x_1 - x_{g,1}(t))}{L_g^{x_1}} \right). \end{aligned} \quad (A.8)$$
<table><caption>**Table B.6**
Fine and coarsened grids information and the relative error of the maximum lift coefficient.</caption><thead><tr><th>Gust wavelength</th><th>Grid</th><th>Total no. of CVs</th><th>C<sub>L,max</sub></th><th>Relative error %</th></tr></thead><tbody><tr><td rowspan="3">2 c</td><td>Fine</td><td>42,883,200</td><td>2.6115</td><td>–</td></tr><tr><td>xy-coarsened</td><td>10,720,800</td><td>2.4910</td><td>4.6</td></tr><tr><td>xyz-coarsened</td><td>5,360,400</td><td>2.4962</td><td>4.4</td></tr><tr><td rowspan="3">4 c</td><td>Fine</td><td>42,883,200</td><td>3.0178</td><td>–</td></tr><tr><td>xy-coarsened</td><td>10,720,800</td><td>3.1153</td><td>3.2</td></tr><tr><td>xyz-coarsened</td><td>5,360,400</td><td>3.1029</td><td>2.8</td></tr></tbody></table>

<figure><img src="image_20.png" /><figcaption>**Fig. A.17.** Vertical gust passing around the rounded flat plate: Instantaneous lift coefficients recorded on different grids for a grid sensitivity analysis. The gust shape is ECG with a gust amplitude of $A_g = 0.8$.</figcaption></figure>

The derivatives are again nearly identical and of opposite signs as in Appendix A.1. However, this time the spatial derivative is multiplied by the total velocity $u_1 = \bar{u}_1 + u_{g,1}$ which is different from $u_\infty$ even in the far-field region. Consequently, one could clearly observe the role that the source term (Eq. (A.5)) plays in the early stages of the gust injection by its strong effect on the pressure field (see Fig. 12). As a result, major deviations between FVM and SVM are clearly visible for the horizontal gust configuration.

## Appendix B. Grid sensitivity analysis

In order to examine the dependency of the results obtained on the grid used for the simulations, a grid sensitivity study is carried out. Two coarser grids are generated from the fine grid described in Section 5.1. Coarsening along a specific coordinate direction is performed by skipping every second grid line. The first grid is created by coarsening the blocks of the grid in the x- and y-directions only. This results in a computational grid containing 10,720,800 CVs. The second grid is generated by coarsening the blocks in all three coordinate directions resulting in a grid of 5,360,400 CVs. SVM simulations for a vertical gust configuration of different wavelengths (2c and 4c) are performed on the coarsened grids. The resulting lift coefficients obtained are plotted against those corresponding to the fine grid as a function of the convective time $t^*$ shown in Fig. A.17. It could be observed that no significant variations among the lift coefficients are visible and the curves are, almost all times, coincident. The maximum lift coefficient recorded is chosen as a parameter for quantifying the error between the coarsened grids and the fine grid (see Table B.6). Although the difference between

the total number of control volumes of the grids is significantly large (compared to the fine grid, the *xy*-coarsened grid contains 25% of the total number of control volumes and the *xyz*-coarsened grid has only 12.5%), the relative error of the lift coefficient does not exceed 5%. This indicates that the flow fields are overall well resolved in the vicinity of the plate even upon the coarsening. Hence, one could conclude that the use of the fine grid is appropriate for the purpose of this study.

Note that the oscillations observed in the time histories of the lift coefficient can be physically explained. Looking at Figs. 6 to 8 the cause of these fluctuations becomes obvious. When the vertical gust impacts the flat plate, vortical structures are shed especially from the leading edge but also from the trailing edge. As visible in the pressure distributions the center of these vortices have low pressures influencing the local pressure distribution on the plate and finally the resulting lift coefficients.

## References

* Apostolatos, A., De Nayer, G., Bletzinger, K.-U., Breuer, M., Wüchner, R., 2019. Systematic evaluation of the interface description for fluid-structure interaction simulations using the isogeometric mortar-based mapping. J. Fluids Struct. 86, 368–399.

* Badrya, C., Baeder, J.D., 2019. Numerical study of a flat plate wing response to large transverse gusts at low Reynolds number. AIAA Scitech 2019 Forum, San Diego, California, p. 0638, URL https://arc.aiaa.org/doi/abs/10.2514/6.2019-0638.

* Badrya, C., Biler, H., Jones, A.R., Baeder, J.D., 2021. Effect of gust width on flat-plate response in large transverse gust. AIAA J. 59 (1), 49–64.

* Biler, H., Badrya, C., Jones, A.R., 2019. Experimental and computational investigation of transverse gust encounters. AIAA J. 57 (11), 4608–4622.

* Breuer, M., 2002. Direkte Numerische Simulation und Large-Eddy Simulation turbulenter Strömungen auf Hochleistungsrechnern. In: Habilitationsschrift, Universität Erlangen–Nürnberg, Berichte aus der Strömungstechnik, Shaker Verlag, Aachen, Germany.

* 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. Fluids 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., 2018. Numerical studies on the instantaneous fluid-structure interaction of an air-Inflated flexible membrane in turbulent flow. J. Fluids 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., 2020. A source-term formulation for injecting wind gusts in CFD simulations. J. Wind Eng. Ind. Aerodyn. 207, 104405.

* 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. 25, Springer Nature Switzerland AG 2019, pp. 453–459.

* 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. & Fluids 99, 18–43.

* Demirdžić, I., Perić, M., 1988. Space conservation law in finite-volume calculations of fluid flow. Int. J. Heat Fluid Flow 8 (9), 1037–1050.

* Demirdžić, I., Perić, M., 1990. Finite-volume method for prediction of fluid flow in arbitrarily shaped domains with moving boundaries. Int. J. Numer. Methods Fluids 10 (7), 771–790.

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

* Ferziger, J.H., Perić, M., 2002. Computational Methods for Fluid Dynamics, third ed. Springer Berlin.
Franke, J., Baklanov, A., 2007. Best practice guideline for the CFD simulation of flows in the urban environment. In: COST Action 732 Quality Assurance and Improvement of Microscale Meteorological Models. University of Hamburg, Germany.

Franke, J., Hellsten, A., Schlünzen, K.H., Carissimo, B., 2011. The COST 732 best practice guideline for CFD simulation of flows in the urban environment: A summary. Int. J. Environ. Pollut. 44 (1-4), 419-427.

Frost, W., Long, B.H., Turner, R.E., 1978. Engineering handbook on the atmospheric environmental guidelines for use in wind turbine generator development. Tennessee Univ., Tullahoma (USA). Space Inst.; National Aeronautics and Space Administration, Huntsville, AL (USA). George C. Marshall Space Flight Center.

Ghoreyshi, M., Jirasek, A., Greisz, I., Satchell, M., 2018. Indicial aerodynamic modeling for arbitrary gust responses, in: Applied Aerodynamics Conference, June 25-29, 2018, Atlanta, Georgia, AIAA-2018-3965.

Gordnier, R.E., Visbal, M.R., 2015. Impact of a vortical gust on the aerodynamics of a finite aspect-ratio wing. In: 45th AIAA Fluid Dynamics Conference. pp. 1-18, URL https://arc.aiaa.org/doi/abs/10.2514/6.2015-3072.

Granlund, K., Monnier, B., Ol, M., Williams, D., 2014. Airfoil longitudinal gust response in separated vs. attached flows. Phys. Fluids 26 (2), 027103.

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, May 2, 2009.

Gu, N., Lu, Z., Guo, T., 2012. Simulation of viscous flows around a moving airfoil by field velocity method with viscous flux correction. Adv. Appl. Math. Mech. 4 (3), 294-310.

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: Contributions to the 18th STAB/DGLR Symposium, Stuttgart, Germany, 2012. Springer International Publishing, Cham, pp. 503-511.

Heinrich, R., Reimer, L., 2013. Comparison of different approaches for gust modeling in the CFD code TAU. In: International Forum on Aeroelasticity & Structural Dynamics 2013. URL https://elib.dlr.de/85834/.

Huntley, S., Jones, D., Gaitonde, A., 2016. 2D and 3D gust response using a prescribed velocity method in viscous flows, in: 46th AIAA Fluid Dynamics Conference.

Huntley, S.J., Jones, D., Gaitonde, A., 2017. Aeroelastic gust response of an aircraft using a prescribed velocity method in viscous flows. In: 23rd AIAA Computational Fluid Dynamics Conference. Denver, Colorado, URL https://arc.aiaa.org/doi/abs/10.2514/6.2017-3616.

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

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

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

Krigge, 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. Aerodyn. 153, 46-59.

Küssner, H.G., 1936. Zusammenfassender Bericht über den instationären Auftrieb von Flügeln. Luftfahrtforschung 13 (12), 410-424.

Kwon, D.K., Kareem, A., Butler, K., 2012. Gust-front loading effects on wind turbine tower systems. J. Wind Eng. Ind. Aerodyn. 104-106, 109-115.

Lesoinne, M., Farhat, C., 1996. Geometric conservation laws for flow problems with moving boundaries and deformable meshes, and their impact on aerelastic computations. Comput. Methods Appl. Mech. Engrg. 134 (1-2), 71-90.

Montenegro, P.A., Heleno, R., Carvalho, H., Calcada, R., Baker, C.J., 2020. A comparative study on the running safety of trains subjected to crosswinds simulated with different wind models. J. Wind Eng. Ind. Aerodyn. 207, 104398.

Münsch, M., 2015. Entwicklung und Anwendung eines semi-impliziten Kopplungsverfahrens zur numerischen Berechnung der Fluid-Struktur-Wechselwirkung in turbulenter Strömungen mittels Large-Eddy Simulationen. (Ph.D. thesis). University of Erlangen-Nuremberg, Germany.

Norris, S.E., Cater, J.E., Stol, K.A., Unsworth, C.P., 2010. Wind turbine wake modelling using large-eddy simulation. In: 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. Aerodyn. 162, 1-11.

Parameswaran, V., Baeder, J.D., 1997. Indicial aerodynamics in compressible flow—Direct computational fluid dynamic calculations. J. Aircr. 34 (1), 131-133.

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.

Sankar, L.N., Rou, S.Y., Malone, J.B., 1986. Application of surface transpiration in computational aerodynamics. In: 24th Aerospace Sciences Meeting. Reno, NV, USA, URL https://arc.aiaa.org/doi/abs/10.2514/6.1986-511.

Singh, R., Baeder, J.D., 1996. On the significance of transonic effects on aerodynamics and acoustics of blade-vortex interactions. In: Aeroacoustics Conference. Pennsylvania, USA, URL https://arc.aiaa.org/doi/abs/10.2514/6.1996-1697.

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

Sitaraman, J., 2003. CFD Based unsteady aerodynamic modeling for rotor aeroelastic analysis. (Ph.D. thesis). University of Maryland, College Park, USA.

Sitaraman, J., Baeder, J.D., 2000. Enhanced unsteady airload models using CFD. In: Fluids 2000 Conference and Exhibit. Denver, Colorado, URL https://arc.aiaa.org/doi/abs/10.2514/6.2000-2465.

Sitaraman, J., Baeder, J.D., 2006. Field velocity approach and geometric conservation law for unsteady flow simulations. AIAA J. 44 (9), 2084-2094.

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

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

Wu, Z., Cao, Y., Ismail, M., 2019. Gust loads on aircraft. Aeronaut. J. 123 (1266), 1216-1274.

Zaide, A., Raveh, D.E., 2006. Numerical simulation and reduced-order modeling of airfoil gust response. AIAA J. 44 (8), 1826-1834.