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

HAL
open science

# Combination of two FSI methods and their validation based on artificial wind gusts impacting a flexible T-structure

Michael Breuer, Khaled Boulbrachene, Guillaume de Nayer

► To cite this version:

Michael Breuer, Khaled Boulbrachene, Guillaume de Nayer. Combination of two FSI methods and their validation based on artificial wind gusts impacting a flexible T-structure. Computers and Fluids, 2025, 294, pp.106621. 10.1016/j.compfluid.2025.106621. hal-05027280

HAL Id: hal-05027280

https://hal.science/hal-05027280v1

Submitted on 9 Apr 2025

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_2.png" /></figure>

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

# Combination of two FSI methods and their validation based on artificial wind gusts impacting a flexible T-structure

M. Breuer<sup>DID</sup> *, K. Boulbrachene<sup>DID</sup>, G. De Nayer<sup>DID</sup>

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

## ARTICLE INFO

**Keywords:**

Artificial (horizontal) wind gust
Fluid-structure interaction (FSI)
Coupled numerical FSI simulations
Immersed boundary method
Wind tunnel experiments
Validation

## ABSTRACT

The study focuses on the combination of two numerical approaches that are typically not used together in this manner. The first is a well-established partitioned fluid-structure interaction (FSI) simulation methodology relying on a finite-volume fluid solver for curvilinear, block-structured, body-fitted grids written in the Arbitrary Lagrangian-Eulerian (ALE) formulation, and a finite-element solver for the structural analysis. The second approach is an immersed boundary (IB) method employing a continuous and direct forcing strategy. The IB method, often applied to Cartesian grids, is also referred to as an approach to simulate fluid-structure interactions. In this study, both methods are combined to exploit their respective advantages in simulating a complex flow problem. The coupled FSI problem involves the interaction of a thin, flexible structure deforming under the dynamic load of a wind gust (task 1). The gust itself is generated by an artificial wind gust generator, which includes a *paddle* that partially obstructs the wind tunnel's outlet, thereby defining an FSI problem of its own (task 2). For task 1, the classical partitioned ALE approach is employed, while the IB method is more appropriate for task 2. Using available experimental measurement data for both the fluid flow and the structural deformation, the combined simulation framework is first validated for the case without gust. In a second step, the more challenging FSI problem of discrete gusts impacting the T-structure is thoroughly analyzed and the predicted data are compared with the available measurement data. For both cases without and with gusts, a very good agreement between simulation and experiment is achieved, which justifies the chosen approach.

## 1. Introduction

Extreme weather phenomena characterized by high wind velocities and strong wind gusts have been observed more and more frequently in recent years. Structures exposed to such wind conditions are not only subjected to high static forces. Dynamic loads caused by the short-term wind gusts are also particularly relevant for stability and fatigue strength. This is especially true for constructions made of thin structures and/or lightweight materials, which are vulnerable to wind events.

To better understand the impact of strong wind gusts on such constructions, conducting measurements on real structure may be seen like the most straightforward approach for investigating fluid-structure interactions. However, simultaneously measuring the fluid flow and the structure deformation at real scales is extremely challenging and costly. Additionally, the difficulty exists that the wind speed and gust strength cannot be controlled in nature. Instead, one must set up the measuring equipment and wait for the right weather conditions to occur. To conduct these investigations in a laboratory setting, it is necessary to develop methods to artificially generate wind gusts in a wind tunnel.

Recently, Wood et al. [1] conducted a comprehensive literature review on methods for generating artificial wind gusts in a laboratory environment. Generally speaking, wind gusts can be divided into two classes, vertical and horizontal gusts. Vertical gusts are directed perpendicular to the main flow direction and are particularly relevant in aeronautics leading to rapid changes in lift forces when an airplane is subjected to such a gust, c.f., [2–5]. Contrarily, in the horizontal case the gust is directed in the main flow direction. This gust type is especially of relevance in civil engineering applications, i.e., for the stability of bridges, high-rise buildings or thin-walled constructions. Moreover, energy harvesting wind turbines are also vulnerable to such gusts [6] inducing a variety of studies, c.f., [7–10].

Building on their literature review and concentrating on horizontal wind gusts, Wood et al. [1] developed a novel artificial wind gust generator (WGG) known as the *paddle*. By partially obstructing the outlet of the wind tunnel nozzle with the *paddle* within a short time interval, the flow is significantly accelerated leading to a horizontal gust. In Wood et al. [1] and Wood and Breuer [11] it is demonstrated that the shape and the intensity of the generated wind gust can be

<footnote>* Corresponding author.</footnote>

<footnote>E-mail address: breuer@hsu-hh.de (M. Breuer).</footnote>
effectively controlled achieving high reproducibility and allowing the superposition of these wind gusts onto the free-stream flow.

The possibility to synchronize the wind gust generator with the measurement equipment for both the fluid flow and the structural deformation using a trigger signal was first used in a follow-up study by Breuer and Neumann [12]. The aim of this study was to develop a measurement strategy capable of simultaneously capturing the fluid flow and the structural deformation during the impact of the gust with high temporal resolution. For this purpose, a synchronized setup was utilized, comprising the previously described WGG and two distinct high-speed optical measurement techniques, i.e., particle-image velocimetry (PIV) for the fluid flow analysis and digital-image correlation (DIC) for monitoring the structural deformations. More than 100 gust events with consistent characteristics were recorded, enabling the comparison of different gust realizations and facilitating phase averaging. Designed as a benchmark case, a geometrically simple T-shaped thin structure was selected as the test object. The structure was mounted onto the bottom wall of the test section and subjected to the fluid load of the impacting gusts. The data generated in this study are employed in this paper for comparative analysis.

In addition to experimental investigations, engineers are also interested in describing such coupled FSI problems using advanced numerical simulations. This naturally raises the question of how to accurately reproduce the experimental generation of wind gusts observed in the wind tunnel within a numerical framework. The study by Boulbrachene and Breuer [13] addressed this question in detail. To describe the motion of the paddle that partially blocks the nozzle's outlet and is responsible for the gust generation [1], the immersed boundary (IB) method with a continuous and direct forcing approach was implemented into a finite-volume flow solver for large-eddy simulations. The IB method ensures smooth pressure and velocity fields without any numerically induced oscillations when the paddle moves through the fluid stream. Compared to alternative approaches for generating gusts, which will be discussed in more detail in Section 2.2, the IB method offers the most realistic option for reproducing the gusts generated by the experimental device ('numerical wind tunnel'). The simulation results in Boulbrachene and Breuer [13] demonstrated excellent agreement with the measurements conducted by Wood et al. [1] at specific points in the flow field, validating the robustness and accuracy of the proposed simulation methodology.

The objective of this study is to create a numerical counterpart of the experimental investigations on the T-structure under wind gust load conducted by Breuer and Neumann [12]. To achieve this goal, the IB gust generation methodology described above is combined with the possibilities of coupled high-fidelity simulations for flexible structures in complex turbulent flows. This combination is a novelty. In contrast to the experiment, all field quantities, including velocities, pressure and deformation, will be provided in the entire domain with high spatial and temporal resolution. This measure shall facilitate a deeper analysis of the complex fluid-structure interaction of the thin structure and the gust-induced flow. The T-shaped thin structure mounted in the wind tunnel can be considered as a first step to prove and validate the newly combined methodology. As an initial step, the test case offers several advantages: the structure is geometrically simple, it consists of a thin brass sheet with easily determinable material properties, and, finally, comprehensive experimental measurements are available for validation and comparison.

The paper is organized as follows: Section 2 outlines the entire numerical methodology used for simulating both the fluid flow and the structural deformation, including the coupling procedure between the two. Section 3 elaborates on the definition of the test case, covering the physical properties as well as the numerical setup employed in the simulations. Section 4 validates and analyzes the results obtained in two main substeps: As a first step, Section 4.1 serves as an additional validation of the already established ALE FSI computational framework. The IB method is then activated within the ALE procedure and the entire framework is validated against the experimental results in Section 4.2.

## 2. Numerical methods

Numerically simulating the coupled FSI problem, which involves the turbulent flow including wind gusts and the flexible structure, can be broken down into four key subtasks: Calculating the turbulent flow field, accurately representing the generation of gusts in the wind tunnel, determining the structural deformation, and finally coupling the fluid and structural simulations. For each individual subtask, a general description, the relevant equations, and previous publications thoroughly addressing these topics are provided below. The interactions among these subtasks are illustrated by the flowchart depicted in Fig. 1.

### 2.1. Numerical description of the fluid flow

The complex turbulent flow is predicted using the large-eddy simulation (LES) technique. The underlying computational fluid dynamics (CFD) solver, FASTEST-3D [14], solves the filtered Navier-Stokes equations in an Arbitrary Lagrangian-Eulerian (ALE) formulation [15]. These equations are discretized using the finite-volume technique on a curvilinear, block-structured, body-fitted grid with a collocated variable arrangement. The integral form of the governing equations formulated for a temporally varying domain read:

$$ \frac{d}{dt} \int_{V(t)} \rho_f \mathbf{u} \, dV + \int_{S(t)} \rho_f (\mathbf{u} - \mathbf{u}^{\text{grid}}) \mathbf{u} \cdot \mathbf{n} \, dS = - \int_{S(t)} \tau \cdot \mathbf{n} \, dS \\ - \int_{S(t)} \rho \mathbf{n} \, dS, \quad (2) $$

where $\mathbf{u}$ is the flow velocity, $\rho_f$ is the density of the fluid, $\tau$ is the stress tensor, and $\mathbf{n}$ is the normal unit vector of the surface $S(t)$ belonging to the control volume $V(t)$. Since the grid is deformable, $\mathbf{u}^{\text{grid}}$ denotes the grid velocity with which the surface of a control volume is moving.

Surface and volume integrals appearing in Eqs. (1) and (2) are approximated using the midpoint rule. The pressure and the viscosity are linearly interpolated to the cell faces, resulting in a second-order accurate central scheme. To avoid wiggles in the simulation data, the velocities used for the prediction of the convective fluxes on the cell faces are approximated using flux blending [16], which combines 3% of a first-order accurate upwind scheme with 97% of a second-order accurate central scheme. Furthermore, the momentum interpolation technique of Rhie and Chow [17] is applied to predict the mass fluxes through the cell faces preventing unwanted oscillations on non-staggered grids. To couple the pressure and velocity fields, a semi-implicit predictor-corrector scheme is employed for solving the incompressible Navier-Stokes equations [15]. First, the momentum equations are time-marched using a low-storage multi-stage Runge-Kutta method to derive an intermediate velocity $\mathbf{u}^*$. Then, the corrector step ensures mass conservation, yielding a divergence-free velocity field by solving a Poisson equation for the pressure correction. In the LES technique, the large scales of turbulence are resolved while the small scales are modeled with a subgrid-scale (SGS) model [18]. This study uses the classical Smagorinsky [19] model with a standard parameter $C_s = 0.1$ and the Van-Driest damping function for the near-wall region. Given that the near-wall grid resolution is sufficient to resolve the viscous sublayer, as detailed in Section 3, wall functions are not required, enabling a fully wall-resolved LES approach.

In the FSI simulation framework, the structural deformations resulting from fluid forces $f_{\text{ALE}}$, which are computed based on the fluid pressure and shear stresses acting on the structure, are taken into account at the interface between the fluid and the solid. This necessitates that the ALE formulation accounts for the fluid motion within a temporally varying domain, requiring a body-fitted grid that must be adapted at each time step. Several techniques are available for grid adaptation depending on the magnitude of the structural deformations.
<figure><img src="image_4.png" /><figcaption>Fig. 1. Flowchart of the computational framework.</figcaption></figure>

For small deformations as in the present test case, a very fast algebraic procedure is employed, utilizing a combination of linear and transfinite interpolations (TFI) [20]. This method is efficient, with its computational effort scaling proportionally to the number of volume grid points, and is well-suited for scenarios where the structural displacements do not significantly alter the fluid domain. A measure for the grid quality was introduced by Knupp [21] as the skew quality metric. For LES predictions, this value should not fall below 0.9 during the grid adaptation to guarantee a high grid quality. For alternatives in case of substantial deformations, refer to Sen et al. [22].

## 2.2. Numerical description of the gust generation

### 2.2.1. Overview

In the literature, continuous and discrete gust models are distinguished [23,24]. Continuous gust models, which describe wind conditions in form of multi-correlated wind-velocity time series, are often applied in the context of fatigue evaluation. However, for the determination of extreme loads they are rather rarely employed. For this purpose, discrete models, describing the gust in a simplistic manner via mathematical functions representing the spatial and temporal velocity distributions [6], are more appropriate. Therefore, the brief overview is limited to this category.

A discrete wind gust can be described within a simulation environment using various approaches. One classical method for introducing wind gusts into the computational domain is to superimpose them at the inlet boundary, allowing the gusts to freely propagate through the domain. This technique is commonly known as the far-field boundary condition method. To maintain the global mass balance, especially in simulations involving incompressible fluid flows, this method requires dynamic adjustments of the inlet velocity in regions away from the gust. This measure ensures a constant total mass flow through the inlet [25,26]. Additionally, a highly resolved grid is essential to prevent significant numerical damping of the gust while being convected through the flow domain toward the region of interest.

To address the need for fine grids throughout the computational domain, several alternative approaches have been proposed. One such approach is the pure source-term formulation method, as suggested by De Nayer and Breuer [27]. In this method, the wind gust is injected at any position in the computational domain by adding a source term to the momentum conservation equation. This approach offers the advantage of maintaining the global mass conservation while also accounting for the feedback effect of the flow field on the gust itself. This method circumvents the need for globally fine grids, thereby reducing computational costs.

Another class of approaches, often used in aeronautical applications for evaluating individual responses of airplanes, is referred to as the prescribed velocity method. Unlike the previously mentioned methods, which allow gusts to propagate freely after the injection phase, the prescribed velocity methods define the position of the imposed field at each instant in time. Within this category, two distinct approaches are recognized: the field velocity method (FVM) and the split velocity method (SVM). The field velocity method was initially proposed by Parameswaran and Baeder [28] and Singh and Baeder [29]. The key idea behind FVM is to modify the convective fluxes of the mass and momentum conservation equations by superimposing a pre-defined additional velocity that varies in space and time according to a gust profile. To account for the missing feedback effect of the surrounding flow on the gust itself, the split velocity method was introduced, as originally proposed by Wales et al. [4]. In SVM, the velocity field is decomposed into a prescribed gust velocity and the remaining background velocity, as denoted by Huntley et al. [30]. The momentum equations are effectively rearranged to solve for the background velocity, with additional source terms included to account for the feedback effect. Typical applications of this method are found in the works of Biler et al. [31], Badrya and Baeder [32], Badrya et al. [5] and Boulbrachene et al. [33].

All discrete gust approaches described so far rely on the assumption that the velocity distribution within the gust is both known and can be explicitly prescribed. Typically, deterministic functions are employed
for the spatial and temporal distributions, such as the Extreme Coherent Gust (ECG) or the Extreme Operating Gust (EOG) as defined by IEC-Standard [6]. However, in the current study, where the gust is generated within the wind tunnel using the WGG, these velocity distributions are not explicitly known. Although the experimental investigations conducted by Wood et al. [1] aimed to replicate the 1-cosine gust shape characteristic of ECG gusts, this does not mean that the artificial gust generated by the paddle precisely matches the corresponding velocity profile. To realistically reproduce the gust generated in the wind tunnel, the generation process of the gust involving the moving paddle must be taken into account. Unfortunately, none of the previously mentioned methods is capable of accurately capturing this process. Consequently, an alternative approach, as suggested by Boulbrachene and Breuer [13], must be applied to achieve a more realistic simulation. This alternative approach will be explained next.

### 2.2.2. Immersed boundary gust generator

The working principle of the wind gust generator installed in a Göttingen-type wind tunnel with an open test section is straightforward (see Fig. 2): Initially, the nozzle of the wind tunnel is fully open. The fluid flows through the measurement section between the nozzle outlet and the receiver. The gust is generated by partially blocking the outlet of the nozzle for a short time period. This closure operation is achieved by a rigid plate (referred to as the paddle) mounted on a movable saddle of a tooth belt axis. All key kinematic parameters of the tooth belt axis – such as stroke length, velocity and acceleration – can be fully controlled by the software allowing to generate various gust shapes, as detailed in Wood et al. [1].

The finite-volume method described in Section 2.1 in the Eulerian (E) frame of reference forms the basis for calculating the flow in the cell centers $x_{i,j,k}$ of the CFD grid $g_h$. The paddle is modeled as a moving boundary employing the immersed boundary (IB) method with a direct forcing approach, as proposed by Uhlmann [34]. The effect of the immersed boundary is incorporated into the governing equations in the ALE formulation by adding a forcing term to the right-hand side of the momentum equation:

$$ \frac{d}{dt} \left( \int_{V(i)} \rho_f \mathbf{u} \, dV + \int_{S(i)} \rho_f (\mathbf{u} - \mathbf{u}_{\text{grid}}) \mathbf{u} \cdot \mathbf{n} \, dS \right) = - \int_{S(i)} \bar{\tau} \cdot \mathbf{n} \, dS - \int_{S(i)} \rho_f \mathbf{n} \, dS + \underbrace{\int_{V(i)} \mathbf{f}^{\text{IB}} \, dV}_{\text{IB forcing term}} \quad \text{N.-S. in ALE formulation} \quad (3) $$

This volume force $\mathbf{f}^{\text{IB}}$, applied at the cell centers of the Eulerian grid, represents the influence of the solid paddle on the surrounding fluid. The immersed boundary $\Gamma^{\text{IB}}$ is represented by a set of Lagrangian (L) markers with coordinates denoted as $\mathbf{X}_i$ (i.e., $\mathbf{X}_i \in \Gamma \ \forall 1 \le i \le N_L$, where $N_L$ represents the total number of Lagrangian markers). In the direct forcing scheme [34], the boundary force is first evaluated at the Lagrangian force points as:

$$ \mathbf{F}_i^{\text{IB}} = \rho_f \frac{\mathbf{U}_i^d - \mathbf{U}_i^*}{\Delta t} \quad \forall \mathbf{X}_i, \quad (4) $$

where $\mathbf{U}_i^d$ is the desired velocity at the fluid-solid interface and $\mathbf{U}_i^*$ is the intermediate velocity after the predictor step. Due to the movement of the paddle and the resulting time-dependent immersed boundary, the fluid-solid interface represented by the Lagrangian markers is not aligned with the cell centers of the Eulerian fluid grid. To effectively couple the fluid and the solid phases, it is necessary to transfer quantities between the Eulerian and the Lagrangian frames of reference. Therefore, the original IB method introduced by Peskin [35] depending on the use of a regularized delta function $\delta_h$ [36] for the transfer of quantities is employed. The coupling is achieved through the interpolation operation $I(\mathbf{u}^*)$ for the velocity ($E \to L$) and the spreading operation $S(\mathbf{F}^{\text{IB}})$ for the forces ($L \to E$) defined as:

$$ \mathbf{u}^*(\mathbf{X}_i) = \sum_{i,j,k \in g_h} \mathbf{u}^*(\mathbf{x}_{i,j,k}) \delta_h(\mathbf{x}_{i,j,k} - \mathbf{X}_i) \Delta v_{i,j,k} \quad \forall 1 \le i \le N_L, \quad (5) $$

$$ \mathbf{f}^{\text{IB}}(\mathbf{x}_{i,j,k}) = \sum_{l=1}^{N_L} \mathbf{F}^{\text{IB}}(\mathbf{X}_l) \delta_h(\mathbf{x}_{i,j,k} - \mathbf{X}_l) \Delta V_l \quad \forall \mathbf{x}_{i,j,k} \in g_h. \quad (6) $$

Hence, the steps of the predictor-corrector method accounting for the immersed boundary can be summarized by the following equations:

$$ 1. \quad \mathbf{u}^* = \mathbf{u}^n + \frac{\Delta t}{\rho_f} \mathbf{r} \mathbf{h} \mathbf{s}^n \quad \forall \mathbf{x}_{i,j,k} \in g_h, \quad (7) $$

$$ 2. \quad \mathbf{U}^* = I(\mathbf{u}^*) \quad \forall \mathbf{X}_l \in \Gamma, \quad (8) $$

$$ 3. \quad \mathbf{F}^{\text{IB}} = \rho_f \frac{\mathbf{U}_i^d - \mathbf{U}_i^*}{\Delta t} \quad \forall \mathbf{X}_l \in \Gamma, \quad (9) $$

$$ 4. \quad \mathbf{f}^{\text{IB}} = S(\mathbf{F}^{\text{IB}}) \quad \forall \mathbf{x}_{i,j,k} \in g_h, \quad (10) $$

$$ 5. \quad \mathbf{u}^{**} = \mathbf{u}^* + \frac{\Delta t}{\rho_f} \mathbf{f}^{\text{IB}} \quad \forall \mathbf{x}_{i,j,k} \in g_h, \quad (11) $$

$$ 6. \quad \Delta p' = \frac{\rho_f}{\Delta t} \nabla \cdot \mathbf{u}^{**} \quad \forall \mathbf{x}_{i,j,k} \in g_h, \quad (12) $$

$$ 7. \quad \mathbf{u}^{n+1} = \mathbf{u}^n - \frac{\Delta t}{\rho_f} \nabla p' \quad \forall \mathbf{x}_{i,j,k} \in g_h. \quad (13) $$

The flow field around the T-structure in the test section simulated in this work necessitates a stretched grid that clusters grid points near the walls to accurately capture the boundary layers. As a result, a general curvilinear non-uniform grid is employed. To account for the non-uniformity of the grid, a modified window function is defined according to Pinelli et al. [37]. This window function includes 10 unknown polynomial coefficients, which are determined by imposing a set of discrete modified moment conditions to ensure the integral conservation of both force and moment. The number of polynomial coefficients arises from the three-dimensionality of the problem and the second-order accurate spatial discretization scheme used [37]. The resulting linear system of equations is solved using the linear algebra package LAPACK [38]. For further details, readers are referred to Boulbrachene and Breuer [13].

The tracking of Lagrangian markers is performed in the computational space (c-space), which is derived by transforming the structured grid connecting the cell centers to an orthonormal grid, as described by Schäfer and Breuer [39]. By using this transformation, the time-consuming search process in the physical space (p-space) associated with finding the new computational cell containing the Lagrangian marker is avoided [40]. In c-space, the location of the Lagrangian marker is defined as a combination of an index triple, representing the cell containing the marker and a fractional offset within that cell. This approach not only facilitates efficient searching of the nearest cell center, necessary for determining the bounds of the regularized delta function used in the IB method, but also allows for the conversion of the marker's position back to p-space. This transformation from c-space to p-space is accomplished using a trilinear interpolation based on the cell-centered coordinates [39]. In addition to its role in evaluating the precise location of the Lagrangian markers in physical coordinates, this interpolation is also crucial for visualizing the immersed boundary in p-space.

The immersed boundary (i.e., the paddle) is set into motion according to a prescribed velocity defined in the physical space. To accurately update the marker's location in c-space, it is first necessary to transform the prescribed velocity components from p-space to c-space. This transformation is accomplished using the spatial metric coefficients followed by a temporal integration. The metric coefficients are the entries of the Jacobian matrix of the p-space to c-space coordinate transformation and are evaluated at the Lagrangian location. For this purpose, the CZ+ tracking scheme proposed by Schäfer and Breuer [39] is employed. Further details can be found in Boulbrachene and Breuer [13].

### 2.3. Numerical description of the structural deformation

The deformation of the T-shaped structure due to fluid loads at the FSI interface is predicted using the Carat++ solver. This software specializes in the finite-element and isogeometric analysis of structures, with focus on predicting the geometrically non-linear behavior of thin-walled components, such as shells and membranes [41,42]. In the
current study, the structure is assumed to behave in a linear-elastic and isotropic manner and a St. Venant-Kirchhoff material law is employed to link the second Piola-Kirchhoff stress tensor S with the Green-Lagrange strain tensor E [43]. The dynamic equilibrium of the structure is described by the momentum equation (Cauchy's first equation of motion) written in a Lagrangian frame of reference and expressed with the second Piola-Kirchhoff stress tensor as:

$$ \rho_s \frac{\partial^2 \mathbf{d}}{\partial t^2} = \nabla \cdot (\mathbf{F} \cdot \mathbf{S}) + \rho_s \mathbf{b}. \quad (14) $$

Here, $\rho_s$ represents the density of the structure, $\mathbf{d}$ denotes the displacement of the material point in space, $\mathbf{F}$ is the material deformation gradient, and $\mathbf{b}$ represents the specific volume forces. To solve the problem in the finite-element context numerically, the weak form of Eq. (14) is formulated based on the principle of the virtual work. To derive the principle of virtual work, Eq. (14) is multiplied by an arbitrary test function (chosen here as the variation of the displacement field $\delta \mathbf{d}$) and integrated over the volume of the structural domain $\Omega_s$. In a state of equilibrium, the resulting virtual work $\delta W$ that is performed by the inertial, internal, and external forces of the system due to the virtual displacement field $\delta \mathbf{d}$ shall be zero. The weak form then reads:

$$ \delta W = \underbrace{\int_{\Omega_s} \rho_s \frac{\partial^2 \mathbf{d}}{\partial t^2} \cdot \delta \mathbf{d} \, d\Omega_s}_{\delta W_{\text{dyn}}} + \underbrace{\int_{\Omega_s} \delta \mathbf{E} : \mathbf{S} \, d\Omega_s}_{\delta W_{\text{int}}} - \underbrace{\left( \int_{\Gamma_s} \hat{\mathbf{T}} \cdot \delta \mathbf{d} \, d\Gamma_s + \int_{\Omega_s} \rho_s \mathbf{b}_s \cdot \delta \mathbf{d} \, d\Omega_s \right)}_{\delta W_{\text{ext}}} = 0, \quad (15) $$

where $\delta \mathbf{E}$ is the virtual strain and $\hat{\mathbf{T}} = \mathbf{n} \cdot (\mathbf{F} \cdot \mathbf{S})$ is the traction on the boundary $\Gamma_s$. Given the weak form, a spatial and a temporal discretization is conducted to construct a non-linear system of equations that can be solved in an iterative manner. The continuous displacement field $\mathbf{d}(\mathbf{x})$ is approximated by interpolating the discrete nodal displacements $\bar{\mathbf{d}}$ with the help of shape functions $N(\theta^1, \theta^2, \theta^3)$ (where $\theta$ is the curvilinear convective coordinate):

$$ \mathbf{d}(\theta^1, \theta^2, \theta^3) \approx \sum_{i=1}^{n_{\text{nod}}} N_i(\theta^1, \theta^2, \theta^3) \bar{\mathbf{d}}_i. \quad (16) $$

Here, $n_{\text{nod}}$ is the number of nodes within the finite element. Inserting Eq. (16) in Eq. (15) yields the semi-discrete form:

$$ \mathbf{M} \ddot{\mathbf{d}} + \mathbf{f}^{\text{int}}(\bar{\mathbf{d}}) = \mathbf{f}^{\text{ext}}, \quad (17) $$

where $\mathbf{M}$ is the mass matrix, $\mathbf{f}^{\text{int}}$ is the internal force (non-linear function of $\bar{\mathbf{d}}$), and $\mathbf{f}^{\text{ext}}$ is the external force. The damping effect in the structure is modeled by adding a damping matrix $\mathbf{C}$ in Eq. (17) as:

$$ \mathbf{M} \ddot{\mathbf{d}} + \mathbf{C} \dot{\mathbf{d}} + \mathbf{f}^{\text{int}}(\bar{\mathbf{d}}) = \mathbf{f}^{\text{ext}}. \quad (18) $$

Given the significance of internal structural damping in such scenarios, the well-established Rayleigh damping method is applied. The constant mass- and stiffness-proportional parameters ($\alpha_r, \beta_r$) of the Rayleigh damping are evaluated based on additional experimental structural tests, as detailed in Appendix A.

For time integration, the structural response is discretized using the standard second-order Newmark scheme. This method ensures accurate prediction of dynamic structural responses under fluid-induced loads.

## 2.4. Coupling procedure

The simulation framework for addressing the FSI problem is based on a partitioned two-way coupling procedure [15], as sketched in Fig. 1. This method integrates the fluid solver (FASTEST-3D, discussed in Section 2.1) with the computational structure dynamics (CSD) solver (Carat++, Section 2.3). The coupling between these two solvers, as well as the mapping of data between the non-matching CFD and CSD surface grids at the FSI interface, is managed by the open-source software

EMPIRE (Enhanced Multi Physics Interface Research Engine, [44]). The mapping process, critical for accurate data transfer between the CFD and CSD solvers, utilizes a mortar mapping technique suitable for the FEM discretizations. Comprehensive details on the mortar mapping approach can be found in the works of Wang et al. [45] and Apos-tolatos et al. [46]. The framework supports both standard loose and strong coupling schemes, ensuring flexibility in handling various FSI scenarios. All data exchange between the fluid and structural solvers relies on message passing interface (MPI) communications, which ensures efficient parallel processing and synchronization of the coupled simulations.

This FSI simulation framework without IBM has been validated using a variety of test cases, including experimental FSI benchmarks, i.e., De Nayer et al. [47], De Nayer et al. [48], and De Nayer and Breuer [49] to cite just the latest studies. These studies provide a robust basis for the reliability and accuracy of the FSI framework. However, it is important to note that the IB technique, integrated to the methodology in this study, was validated separately in Boulbrachene and Breuer [13] and Boulbrachene and Breuer [50]. What makes the current study distinctive is the integration of the FSI simulation framework with the IB technique – a combination that has not been previously applied. This novel integration is a special feature of this research, allowing for the accurate simulation of complex fluid-structure interactions where additional moving boundaries are involved, as in the case of the WGG discussed earlier. The successful application of this combined methodology (see Fig. 1) demonstrates its potential for advancing the simulation of dynamic FSI problems.

## 3. Test case configuration

The current numerical study is based on two aforementioned experimental investigations: the implementation and testing of the wind gust generator in the wind tunnel conducted by Wood et al. [1], and the synchronized measurements of the T-structure under wind gust load in the wind tunnel by Breuer and Neumann [12] (see Fig. 2). These studies serve as the experimental basis for validating the numerical simulations in this work. The relevant physical parameters<fnref n="1" /> will be described next.

### 3.1. Physical parameters

The subsonic wind tunnel with an open test section of the length $l$ has a rectangular nozzle outlet (width $w \times$ height $h$). The parameters relevant to the fluid flow and the gust generation are summarized in Table 1. The wind tunnel blower is operated at a low speed, resulting in a low free-stream velocity $u_\infty$. To generate a horizontal gust, the paddle of the WGG performs a rapid downward and upward movement, thus obstructing the nozzle outlet for a short period of time $\Delta t_{\text{gust}}$, denoted the gust generation time. When the paddle reaches its maximum displacement $\Delta y_{\text{max}}$, the maximum blocking ratio $BR$ of the nozzle is attained, which approximately coincides with the instant at which the maximum velocity of the artificially generated horizontal gust is achieved. Wood et al. [1] investigated five different motion patterns of the paddle and analyzed their effect on the gusts generated. The present study focuses on gust type 1, which features a fast and symmetric motion pattern designed to mimic the 1-cosine shape of the ECG gust as defined by IEC-Standard [6]. The characteristic parameters of the gust generated by this motion pattern are outlined in Table 1. It is important to note that the chosen gust case leads to a strong variation of the velocity at the outlet of the nozzle with a maximum gust velocity, which is about 63 % higher than the free-stream velocity [1]. For detailed

<footnote n="1">¹ Note that the experimental study by Breuer and Neumann [12] and the numerical study by Boulbrachene and Breuer [13] use different coordinate systems. In this study, we follow the latter, where the direction normal to the ground plate is defined as the y-axis.</footnote>
<figure><img src="image_5.png" /><figcaption>Fig. 2. Sketch of the wind tunnel with the open test section, the wind gust generator, and the T-structure mounted on the ground plate.</figcaption></figure>

<table><caption>**Table 1**
Fluid flow and gust generation parameters.</caption><tr><td>Free-stream velocity</td><td>u<sub>∞</sub></td><td>5.14</td><td>m/s</td></tr><tr><td>Density (ρ = 20 °C)</td><td>ρ<sub>f</sub> = ρ<sub>air</sub></td><td>1.225</td><td>kg/m<sup>3</sup></td></tr><tr><td>Dynamic viscosity</td><td>μ<sub>air</sub></td><td>18.27 × 10<sup>-6</sup></td><td>kg/(m s)</td></tr><tr><td>Width of the nozzle outlet</td><td>ω</td><td>0.5</td><td>m</td></tr><tr><td>Height of the nozzle outlet</td><td>h</td><td>0.375</td><td>m</td></tr><tr><td>Test section length in wind tunnel</td><td>l</td><td>0.8</td><td>m</td></tr><tr><td>Maximum displacement of paddle in flow</td><td>Δy<sub>max</sub></td><td>0.1437 (0.145)<sup>a</sup></td><td>m</td></tr><tr><td>Maximum blocking ratio</td><td>BR = Δy<sub>max</sub>/h</td><td>0.387</td><td>—</td></tr><tr><td>Duration of paddle in flow</td><td>Δt<sub>gust</sub></td><td>0.352 (0.356)<sup>a</sup></td><td>s</td></tr><tr><td>Extremal downward velocity</td><td>u<sup>d</sup><sub>WGG</sub></td><td>-1.05</td><td>m/s</td></tr><tr><td>Extremal upward velocity</td><td>u<sup>u</sup><sub>WGG</sub></td><td>+0.95</td><td>m/s</td></tr><tr><td>Extremal downward acceleration</td><td>Δd<sup>d</sup><sub>WGG</sub></td><td>-20</td><td>m/s<sup>2</sup></td></tr><tr><td>Extremal upward acceleration</td><td>Δu<sup>u</sup><sub>WGG</sub></td><td>+15</td><td>m/s<sup>2</sup></td></tr></table>

time histories of the paddle's kinematics – displacement, velocity and acceleration – please refer to Fig. 4 in Wood et al. [1].

Designed as a benchmark case, a geometrically simple T-shaped structure is mounted perpendicularly on the flat and smooth ground plate, which serves as the bottom boundary of the wind tunnel's test section, as sketched in Fig. 2. The structure consists of a thin sheet of brass (alloy CuZn37) with the thickness $t_s$. The dimensions of the rectangular thin rod ($b_1, h_1$) and the shield-shaped part at the top ($b_2, h_2$) are given in Table 2. The structure is mounted in the symmetry plane of the test section at a distance of $d_s$ from the nozzle's outlet. Gravitational acceleration with $|g| = 9.81 \text{ m/s}^2$ in negative y-direction is considered. The material properties such as the density $ρ_s$, the Young's modulus $E_s$ and the Poisson's ratio $v_s$ given in Breuer and Neumann [12] were estimations based on literature values. Thorough investigations on the material used in the experiment are detailed in Appendix A. These tests showed that the value of the Young's modulus assumed in Breuer and Neumann [12] was too high. The revised value is also provided in Table 2. It is worth noting that this initial false estimation of the Young's modulus has no influence on the measurement data [12] of the flow field and the structural deformation used for comparison purpose in the present study.

The Reynolds number for the undisturbed free-stream flow around the T-structure, calculated based on the entire height $h_s = h_1 + h_2$ and the properties provided in Table 1, is approximately $Re = 41,360$.

<table><caption>**Table 2**
Properties of the T-structure.</caption><tr><td>Width of rectangular thin rod</td><td>b<sub>1</sub></td><td>9.4 × 10<sup>-3</sup></td><td>m</td></tr><tr><td>Height of rectangular thin rod</td><td>h<sub>1</sub></td><td>80 × 10<sup>-3</sup></td><td>m</td></tr><tr><td>Width of shield-shaped part</td><td>b<sub>2</sub></td><td>80 × 10<sup>-3</sup></td><td>m</td></tr><tr><td>Height of shield-shaped part</td><td>h<sub>2</sub></td><td>40 × 10<sup>-3</sup></td><td>m</td></tr><tr><td>Total height of the structure</td><td>h<sub>s</sub> = h<sub>1</sub> + h<sub>2</sub></td><td>0.12</td><td>m</td></tr><tr><td>Thickness of structure</td><td>t<sub>s</sub></td><td>0.5 × 10<sup>-3</sup></td><td>m</td></tr><tr><td>Distance between structure and nozzle</td><td>d<sub>s</sub></td><td>0.25</td><td>m</td></tr><tr><td>Density (ρ = 20 °C)</td><td>ρ<sub>s</sub></td><td>8.435 × 10<sup>3</sup></td><td>kg/m<sup>3</sup></td></tr><tr><td>Young's modulus</td><td>E<sub>s</sub></td><td>93.7 × 10<sup>9</sup></td><td>N/m<sup>2</sup></td></tr><tr><td>Poisson's ratio</td><td>ν<sub>s</sub></td><td>0.34</td><td>—</td></tr></table>

## 3.2. Numerical parameters

### 3.2.1. Computational domain and boundary conditions

The complex wind tunnel geometry depicted in Fig. 2 is simplified to model the nozzle, the ground plate and the moving paddle as shown in Fig. 3. For this purpose, the geometry of the nozzle is reduced to a straight duct of length 1.5 m neglecting the nozzle's contraction ratio. The rectangular cross-section of the duct matches the nozzle's outlet in the experimental setup. A uniform velocity distribution with $u_{\text{inlet}}^{\text{steady}} = 4.914 \text{ m/s}$ is applied at the inlet patch. As the boundary layers develop along the walls of the duct, the flow reaches the desired free-stream

<footnote>a The values in brackets are the original values provided in Wood et al. [1] based on the prescribed values of the tooth belt axis in the experimental setup. They slightly deviate from the values adjusted in Boulbrachene and Breuer [13] and employed in this study to better match the real measured values in the experiment.</footnote>
<figure><img src="image_6.png" /><figcaption>Fig. 3. Geometry of the test case and definition of the boundary conditions.</figcaption></figure>

Source: Schematic sketch adapted from Boulbrachene and Breuer [13].

velocity $u_\infty = 5.14$ m/s at the outlet, consistent with the experimental conditions.

This description of the inlet velocity applies to the scenario where the wind tunnel is operated without WGG. However, Wood et al. [1] observed an undershoot of the streamwise velocity below the free-stream velocity at the nozzle's outlet during and after the paddle encounter. This drop in the inlet velocity is attributed to the paddle acting as a bluff body, obstructing the fluid flow and thus causing significant pressure losses. These increased pressure losses cannot be compensated by the blower in the short term. Consequently, less air is blown through the wind tunnel manifesting as an undershoot below the free-stream velocity. Further experimental tests in Wood et al. [1] demonstrated that the recovery time for the free-stream velocity to be reestablished is much longer than $\Delta t_{\text{gust}}$. To account for this inlet velocity drop, Constant Temperature Anemometry (CTA) measurements recorded inside the nozzle during the gust generation process were utilized in Boulbrachene and Breuer [13] to define a time-dependent inlet boundary condition. The specifics of the inlet velocity function including the CTA measurements on which it is based are detailed in Boulbrachene and Breuer [13]. Additionally, no-slip and impermeability boundary conditions are imposed on the walls of the duct, the ground plate, and the T-structure.

The duct is connected to a large container that extends 2 m in all directions from the duct's outlet. This container, measuring 2 m × 4.375 m × 4.5 m, includes the test section, which contains the ground plate and the mounted T-structure, as well as the immediate surrounding area of the open measuring section. To simulate the real wind tunnel environment, no-slip and impermeability boundary conditions are applied at the far-field walls of this container, except at a small outlet patch that models the receiver of the wind tunnel. At this outlet, a convective boundary condition is employed, with the convective velocity set to $u_{\text{conv}} = u_{\text{inlet}}^{\text{steady}}$, ensuring the proper outflow of air from the domain.

### 3.2.2. Computational grids

Given the geometric simplification of the wind tunnel's nozzle and the specific T-structure under consideration, the computational domain can be represented using a Cartesian grid, as long as the flexible structure remains undeformed by the flow field. This allows for the use of an initial Eulerian grid, which is a block-structured Cartesian grid consisting of approximately $117.5 \times 10^6$ control volumes (CVs). Compared to the preliminary validation study [13] of the flow in the empty test section without the T-structure, which showed good agreement between the data predicted by the immersed boundary method and the experimental data [1], the grid in the current study

comprises a factor of about 3 more CVs to account for an appropriate resolution of the flow around the T-structure.

The grid design and distribution of grid points across the computational domain in the current case are motivated by the following considerations. The geometric block-structure comprises 120 blocks, which are further subdivided into a parallel block-structure consisting of 942 blocks. This subdivision enables parallel execution across 942 computational cores. Grid refinement near relevant rigid walls (i.e., duct walls, ground plate and T-structure) is essential to capture the significant viscous effects leading to the formation of thin boundary layers. The goal is to position at least some nodes (i.e., the centers of the cells) within the viscous sublayer ($y^+ \le 5$). This careful refinement allows to accurately resolve the near-wall turbulence and boundary layer effects, which is critical for the fidelity of the simulation results. To ensure this, the following considerations were made.

The hydraulic diameter of the rectangular duct, calculated as $D_h = 4w_h/(2(w_h + h)) = 0.4286$ m, serves as the characteristic length for determining the Reynolds number $Re_{D_h} = \rho_{\text{air}} D_h u_\infty / \mu_{\text{air}} = 147,701$. Assuming a smooth pipe, the Darcy friction factor can be estimated based on the relation provided by Moody [51], yielding a skin friction coefficient $c_f$ of approximately $3.97 \times 10^{-3}$. To accurately resolve the near-wall region, the height of the first grid cell off the wall ($\Delta y$) is crucial. The grid spacing must ensure that the dimensionless wall distance $y_{cc}^+$ at the cell center is approximately 2, a value suitable for capturing near-wall effects. This can be calculated using the relation $\Delta y/D_h = (4/\text{Re}_{D_h}) \cdot \sqrt{2/c_f} \approx 6 \times 10^{-4}$, which results in a required grid spacing of $\Delta y \approx 2.6 \times 10^{-4}$ m. This grid spacing is applied along the walls of the duct, the ground plate, and in the vicinity of the T-structure (see Figs. 4(a) and (b)). This careful design of the grid allows the simulation to accurately capture the complex interactions between the fluid and the structure, particularly in regions of high shear stress and strong velocity gradients near solid boundaries.

To reduce the computational costs, the grid is mildly stretched in the off-wall directions with an expansion factor of about 1.05. Nearly the same expansion factor is also applied in the streamwise direction. Starting with the finest grid resolution at the outlet of the nozzle, where the paddle is located, the grid is mildly coarsened in upstream and downstream direction. Furthermore, in the vicinity of the flexible structure the same grid refinement is applied in all directions. Note that the current grid was already used in the preliminary study by Boulbrachene and Breuer [50] and the predicted results for the fluid flow induced by the gust generator but without considering the T-structure were validated against the experimental data by Wood et al. [1]. Beside the validation study [13], this is a further evidence that the presently applied grid is appropriate.
<figure><img src="image_7.png" /><figcaption>(a) Zoom of the nozzle outlet with the Eulerian grid and the uniformly distributed Lagrangian markers sketched in red (left). Zoom of the Eulerian grid (x-y-plane) surrounding the T-structure sketched in orange (right).</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>(b) Eulerian fluid grid on and surrounding the T-structure (y-z-plane).</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>(c) T-structure FE grid.</figcaption></figure>

<caption>Fig. 4. Details on the fluid and structure grids.</caption>

Based on the free-stream velocity, the stability of the explicit time integration scheme is ensured by choosing a time step size of $\Delta t = 10^{-5}$ s leading to a Courant–Friedrichs–Lewy number of less than 0.43. The same time step is set for the CSD solver.

For the structural dynamics simulation, the thin, flexible T-structure is discretized by 980 quadrilateral four-node shell elements, resulting in a total of 1080 nodes for the standard grid, as visible in Fig. 4(c). This choice is based on a grid sensitivity study detailed in Appendix B considering a static and a dynamic test case. Taking two finer grid levels into account, this investigation confirms that the standard grid leads to accurate results making further refinement and the associated increase of the already high computational effort superfluous. The elements used are special seven-parameters shell elements [52]. To avoid locking phenomena, the Enhanced Assumed Strain (EAS) method as proposed by Bischoff et al. [53] is employed. Specifically, the EAS parameters are set to (M, B, T, Q, S) = (0, 0, 4, 4, 0). The Rayleigh damping values mentioned in Section 2.3 are set to $\alpha_r = 0.47383$ Hz and $\beta_r = 2.1602 \times 10^{-8}$ s based on preliminary studies detailed in Appendix A.

### 3.2.3. Lagrangian markers

The geometry of the thin moving paddle, with a thickness of $3 \times 10^{-3}$ m, is simplified to a parallelepiped with right angles, measuring 0.3 m × 0.5 m. This rectangular cross-section penetrates the flow exiting the nozzle and is positioned 0.015 m away from the nozzle outlet in accordance with the experiment. The spacing between the Lagrangian markers that describe the paddle is carefully chosen based on the minimum local Eulerian grid spacing encountered by the paddle

during its closing and opening stroke (see Fig. 4(a)). This strategy is crucial to prevent any fluid from penetration through the surface of the paddle due to insufficient Eulerian cells coverage by the Lagrangian markers’ window function. In each coordinate direction, the minimum local Eulerian grid spacing is determined leading to the required number of uniformly distributed Lagrangian markers. As shown in Boulbrachene and Breuer [13], this procedure leads to a very high number of required Lagrangian markers, significantly slowing down the computations. This is because only a limited number of blocks contain Lagrangian markers and have to conduct the interpolation, spreading and tracking processes. A one-dimensional non-uniform distribution in spanwise direction allows to reduce the number of Lagrangian points from 6.7 million (uniform distribution) to about 1.57 million.

In this enhanced variant, the Lagrangian markers maintain a uniform distribution along the x- and y-axes. However, along the z-axis, the distribution is defined in such a way that for each cell center a Lagrangian marker is assigned, following the grid stretching employed along the z-axis. This strategy is chosen because, during the paddle’s opening and closing stroke, the distribution of Eulerian points along the z-axis does not change. By matching the Lagrangian marker distribution with the Eulerian grid stretching, the simulation ensures accurate interaction between the fluid and the paddle while minimizing unnecessary computational costs. This tailored distribution of markers optimizes the computational efficiency without compromising the fidelity of the simulation.

In addition to specifying the initial positions of the markers, the vertical velocity of the paddle as a function of time must also be
<figure><img src="image_10.png" /><figcaption>Fig. 5. Kinematic characteristics of the fitted gust type 1 paddle movement pattern.</figcaption></figure>

Source: Adapted from Boulbrachene and Breuer [13].

defined to describe its movement. For the gust type 1 considered in this study, Boulbrachene and Breuer [13] fitted a smooth (differentiable) analytical function to the vertical velocity curve derived from the experimental data, as shown in Fig. 5. In accordance with the experiment, an offset displacement of $\Delta y_{os} = 0.045$ m is included to account for the distance the paddle must travel from its idle position to start blocking the nozzle's outlet. This measure prevents unwanted initial interactions between the paddle and the undisturbed free-stream flow. Additionally, it facilitates a gradual buildup of the paddle's velocity prior to the interaction. A bell-shaped displacement pattern is prescribed, allowing the paddle to block a maximum of $\Delta y_{max} = 0.1437$ m of the nozzle's outlet, with a total interaction time of $\Delta t_{gust} = 0.352$ s. Following this, the paddle returns to its idle position. The present values of $\Delta y_{max}$ and $\Delta t_{gust}$ differ slightly from those reported in Wood et al. [1]. The reason is that the latter correspond to the prescribed values for the tooth belt axis in the experimental setup, whereas the present values more accurately reflect the actual measured motion pattern [13]. As described in Section 2.1, the displacement of the paddle is a derived quantity. Consequently, in the simulation, it is sufficient to prescribe the vertical velocity of the paddle.

# 4. Results and discussion

For the discussion of the results, the Cartesian coordinate system introduced in the zoomed view of Fig. 3 is used. The origin is located in the symmetry plane at the clamping point of the T-structure on the ground plate. The coordinates are normalized by the height $h_s$ of the entire T-structure, and the velocities by the free-stream velocity $u_\infty$. The time $t$ is made dimensionless by the gust generation time $\Delta t_{gust}$, as detailed in Section 3.1. The time instant, at which the paddle starts to obstruct the outlet of the nozzle initiating the gust generation process, serves as the reference point in time defining $t = 0$ for the subsequent analysis.

The discussion of the results is organized as follows: To understand the flow around the T-structure, Section 4.1 begins with an examination of the flow field and the structural deflections in the case without a gust as an additional validation of the FSI framework. The analysis of the wind gust and its impact on the structure throughout the entire process is presented in Section 4.2. For both scenarios, a comparison with the available experimental data is carried out.

## 4.1. Flow field and structural displacements without gust

In Breuer and Neumann [12], time-resolved measurements were repeated 104 times under identical operational conditions of the wind tunnel and wind gust generator, enabling phase-averaging of both the flow and the structural data. The data recorded prior to the gust injection, i.e., $t/\Delta t_{gust} \le 0$, are used to analyze the flow field without the influence of the gust. This part of the study can also be considered as an additional validation study for the ALE framework complementing the previously mentioned studies [13,50] for the case including the flexible T-structure.

As described in detail below, the deflection of the T-structure, and particularly the amplitude of its oscillations, is minimal in the absence of a gust. This allows the flow field to be time-averaged in the numerical simulation, facilitating comparison with the phase-averaged experimental data. For this purpose, the flow field was averaged over 1 million time steps, corresponding to a time interval of 10 s. Assuming that the T-structure approximately oscillates at its first dominant eigenfrequency $f_l^{num} = 6.826$ Hz, as determined for the bending mode in Appendix A, this time interval corresponds to approximately 68 oscillations, which is sufficient to ensure statistical convergence.

Fig. 6 presents contours of the averaged streamwise and vertical velocity fields in the geometric symmetry plane ($z = 0$). Additionally, streamlines calculated for the averaged flow in this plane are depicted for both the measured and the predicted data. The flow decelerates in front of the structure and deflects both in upward and downward direction, as evident in the contours of the vertical velocity component. The resulting stagnation region is clearly visible in the streamwise velocity contours. Both the rectangular thin rod and the shield-shaped section act as bluff bodies, disrupting the flow. Flow separation occurs at the sharp corners of the T-structure, generating a large recirculation region, which is especially prominent behind the shield-shaped part. A black contour line with the streamwise velocity $\bar{u} = 0$ marks the recirculation region in the symmetry plane, further highlighted by the streamlines of the time-averaged flow. The length of this backflow region is approximately 1.25 times the total height $h_s$ of the structure. Comparing the numerical and the experimental data depicted in Fig. 6, the predicted streamwise and vertical velocity components generally align well with the measured data. However, due to reflections of the laser light used for the particle-image velocimetry, velocities in the immediate vicinity of the ground wall could not be determined. Furthermore, the comparison clearly reveals that the statistics of the
<figure><img src="image_11.png" /><figcaption>**Fig. 6.** Comparison between measured and predicted data of the averaged flow field in the symmetry plane: streamwise velocity, vertical velocity, and streamlines.
*Source:* Experimental data taken from Breuer and Neumann [12].</figcaption></figure>

simulation data have converged better than the measurements. Considering that the experimentally determined mean values are based on only 104 independent measurements, while in the simulations the flow field was continuously averaged over a time interval of 10 s including 10⁵ samples, it is obvious why the predicted distributions are smoother than the measured data. This minor discrepancy arises due to the fact that the experiment was not designed to determine time-averaged quantities, but rather aimed to capture the instantaneous and statistically unsteady processes of wind gusts. Thus, small differences in the streamlines based on the averaged data can be expected. Despite this issue, the shape and size of the recirculation region in both data sets show good agreement.

The three-dimensional extension of the recirculation region can be
clearly seen in Fig. 7(a), which shows an iso-surface of the predicted
flow field based on $\bar{u} = 0$. Additionally, Figs. 7(b) and 7(c) depict
contours of the streamwise velocity and the pressure at three y-z-slices
in the wake ($x/h_s = 1/6, 2/3$, and $4/3$). At the slice $x/h_s = 2/3$, located
approximately in the middle of the backflow region, both the velocity
and pressure contours display a nearly circular shape, and the region
extends well beyond the width of the T-structure in the spanwise z-
direction. Also notable are the nearly circular wall streamlines around
the T-structure on the ground plate visible in Fig. 7, which form due
to the blockage of the flow field. A backflow area near the ground
plate, directly in front of the rectangular section, can also be observed.
This is attributed to the horseshoe vortex system forming in front of
bluff bodies, a phenomenon commonly observed in flows around wall-
mounted obstacles such as cubes [e.g., 18] or hemispheres [e.g., 54].
Surprisingly, the horseshoe vortex system forms relatively far upstream
of the T-structure, despite the thin profile of the rod and the minimal
blockage of the flow field near the bottom wall.

The good agreement between experiment and simulation is also
evident in the velocity profiles extracted along a vertical line in front
of the T-structure ($x/h_s = -1/6$) and three vertical lines in the wake
($x/h_s = 1/6, 2/3$, and $4/3$). Fig. 8 shows the time-averaged streamwise
and vertical velocity components in the symmetry plane. Four dotted
lines denote the four positions in the flow field where the velocities
are compared. For reference, the T-structure is drawn as a red line.
Each profile is plotted at its respective position which explains the
shift $+x/h_s$ used in the caption of Fig. 8. Due to the different ranges
of the two velocity components, the profiles are scaled with different
factors for improved visibility, e.g., 0.3 for the streamwise component
and 0.8 for the vertical component. Minor deviations are observed,
especially for the vertical component, which is significantly smaller
than the streamwise velocity.

The predicted Reynolds stresses in the symmetry plane are depicted
in Fig. 9, showing all three normal components and the most significant
shear stress component $\overline{u'u'}$. The strong shear layers forming in the
wake of the shield-shaped part of the T-structure are clearly visible,
particularly in the normal Reynolds stress component $\overline{u'u'}$ and the
Reynolds shear stress $\overline{u'v'}$. The normal vertical Reynolds stress $\overline{v'v'}$ is
notably large in the region where the end of the recirculation area is
located. In Breuer and Neumann [12], the Reynolds stresses presented
were limited to the two normal components ($\overline{u'u'}$, $\overline{v'v'}$) and the
shear stress component ($\overline{u'v'}$), as the 2D PIV measurements could only cap-
ture velocities within the plane of the laser light sheet.<fnref n="2" /> Furthermore,
the distributions were not very smooth due to the 104 measurements
available being insufficient for fully converged higher-order moments.

<footnote n="2">² Note the different coordinate system in Breuer and Neumann [12].</footnote>
<figure><img src="image_12.png" /><figcaption>Fig. 7. Predicted time-averaged three-dimensional flow field in the wake of the T-structure: Slices at $x/h_s = 1/6$, 2/3, and 4/3. Wall streamlines at the ground plate.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Fig. 8. Comparison of predicted and measured time-averaged velocity profiles at four positions in the flow field in front of and behind the T-structure depicted as a red line: $x/h_s = -1/6$, 1/6, 2/3, and 4/3.</figcaption></figure>

Source: Experimental data by Breuer and Neumann [12].

Consequently, a direct comparison between the experimental and simulated results is not meaningful. However, qualitatively, there is a high level of agreement between the simulation and the measurement data for the Reynolds stresses.

Finally, a comparison of the mean deformation between the experiment and the simulation is carried out. Fig. 10(a) shows the averaged deformation of the T-structure in the symmetry plane. The bending of the structure under the fluid load is clearly visible. Due to the small thickness of the T-structure, shear forces play a marginal role for the resulting fluid force in streamwise direction. In contrast, the stagnation pressure with a pressure coefficient $c_p$ close to one in the majority of the front face (Fig. 10(b)) and the low-pressure recirculation region with negative $c_p$ values in the back (Fig. 10(c)) lead to the dominance of the pressure forces. Nevertheless, the deflections remain overall small. For instance, the center of the shield-shaped part of the T-structure, denoted as point 1 in [12], deflects by only about 2.45 % of $h_s$, and for the upper edge of the structure, the deflection is approximately 3.2 %. As visible in Fig. 10(a) the predicted data agree well with the measurements. Note the different ranges of values assigned to the axes, which makes the

minor deviations appear much larger than they actually are, i.e., in the order of a few micrometers.

To analyze the deflection behavior of the T-structure for the flow case without gust, Fig. 11(a) presents selected excerpts from the time histories of the structural displacements of point 1 for both the measurement and the simulation. The experimental data clearly reveal two oscillation modes between which the coupled system switches back and forth. In the time period shown, a phase of about 4 s is initially visible, characterized by relatively large oscillation amplitudes. This is followed by a phase of about 6 s, in which the amplitudes have noticeably decreased. This phenomenon is succeeded by another and even longer phase with large amplitudes. Such an irregular switching between the modes with distinctly different amplitudes is also observed in the simulations. A thorough investigation of the deflection behavior, using animations of the predicted results, shows very clearly that there are two different vibration modes. Evidence for this statement is provided by Fast Fourier transformations of the deflections discussed below. If predominantly large vibration amplitudes are observed, the structure is in the bending-dominated mode. If, on the other hand, the amplitudes become small, the twisting mode becomes dominant,
<figure><img src="image_14.png" /><figcaption>Fig. 9. Predicted Reynolds stress components in the symmetry plane.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Fig. 10. (a) Comparison of predicted and measured time-averaged deformation of the T-structure in the symmetry plane; Time-averaged distribution of the pressure coefficient $c_p = (p - p_\infty)/(0.5 \rho_{\text{air}} v_\infty^2)$ on the front side (b) and the back side (c) of the structure including wall streamlines.</figcaption></figure>

Source: Experimental data from Breuer and Neumann [12].

though bending does not disappear completely. The oscillation mode changes in an unpredictable sequence, with both modes occasionally occurring simultaneously, as evident in the animations of the structure.

Fast Fourier transformations of the deflections for both simulated and measured data are depicted in Fig. 11(b). Two dominant frequencies are found ($f = 6.34$ Hz and $6.68$ Hz) in both data sets. These frequencies are close to the natural bending frequency of approximately 6.82 Hz, as determined in Appendix A. A detailed analysis of the predicted data reveals that the stronger peak at $f = 6.68$ Hz corresponds to the oscillation mode dominated by the bending of the structure. The slight reduction in frequency compared to the pure bending frequency in vacuum is attributed to the additional damping caused by the surrounding fluid. The second peak in the FFT at $f = 6.34$ Hz is associated to the state when the amplitudes of the oscillations are found to be significantly smaller and the twisting mode plays a prominent role in the kinematics of the T-structure. In this case, a FFT of the deflections of a point on the T-structure outside the symmetry plane depicted in Fig. 12 also shows a perceptible peak at the natural frequency of the second mode (see $f_2^{\text{num}}$ in Appendix A), i.e., the twisting mode.

It can be observed that the amplitudes of the measured oscillations are larger than those achieved in the simulations, likely due to the absence of inflow turbulence in the simulations. In principle, this could have been incorporated in the coupled predictions by using a synthetic turbulence generator, as done in previous studies such as Breuer [55] and De Nayer et al. [47,48,56]. However, implementing this option would require the knowledge of the characteristic properties of the incoming turbulence, such as the relevant length scales, which were not available for this setup. Consequently, this aspect was excluded in the initial approach.

## 4.2. Flow field and structural displacements with gust

The analysis of the flow field and structural displacements for the case without gust carried out in the previous section provides a description of the initial state before the gust impacts the structure. A large recirculation region is found behind the shield-shaped part, extending approximately $1.25 h_s$, and the maximum deflections remain quite small.
<figure><img src="image_16.png" /><figcaption>Fig. 11. Time histories and fast Fourier transformation of the total displacements of the center of the upper rectangular part of the T-structure. Predicted data and experimental data by Breuer and Neumann [12] for the case without gust.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Fig. 12. Fast Fourier transformation of the total displacements of a point (red circle) on the T-structure outside the symmetry plane.</figcaption></figure>

In the subsequent section, the more interesting case of the gust impact is analyzed. This scenario significantly differs from the previous one. In the former case, the incoming flow was steady, except for the very low turbulence intensity of approximately 0.2 % found in the wind tunnel by Wood et al. [54] and assumed to be valid in Breuer and Neumann [12]. This facilitated the analysis of the data based on time averaging. However, this is no longer possible in case of an active gust generator tackled in this section. The time-averaging approach is deemed invalid as the flow is statistically unsteady. Therefore, a phase-averaging procedure is required to obtain comparable numerical and experimental data.

In the experiment [12], this phase-averaging procedure relied on the synchronization of the measurement devices (PIV and DIC cameras) and the gust generator. The trigger signal for the PIV cameras was used to also trigger the paddle and the DIC cameras. Consequently, each recorded time step was comparable to the respective time step of subsequent gusts. The start of the paddle's motion was defined as $t=0$ for both PIV and DIC recordings. The same definition is applied in the numerical simulation. Thus, the temporal reference for phase averaging is the position of the paddle during the gust generation.

As mentioned in Section 4.1, the phase averaging in the experiments relied on 104 realizations. Indeed, cycle-to-cycle variations with standard deviations of the deflections of about 5 to 10 % of the phase-averaged total deflections were observed across different repetitions of the wind gust. Besides the (low) turbulence level of the flow within the wind tunnel, the timing of the gust initiation relative to the respective oscillation state is mainly made responsible for these variations. Due to the very high numerical effort involved in the coupled simulations, solely 32 predictions are used to determine the phase averages. For this purpose, the gust generation process is initiated starting from different,

randomly chosen time instants of the flow and structural deformation predicted by the simulations without gust.

At $t=0$, the paddle begins to obstruct the outlet of the nozzle, initiating the gust generation process. In the following, the effect of the gust on the flow field and the structural deformation is analyzed at different instants in time. A comparison between the simulation results and the measured data is presented for $t/\Delta t_{gust} = 0, 0.25, 0.5, 0.75, 1.0, 1.5, 2.0$, and $2.5$ (see Fig. 5 for the definition of the time instants). This time interval encompasses the gust generation time, the direct interaction time between the gust and the T-structure, and the phase when the gust exits the domain of interest.

Fig. 13 depicts a comparison between the measured data [12] and predicted data based on the phase-averaged streamwise velocity in the geometric symmetry plane. The flow field at $t/\Delta t_{gust} = 0$ is presented as a reference, using the same contour levels as in the subsequent figures. Furthermore, Fig. 14(a) shows the corresponding displacements at this time instant. At $t/\Delta t_{gust} = 0.25$, first alterations of the flow field due to the growing wind gust become visible in the region upstream and above the T-structure, whereas the region behind the obstacle is only marginally affected. The structure is also only minimally deflected at this stage (see Fig. 14(b)). At the third instant ($t/\Delta t_{gust} = 0.5$) when the paddle reaches its reversal point, the peak gust velocity is found in the region beneath the paddle. Consequently, significant changes occur for both the flow field (Fig. 13c) and the deflection of the structure (Fig. 14(c)). The size of the recirculation region behind the shield-shaped part and the structural deflections of the entire structure increase noticeably compared to the previous snapshot. The maximum deformation is, however, reached slightly later at $t/\Delta t_{gust} \approx 0.6$, when the upper edge of the T-structure deflects by approximately $0.1 h_s$. At the time instants when the paddle gradually unblocks the nozzle outlet, i.e., $t/\Delta t_{gust} = 0.75$ and $1.0$ depicted in Figs. 13d/13e and 14(d)/14(e), the previously described trends reverse. The recirculation region starts to shrink and the flow velocities gradually recover as the structure swings back. At approximately $t/\Delta t_{gust} = 1.0$ the recirculation region has the smallest size with a length of less than $1 h_s$. Afterwards, it increases again and slowly reaches its pre-gust level again. The same applies to the deflections, which aside from some weak oscillations, return to their pre-gust level. At about $t/\Delta t_{gust} = 2.5$, the original flow conditions and structural deflections are almost completely restored as visible in Figs. 13h and 14(h).

To further analyze the structural deformation during the gust impact, Fig. 15(a) presents time histories of the deflection of the center of the shield-shaped part of the T-structure (point 1). On the one hand, the instantaneous deformations predicted for ten individual gusts are shown. On the other hand, the corresponding phase-averaged data based on all predicted realizations (i.e., 32) are determined. Again $t=0$ marks the initialization of the gust generation, and the time
<figure><img src="image_18.png" /><figcaption>(a) $t/\Delta t_{\text{gust}} = 0$: Experiment (left) vs. Simulation (right).</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>(b) $t/\Delta t_{\text{gust}} = 0.25$: Experiment (left) vs. Simulation (right).</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>(c) $t/\Delta t_{\text{gust}} = 0.5$: Experiment (left) vs. Simulation (right).</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>(d) $t/\Delta t_{\text{gust}} = 0.75$: Experiment (left) vs. Simulation (right).</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>(e) $t/\Delta t_{\text{gust}} = 1.0$: Experiment (left) vs. Simulation (right).</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>(f) $t/\Delta t_{\text{gust}} = 1.5$: Experiment (left) vs. Simulation (right).</figcaption></figure>

<caption>**Fig. 13.** Comparison between measured and predicted data of the phase-averaged streamwise velocity in the symmetry plane at different time instants of the gust encounter.

*Source:* Experimental data taken from Breuer and Neumann [12].</caption>
<figure><img src="image_24.png" /><figcaption>Fig. 13. (continued).</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>Fig. 14. Comparison between measured and predicted data of the phase-averaged displacements in the symmetry plane at different time instants of the gust encounter.
Source: Experimental data taken from Breuer and Neumann [12].</figcaption></figure>

axis is normalized by the gust generation time $\Delta t_{\text{gust}}$. For comparison purposes, Fig. 15(b) depicts the corresponding experimental data. Note that the phase-averaged data in this case are based on a larger sample of individual gusts (i.e., 104).

Before the paddle enters the fluid stream ($t < 0$), minor flow-induced oscillations of the structure are visible. As discussed in Section 4.1 for the case without gust, the amplitudes of the measured oscillations are slightly larger than those observed in the simulations. The most likely cause for this deviation is the lack of inflow turbulence in the simulation as explained above.

Between the initiation of the gust and its first noticeable effect on the structure, a time delay of approximately 0.2 $\Delta t_{\text{gust}}$ is evident in the predicted data. This phenomenon was also noted in the measurements by Breuer and Neumann [12] and is attributed to the distance between

the nozzle's outlet and the T-structure, as well as the gust development time. Following this, the deflection increases rapidly, reaching its maximum value of about 0.07 $h_s$ at approximately $t/\Delta t_{\text{gust}} = 0.6$. Obviously, cycle-to-cycle variations are evident affecting the maximum deflection reached in each gust. In the experiment, slightly larger maximum deflections of approximately 40 µm are observed, leading to a marginally higher peak deflection in the phase-averaged data.

After reaching its maximum deflection, the structure swings back in the opposite direction. At $t/\Delta t_{\text{gust}} \approx 0.95$, the primary effect of the gust is largely alleviated, and the monitoring point returns to a deflection level similar to its pre-gust state. Due to the structure's elasticity and inertia, it oscillates beyond this state, reaching a minimum deformation of approximately 0.01 $h_s$ at about $t/\Delta t_{\text{gust}} = 1.3$. An interesting observation is that most gusts pass through a local maximum
<figure><img src="image_26.png" /><figcaption>**Fig. 15.** Time histories of the total displacements of the center of the upper rectangular part of the T-structure: Phase-averaged data (black thick line) and 10 singular gusts (colored lines).</figcaption></figure>

before reaching this minimum, with considerable variations from gust to gust. In the corresponding time interval of $0.9 \le t/\Delta t_{\text{gust}} \le 1.3$, the largest deviations are visible between the phase-averaged predicted displacement and its measured counterpart. Overall, during the initial rebound smaller deflections are seen in the simulation compared to the experiment. This discrepancy can be attributed to the absence of turbulence intensity in the inflow of the current computational framework. Additional preliminary predictions taking a certain level of background turbulence into account, confirm its influence on the amplitude of the T-structure's response. Nonetheless, the finally achieved minimum deflections in the simulation and the experiment are very similar.

Afterwards, at $t/\Delta t_{\text{gust}} > 1.5$, the oscillations gradually return to the status before the gust impact. The oscillation frequencies found in this time interval are similar to the dominant bending frequency of approximately $f = 6.68$ Hz observed in Section 4.1 for the case without gust. Therefore, the gust can also be seen as a kind of an external excitation.

Overall, the predicted time histories exhibit significant cycle-to-cycle variations, similar to those observed in the measurements. These

variations are primarily attributed to the pre-deflected position of the structure at the moment of the gust impact.

### 4.3. *Three-dimensional flow field induced by the gust*

To visualize the three-dimensional representation of the gust phenomenon, Fig. 16 illustrates pressure differences obtained by subtracting the pressure of the time-averaged flow results without gust from the phase-averaged pressure of the gust case. A negative iso-surface of this pressure difference highlights additional vortex cores induced by the gust.

During the paddle's downstroke a large horizontal vortex forms directly behind its base, where the gust is generated (see Fig. 16(a)). Due to the velocity gradient in the shear layer (high velocity at the nozzle's outlet and low velocity behind the paddle), the central region of this vortex rotates around the positive z-axis. A similar velocity gradient on either side of the nozzle outlet leads to the formation of vertical vortices, which later evolve into the 'feet' of a large arc-shaped vortex. Over time, this arc-type structure is convected downstream as it grows in size. At $t/\Delta t_{\text{gust}} = 0.5$, when the paddle reaches its maximum
blocking ratio, the arc-shaped vortex, marking the gust front, passes beyond the T-structure (see Fig. 16(b)). The feet of this structure, symmetrically positioned on both sides of the bottom plate, rotate around the y-axis. This large arc-like structure is connected to the paddle's sides through two prominent side vortices formed at the paddle's lower corners and rotating around the x-axis.

As expected, the shear layer behind the paddle and the gust front contains a series of Kelvin–Helmholtz vortices. Furthermore, the symmetric vortices at the top of the strong front vortex adopt a horn-like shape. These vortices are likely resulting from the interaction between the y-axis, the x-axis and the z-axis rotating vortices (see Fig. 16(b)). Fluid particles flowing close to the nozzle's sides are captured by the y-axis rotating vortices and redirected toward the center of the test section. Some of these fluid particles, located near the x-axis rotating vortices, are accelerated upward and drawn into the counter-rotating flow generated by the large z-axis rotating vortex. The clockwise rotation of the large front vortex pulls the fluid downward toward the floor, and, in combination with the main flow, funnels it toward the center of the test section, giving rise to the elongated horn-like vortices.

As the paddle moves upward, returning to its original position, the arc-shaped vortical structure descends into the wake while remaining connected to the paddle via the side vortices, which gradually weaken. The large z-rotating vortex defining the gust front is initially quasi-2D (see Fig. 16(a)). At later stages (see Figs. 16(b) and 16(c)), the two-dimensionality is disrupted. This is attributed to the deceleration effect that the stationary fluid surrounding the ground plate has on the generated gust. Therefore, the center of the z-rotating vortex (undergoing less resistance) progresses faster. Moreover, the horn-like vortices above the arc-shaped vortex begin to dissipate around $t/\Delta t_{\text{gust}} = 0.75$. Overall, Fig. 16 illustrates the complexity of the flow dynamics induced by the wind gust.

## 5. Conclusions and outlook

Combining different numerical methods in an appropriate way, which at first glance do not seem to belong together, can be useful for simulating complex coupled fluid-structure interactions. In this case, the focus is on investigating a flexible structure in a wind tunnel, subjected to highly dynamic loads generated by a custom-designed wind gust generator. The simulation of the turbulent flow around the T-structure and its dynamic deformation is carried out using a partitioned method based on the ALE approach on general, curvilinear, body-fitted grids, a technique that has been widely applied to similar problems. The gust generation process, which involves a rigid paddle moving through the flow field, is modeled using the IB method with direct forcing. Despite the paddle causing significant disruption to the free-stream flow in the wind tunnel and inducing a rapid, short-term increase in the main flow velocity, the IB method ensures smooth pressure and velocity fields without any numerically induced oscillations. This enables a precise reproduction of the real gust generator in the wind tunnel and a thorough analysis of the actual FSI problem including three-dimensional visualizations of the entire flow field not available in the experiments. Furthermore, it allows for a meaningful comparison of the simulation results with the experimental data reported by Breuer and Neumann [12]. The main findings can be summarized as follows:

* In the case without gust, the deformation of the T-structure remains small. Furthermore, marginal bending and even weaker torsional oscillations of the structure were observed, which allow averaging of the fluid flow and structure in time. The predicted flow field is found to be in close agreement with the PIV measurements of Breuer and Neumann [12]. The same applies to the predicted structural deformations and the DIC measurements.

* Both simulation and experiment exhibit a large recirculation region behind the shield-shaped section and strong shear layers

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

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

<figure><img src="image_29.png" /><figcaption>**Fig. 16.** Three-dimensional visualization of the gust effect on the flow field: Iso-surfaces of the pressure differences between the time-averaged pressure without the gust and the phase-averaged pressure in case of the gust.</figcaption></figure>

occurring around this backflow area. The maximal structural deformations are in the range of 3 to 4 mm.

* Both measurements and simulations exhibit two distinct oscillation modes of the structure, between which the coupled system alternates irregularly. Relatively large amplitudes are observed in the bending-dominated mode, while noticeably smaller amplitudes occur in the mode where twisting also plays a role. Both corresponding bending frequencies are close to the natural bending frequency in vacuum.

* Wind gusts significantly change the flow field around the structure. In close agreement between simulation and experiment,
the recirculation region behind the shield-shaped part increases noticeably to about 1.45 $h_s$ when the velocity of the approaching flow rises due to the wind gust. A maximum deformation of the upper edge of the T-structure of about 0.1 $h_s$ is approximately reached at $t/\Delta t_{\text{gust}} = 0.6$.

* When the paddle unblocks the nozzle's outlet and the gust starts to diminish, the recirculation region shrinks significantly reaching a minimum size of about 0.9 $h_s$ at approximately $t/\Delta t_{\text{gust}} = 1.0$. Simultaneously, the structure moves against the main flow direction. After a few bending oscillations the flow field as well as the structural behavior of the pre-gust levels are recovered.

* Cycle-to-cycle variations play a role for both cases, with and without wind gusts. For the gust case, these variations are mainly attributed to the pre-deflected position of the structure at the moment of the gust impact. A second factor relevant to the gust and the no-gust case is the turbulence intensity of the flow exiting the wind tunnel nozzle. Since the inflow turbulence was not taken into account in the simulations for simplification, the oscillation amplitudes are found to be slightly smaller than those in the experiments. Additional preliminary simulations taking the inflow turbulence into account confirm that the amplitudes of the structural oscillations rise with an increasing turbulence intensity level. Further studies in the future shall quantify this effect.

The study demonstrates that the simulation methodology developed works well, successfully replicating both the gust generation process in the wind tunnel induced by the moving paddle and the coupled FSI problem of the thin structure interacting with the strongly varying turbulent fluid flow with maximum inlet velocities rising by more than 60 %. In future studies, this methodology will enable comprehensive investigations of complex fluid-structure interaction problems, especially those involving membrane structures subjected to predefined wind gust loads.

## CRediT authorship contribution statement

M. Breuer: Writing - review & editing, Writing - original draft, Validation, Supervision, Resources, Project administration, Investigation, Funding acquisition, Conceptualization. K. Boulbrachene: Validation, Software, Methodology, Investigation. G. De Nayer: Visualization, Validation, Software, Methodology, Investigation, Formal analysis, Data curation.

## 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

Computational resources on the HPC cluster HSUper have been provided by the project hpc.bw, funded by dtec.bw - Digitalization and Technology Research Center of the Bundeswehr. dtec.bw is funded by the European Union - NextGenerationEU. Furthermore, we would like to thank Dr. J.N. Wood (HSU Hamburg) for providing the displacement measurements by the laser vibrometer mentioned in Appendix A.

## Appendix A. Determination of the Young's modulus and material damping parameters

To determine the parameters for the material model used in this study, two simple structural test cases were conducted in the laboratory and numerically simulated using the Carat++ solver (see Fig. A.17). The first test involved bending the T-shaped structure and then releasing it in a vacuum environment. The second considered the same

but for the twisting mode. The displacements at point 1, as defined in [12], were measured using a laser vibrometer (see Fig. A.17(a)). This measurement point is located at the center of the upper rectangular part of the T-structure.

The structural simulation employed the finite-element (FE) mesh described in Section 3.2.2, with the dimensions of the T-shaped structure corresponding to those used in the experimental study by Breuer and Neumann [12]. The material density $\rho_s$ was set to $\rho_s = 8435 \text{ kg/m}^3$, as provided in [12]. The brass material (alloy CuZn37) was modeled using the St. Venant-Kirchhoff law, with a typical Poisson's ratio of 0.34 for brass. Using a Young's modulus of $E_s = 110 \text{ GPa}$, as initially suggested by Breuer and Neumann [12], led to an overestimation of the first dominant eigenfrequency associated with the bending mode of the T-shaped structure. To achieve a more accurate match, a reduced Young's modulus of $E_s = 93.7 \text{ GPa}$ was determined through a trial-and-error process based on the dynamic bending test. In these tests the gravitational acceleration was taken into account but the structural damping was first of all ignored. The adjustment of $E_s$ resulted in a natural bending frequency of $f_1^{\text{num}} = 6.814 \text{ Hz}$, closely aligning with the experimental bending frequency of $f_1^{\text{exp}} = 6.82 \text{ Hz}$. The frequency of the second mode, determined by the pure twisting case around the vertical y-axis, also agrees reasonably well between the simulation ($f_2^{\text{num}} = 23.75 \text{ Hz}$) and the experiment ($f_2^{\text{exp}} = 22.45 \text{ Hz}$).

As visible in the experimental data in Fig. A.17(a), the structure exhibits structural damping. To account for the damping, the Rayleigh damping model is employed in this study. The mass- and stiffness-related Rayleigh damping coefficients, $\alpha_r$ and $\beta_r$, were determined using an automated search process with the surrogate-based MATLAB optimization function surrogateopt. The objective was to minimize the $L^2$-norm of the error between the logarithmic decrement of the experimental and numerical deflection data at point 1. To isolate the pure bending mode in the experimental data, a discrete Fourier transformation was performed on the temporal response of the measured data. A rectangular window function was then applied to isolate the frequency response associated with the bending frequency. Finally, a clean signal in the time domain (depicted by the red curve in Fig. A.17(b)) was obtained via an inverse Fourier transformation. The numerical optimization process was configured to perform a maximum of 200 evaluations, with 12 evaluations running in parallel on a computing cluster. The bounds for the optimization were set to $0.4 \text{ Hz} \le \alpha_r \le 0.5 \text{ Hz}$ and $10^{-8} \text{ s} \le \beta_r \le 10^{-6} \text{ s}$, respectively. Each evaluation involved setting up the simulation, running the structural solver, and post-processing the results to evaluate the error in the logarithmic decrement. The procedure yielded several sets of $\alpha_r$ and $\beta_r$ values, that resulted in a very good agreement between the predicted and the measured data of the bending deflections. For example, Fig. A.17(b) shows the results achieved with $\alpha_r = 0.47383 \text{ Hz}$ and $\beta_r = 2.1602 \times 10^{-8} \text{ s}$. The structural damping in these cases has only a marginal influence on the bending frequency of the dynamic simulation. Thus, the simulated value of $f_1^{\text{num}} = 6.826 \text{ Hz}$ remains quasi unchanged and in close agreement with the experimental value of $f_1^{\text{exp}} = 6.82 \text{ Hz}$.

To obtain the final values from the possible parameter combinations ($\alpha_r, \beta_r$) resulting from the bending test, the twisting behavior was also taken into account. Although twisting plays a less important role than bending in the present case, it should still be considered in order to represent the physics as good as possible. In the experiment and in the simulation, this second vibration mode was therefore also excited and its decay behavior was analyzed. For parameter combinations that had shown almost identical results for the first mode, there were noticeable differences in the decay behavior of the second mode. The final selection from this pool of best-fit parameter combinations for bending was made based on the best fit with the experimental twisting decay curve resulting in $\alpha_r = 0.47383 \text{ Hz}$ and $\beta_r = 2.1602 \times 10^{-8} \text{ s}$.
<figure><img src="image_30.png" /></figure>

(a) Experimental displacements measured by the laser vibrometer.

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

(b) Displacements predicted by the FEM simulation normalized by the initial displacement and comparison with the experimental signal limited to the bending mode.

<caption>Fig. A.17. Free-oscillation test in vacuum: Determination of the material parameters of the model.</caption>

<table><caption>**Table B.3**

Static grid sensitivity study: Total displacements *d* of the center of the shield-shaped part of the T-structure for different grids under increasing wind load.</caption><thead><tr><th colspan="2">Grid</th><th colspan="2">Standard</th><th colspan="2">Fine</th><th>Reference</th></tr><tr><th>Wind velocity [m/s]</th><th>Wind load [N/m²]</th><th>d [m]</th><th>Error [%]</th><th>d [m]</th><th>Error [%]</th><th>d [m]</th></tr></thead><tbody><tr><td>5.14</td><td>16</td><td>2.0593 × 10<sup>-3</sup></td><td>0.42</td><td>2.0608 × 10<sup>-3</sup></td><td>0.35</td><td>2.0680 × 10<sup>-3</sup></td></tr><tr><td>8.40</td><td>43</td><td>5.5212 × 10<sup>-3</sup></td><td>0.42</td><td>5.5252 × 10<sup>-3</sup></td><td>0.35</td><td>5.5446 × 10<sup>-3</sup></td></tr><tr><td>12.8</td><td>100</td><td>12.6886 × 10<sup>-3</sup></td><td>0.42</td><td>12.6984 × 10<sup>-3</sup></td><td>0.34</td><td>12.7422 × 10<sup>-3</sup></td></tr></tbody></table>

## Appendix B. Grid sensitivity study for the structure model

A special seven-parameters shell element is used for the finite-element discretization of the T-structure. The material parameters determined in Appendix A are applied. Three levels of refinement are investigated under uniformly applied wind loads of 16, 43 and 100 N/m², resulting from wind velocities of $u_{inflow} = u_\infty = 5.14$ m/s, $1.6 u_\infty = 8.40$ m/s and $2.5 u_\infty = 12.8$ m/s. These exemplary loads are calculated based on the assumption that the static pressure at the front of the body is approximately $0.5 \rho_{air} u_{inflow}^2$. The resulting total displacements *d* at point 1, i.e., the center of the shield-shaped part of the T-structure, are given in Table B.3. The load-controlled geometrically non-linear analysis yields very similar results for the standard, fine, and reference grids, consisting of 980, 4000, and 16,000 quadrilateral shell elements, respectively. This table also includes the relative errors with respect to the results obtained on the reference grid. It is evident that spatial convergence has already been achieved with the standard grid.

To assess the influence of the grid on the dynamic response of the flexible T-structure, the structure is first deformed and then released in a vacuum. In the first configuration, the structure oscillates at a frequency of $f_1^{num}$ in its first dominant mode, the bending mode. In the second configuration, it oscillates at a frequency of $f_2^{num}$ in its second mode, the twisting mode. Table B.4 summarizes these frequencies, which were determined by fast Fourier transformations of the total displacement at the center of the shield-shaped part of the T-structure. The convergence of both frequencies is achieved with the standard grid.

### Data availability

Data will be made available on request.

<table><caption>**Table B.4**

Dynamic grid sensitivity study: First and second dominant frequencies associated with the total displacement of the center of the shield-shaped part of the T-structure.</caption><thead><tr><th colspan="2">Grid</th><th colspan="2">Standard</th><th colspan="2">Fine</th><th>Reference</th></tr><tr><th>Frequencies</th><th>f [Hz]</th><th>Error [%]</th><th>f [Hz]</th><th>Error [%]</th><th>f [Hz]</th><th></th></tr></thead><tbody><tr><td>f<sub>1</sub><sup>num</sup></td><td>6.826</td><td>0.18</td><td>6.814</td><td>0.0</td><td>6.814</td><td></td></tr><tr><td>f<sub>2</sub><sup>num</sup></td><td>23.78</td><td>0.93</td><td>23.67</td><td>0.47</td><td>23.56</td><td></td></tr></tbody></table>

### References

[1] Wood JN, Breuer M, Neumann T. A novel approach for artificially generating horizontal wind gusts based on a movable plate: The paddle. J Wind Eng Ind Aerodyn 2022;230:105170. http://dx.doi.org/10.1016/j.jweia.2022.105170.

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

[3] Zaide A, Raveh DE. Numerical simulation and reduced-order modeling of airfoil gust response. AIAA J 2006;44(8):1826-34. http://dx.doi.org/10.2514/1.16995.

[4] Wales C, Jones D, Gaitonde A. Prescribed velocity method for simulation of aerofoil gust responses. J Aircr 2014;52(1):64-76. http://dx.doi.org/10.2514/1.C032597.

[5] Badrya C, Biler H, Jones AR, Baeder JD. Effect of gust width on flat-plate response in large transverse gust. AIAA J 2021;59(1):49-64. http://dx.doi.org/10.2514/1.J059678.

[6] IEC-Standard. Measurement and assessment of power quality of grid connected wind turbines. Technical report 61400-21. Tech. rep., 2002.

[7] Bos R. Extreme gusts and their role in wind turbine design (Ph.D. thesis), Delft University of Technology, The Netherlands; 2017, http://dx.doi.org/10.4233/uuid:d6097e3a-1cdd-4845-a71c-90f469d28b7a.

[8] Storey RC, Norris SE, Cater JE. Modelling turbine loads during an extreme coherent gust using large eddy simulation. J Phys: Conf Ser 2014;524:012177. http://dx.doi.org/10.1088/1742-6596/524/1/012177.

[9] Onol AO, Yesilyurt S. Effects of wind gusts on a vertical axis wind turbine with high solidity. J Wind Eng Ind Aerodyn 2017;162:1-11. http://dx.doi.org/10.1016/j.jweia.2017.01.003.
[10] Menegozzo L, Dal Monte A, Benini E, Benato A. Small wind turbines: A numerical study for aerodynamic performance assessment under gust conditions. Renew Energy 2018;121:123-32. http://dx.doi.org/10.1016/j.renene.2017.12.086.

[11] Wood JN, Breuer M. The Paddle: A novel procedure for artificial gust generation in a wind tunnel. In: Kähler CJ, Fuchs T, Hain R, Scharnowski S, Ruck B, Leder A, editors. Proc. der 30. Gala-Fachtagung 'Experimentelle Strömungsmechanik', Sept. 7-9, 2023. München; 2023, 37-1-37-8.

[12] Breuer M, Neumann T. Synchronous high-speed measurements of a flexible structure under wind gust load. Phys Fluids 2024;36(7):075147. http://dx.doi.org/10.1063/5.0215724.

[13] Boulbrachene K, Breuer M. A numerical method to mimic an experimental wind gust generator: The immersed boundary gust generator. Phys Fluids 2024;36(1):015146. http://dx.doi.org/10.1063/5.0187787.

[14] Durst F, Schäfer M. A parallel block-structured multigrid method for the prediction of incompressible flows. Int J Numer Methods Fluids 1996;22(6):549-65. http://dx.doi.org/10.1002/(SICI)1097-0363(19960330)22:6<549::AID-FLD366>3.0.CO;2-7.

[15] Breuer M, De Nayer G, Münsch M, Gallinger T, Wüchner R. Fluid-structure interaction using a partitioned semi-implicit predictor-corrector coupling scheme for the application of large-eddy simulation. J Fluids Struct 2012;29:107-30. http://dx.doi.org/10.1016/j.jfluidstructs.2011.09.003.

[16] Ferziger JH, Perić M. Computational methods for fluid dynamics. 3rd ed. Springer Berlin; 2002.

[17] Rhie CM, Chow WL. Numerical study of the turbulent flow past an airfoil with trailing-edge separation. AIAA J 1983;21(11):1525-32. http://dx.doi.org/10.2514/3.8284.

[18] Breuer M. Direkte Numerische Simulation und Large-eddy Simulation turbulenter Strömungen auf Hochleistungsrechnern. Habilitationsschrift, Universität Erlangen-Nürnberg, Berichte aus der Strömungstechnik, Aachen: Shaker Verlag; 2002.

[19] Smagorinsky J. General circulation experiments with the primitive equations I: The basic experiment. Mon Weather Rev 1963;91(3):99-165. http://dx.doi.org/10.1175/1520-0493(1963)091<0099>GCEWTP>2.3.CO;2.

[20] Thompson JF, Warsi ZUA, Wayne Mastin C. Numerical grid generation: Foundations and applications, vol. 45, North-Holland; 1985.

[21] Knupp PM. Algebraic mesh quality metrics for unstructured initial meshes. Finite Elem Anal Des 2003;39(3):217-41.

[22] Sen S, De Nayer G, Breuer M. A fast and robust hybrid method for block-structured mesh deformation with emphasis on FSI-LES applications. Int J Numer Methods Eng 2017;111(3):273-300. http://dx.doi.org/10.1002/nme.5465.

[23] Hoblit FM. Gust loads on aircraft: Concepts and applications. AIAA Education Series, Washington D.C.; 1988, http://dx.doi.org/10.2514/4.861888.

[24] Wu Z, Cao Y, Ismail M. Gust loads on aircraft. Aeronaut J 2019;123(1266):1216-74. http://dx.doi.org/10.1017/aer.2019.48.

[25] Norris SE, Cater JE, Stol KA, Unsworth CP. Wind turbine wake modelling using large-eddy simulation. In: Proceedings of the 17th Australian Fluid Mechanics Conference. University of Auckland, Australia: Curran Associates, Inc.; 2010.

[26] De Nayer G, Breuer M, Perali P, Grollmann K. Modeling of wind gusts for large-eddy simulations related to fluid-structure interactions. In: Salvetti MV, Armenio V, Fröhlich J, Geurts BJ, Kuerten H, editors. ERCOFTAC Series, Direct and Large-eddy Simulation XI, 11th Int. ERCOFTAC Workshop on Direct and Large-eddy Simulation: DLES-11. Pisa, Italy, May 29-31, 2017. 25, Springer Nature Switzerland AG; 2019, p. 453-9. http://dx.doi.org/10.1007/978-3-030-04915-7_60.

[27] De Nayer G, Breuer M. A source-term formulation for injecting wind gusts in CFD simulations. J Wind Eng Ind Aerodyn 2020;207:104405. http://dx.doi.org/10.1016/j.jweia.2020.104405.

[28] Parameswaran V, Baeder JD. Indicial aerodynamics in compressible flow-Direct computational fluid dynamic calculations. J Aircr 1997;34(1):131-3. http://dx.doi.org/10.2514/2.2146.

[29] Singh R, Baeder JD. Direct calculation of three-dimensional indicial lift response using computational fluid dynamics. J Aircr 1997;34(4):465-71. http://dx.doi.org/10.2514/2.2214.

[30] Huntley SJ, Jones D, Gaitonde A. Aeroelastic gust response of an aircraft using a prescribed velocity method in viscous flows. In: 23rd AIAA Computational Fluid Dynamics Conference. Denver, Colorado; 2017.

[31] Biler H, Badrya C, Jones AR. Experimental and computational investigation of transverse gust encounters. AIAA J 2019;57(11):4608-22. http://dx.doi.org/10.2514/1.J057646.

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

[33] Boulbrachene K, De Nayer G, Breuer M. Assessment of two wind gust injection methods: Field velocity vs. split velocity method. J Wind Eng Ind Aerodyn 2021;218:104790. http://dx.doi.org/10.1016/j.jweia.2021.104790.

[34] Uhlmann M. An immersed boundary method with direct forcing for the simulation of particulate flows. J Comput Phys 2005;209(2):448-76. http://dx.doi.org/10.1016/j.jcp.2005.03.017.

[35] Peskin CS. The immersed boundary method. Acta Numer 2002;11:479-517. http://dx.doi.org/10.1017/S0962492902000077.

[36] Roma AM, Peskin CS, Berger MJ. An adaptive version of the immersed boundary method. J Comput Phys 1999;153(2):509-34. http://dx.doi.org/10.1006/jcp.1999.6293.

[37] Pinelli A, Naqavi I, Piomelli U, Favier J. Immersed-boundary methods for general finite-difference and finite-volume Navier-Stokes solvers. J Comput Phys 2010;229(24):9073-91. http://dx.doi.org/10.1016/j.jcp.2010.08.021.

[38] Anderson E, Bai Z, Bischof C, Blackford S, Demmel J, Dongarra J, Du Croz J, Greenbaum A, Hammarling S, McKenney A, Sorensen D. LAPACK users' guide. 3rd ed., Philadelphia, PA: Society for Industrial and Applied Mathematics; 1999.

[39] Schäfer F, Breuer M. Comparison of c-space and p-space particle tracing schemes on high-performance computers: Accuracy and performance. Int J Numer Methods Fluids 2002;39(4):277-99. http://dx.doi.org/10.1002/fld.324.

[40] Breuer M, Baytekin HT, Matida EA. Prediction of aerosol deposition in 90 degrees bends using LES and an efficient Lagrangian tracking method. J Aerosol Sci 2006;37(11):1407-28. http://dx.doi.org/10.1016/j.jaerosci.2006.01.013.

[41] Breitenberger M, Apostolatos A, Philipp B, Wüchner R, Bletzinger K-U. Analysis in computer aided design: Nonlinear isogeometric B-Rep analysis of shell structures. Comput Methods Appl Mech Engrg 2015;284:401-57. http://dx.doi.org/10.1016/j.cma.2014.09.033.

[42] Philipp B, Wüchner R, Bletzinger K-U. Advances in the form-finding of structural membranes. Procedia Eng 2016;155:332-41. http://dx.doi.org/10.1016/j.proeng.2016.08.036.

[43] Basar Y, Weichert D. Nonlinear continuum mechanics of solids: Fundamental mathematical and physical concepts. Springer Science & Business Media; 2013.

[44] Sicklinger S, Belsky V, Engelmann B, Elmqvist H, Olsson H, Wüchner R, Bletzinger K-U. Interface Jacobian-based co-simulation. Int J Numer Methods Eng 2014;98:418-44. http://dx.doi.org/10.1002/nme.4637.

[45] Wang T, Wüchner R, Sicklinger S, Bletzinger K-U. Assessment and improvement of mapping algorithms for non-matching meshes and geometries in computational FSI. Comput Mech 2016;57(5):793-816. http://dx.doi.org/10.1007/s00466-016-1262-6.

[46] Apostolatos A, De Nayer G, Bletzinger K-U, Breuer M, Wüchner R. Systematic evaluation of the interface description for fluid-structure interaction simulations using the isogeometric mortar-based mapping. J Fluids Struct 1999;86:368-99. http://dx.doi.org/10.1016/j.jfluidstructs.2019.02.012.

[47] De Nayer G, Apostolatos A, Wood JN, Bletzinger KU, Wüchner R, Breuer M. Numerical studies on the instantaneous fluid-structure interaction of an air-inflated flexible membrane in turbulent flow. J Fluids Struct 2018;82:577-609. http://dx.doi.org/10.1016/j.jfluidstructs.2018.08.005.

[48] De Nayer G, Breuer M, Boulbrachene K. FSI simulations of wind gusts impacting an air-inflated flexible membrane at Re=100,000. J Fluids Struct 2022;109:103462. http://dx.doi.org/10.1016/j.jfluidstructs.2021.103462.

[49] De Nayer G, Breuer M. Assessment of discrete wind gust parameters: Towards the worst-case scenario of a FSI application in form of an inflated hemisphere. J Wind Eng Ind Aerodyn 2022;231:105207. http://dx.doi.org/10.1016/j.jweia.2022.105207.

[50] Boulbrachene K, Breuer M. Application of an immersed boundary method to predict the flow in an experimental wind gust generator: An LES study. In: 14th Int. ERCOFTAC Workshop on DNS and LES: DLES-14, Erlangen, Germany, April 10-12, 2024.

[51] Moody LF. An approximate formula for pipe friction factors. Trans ASME 1947;69(12):1005-11.

[52] Bischoff M. Theorie und Numerik einer dreidimensionalen Schalenformulierung (Ph.D. thesis), Institut für Baustatik, Universität Stuttgart, Germany; 1999.

[53] Bischoff M, Wall WA, Bletzinger K-U, Ramm E. Models and finite elements for thin-walled structures. In: Stein E, De Borst R, Hughes TJR, editors. Encyclopedia of computational mechanics. Vol. 2, Chichester: John Wiley & Sons Ltd; 2004, p. 59-138.

[54] Wood JN, De Nayer G, Schmidt S, Breuer M. Experimental investigation and Large-Eddy simulation of the turbulent flow past a smooth and rigid hemisphere. J Flow Turbul Combust 2016;97(1):79-119. http://dx.doi.org/10.1007/s10494-015-9690-5.

[55] Breuer M. Effect of inflow turbulence on an airfoil flow with laminar separation bubble: An LES study. J Flow Turbul Combust 2018;101(2):433-56. http://dx.doi.org/10.1007/s10494-017-9890-2.

[56] De Nayer G, Schmidt S, Wood JN, Breuer M. Enhanced injection method for synthetically generated turbulence within the flow domain of eddy-resolving simulations. Comput Math Appl 2018;75(7):2338-55. http://dx.doi.org/10.1016/j.camwa.2017.12.012.