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

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


Introduction
Singular wind events such as tornados, downbursts or wind gusts are not only of interest for meteorologists who are trying to analyze, understand and foresee such extreme phenomena (e.g., Feser et al., 2015). They have also garnered the attention of engineers due to their devastating consequences on structures. An impressive example of extreme structural damages on wind turbines observed at a wind farm in 2011 was investigated in details by Hawbecker et al. (2017). A wind turbine tower was buckled and blades of several others turbines were stripped away during a storm. Simulations including transient wind gusts were performed and compared to available wind datasets. Especially constructions made of lightweight thin structures (e.g., large umbrellas, tents, stadium roofs) are vulnerable to wind events. The National Institute of Standards and Technology (2020) (NIST, USA) gathers disaster studies related to windstorms in the USA since 1969. A variety of similar reports are available. For example, Yokel et al. (1976) present the wind damages on buildings in the metropolitan area of Washington DC during a storm. The report compiles the different destructions but also the meteorological data. More recently in 2009, the Dallas Cowboys indoor practice facility, a membrane-covered frame structure, collapsed during a severe thunderstorm. The facility including the structural design and analysis, its destruction and the wind environment during the storm are described in detail in Gross et al. (2013). Massive wind destructions also occur in nature, particularly forests leading to typical damage patterns. Consequently, wind tunnel measurements on the flow and the gust dynamics as well as studies on the modeling of these phenomena can be found in the literature (see, e.g., Gromke and Ruck, 2018;Kamimura et al., 2019;Tischmacher and Ruck, 2013).
In order to prevent structural damages due to the wind, design standards are prescribed. For example, the engineering handbook of Frost et al. (1978), the IEC-Standard 61400-21 IEC-Standard (2002) and the wind energy handbook of Burton et al. (2001) provide guidelines for the construction of wind turbines. These are often based on measurements and estimations such as in the work of Kasperski (2007). However, the assumptions made in the guidelines are rather simplistic. Therefore, thorough investigations of the proposed design submitted to different wind conditions either in experimental facilities and/or by the application of sophisticated highly-resolved numerical simulations is a must for a technological advancement in the future.
Due to the extremely wide field the following considerations are restricted on the one hand to numerical simulations and on the other hand to wind gusts occurring in the micro scale, e.g., from meters to one hundred meters (Steyn et al., 1981). Methods to describe these wind gusts can be classified into deterministic and stochastic models.
In the deterministic model-based method the mean shape of the gust is defined by a mathematical function. Two well-known shapes were proposed in the IEC-Standard (2002) and prescribed by the Federal Aviation Regulation: The Extreme Operating Gust (EOG) also called Mexican Hat and the Extreme Coherent Gust (ECG) also denoted 1-Cosine shape. A comparison of selected real gust events with the idealized EOG is available in Emeis (2018). A mean shape in form of a Gaussian distribution was also recently applied by Yawar et al. (2019). Another more realistic (one-and two-dimensional) model was proposed by Knigge and Raasch (2016) derived from turbulent flow fields computed by large-eddy simulations (LES). In general, the deterministic representation of a gust by a prescribed mathematical function is easy to realize computationally, but also simplistic. Thus, it often does not fully represent the reality (Bos, 2017;Branlard, 2009).
The stochastic model-based method tries to represent the gust in a more physical manner, see, e.g., Bierbooms (2005Bierbooms ( , 2006Bierbooms ( , 2007; Bierbooms and Cheng (2002); Bierbooms et al. (1999); Bos (2017); Larsen et al. (2003) and Nielsen et al. (2003). Here the gust is not generated by superposition of a deterministic mean shape to the velocity signal, but the velocity is constrained to a prescribed gust velocity at a pre-defined time, so that the whole signal fulfills the correlations of the isotropic turbulence theory. The resulting constrained stochastic gust model was compared with classical EOG and successfully applied to simulations for the flow around wind turbines.
Since the main objective of the present contribution is to describe a method how wind gusts can be imposed on the fluid field, this aspect is described in more detail in the following. The classical method denoted far-field boundary condition is to inject wind gusts as inlet boundary conditions. Since the gust superimposed at the inlet boundary induces a strong velocity perturbation, for an incompressible fluid it is necessary to dynamically adapt the velocity at the inlet boundary far away from the gust, so that the global inlet mass flow remains constant. Norris et al. (2010) proposed a spatial and temporal correction of the inlet velocity profile to fulfill this requirement in their LES carried out to investigate the three-dimensional flow and turbulence structures in the wake of a wind turbine impacted by an EOG. Later the same group studied the turbulent flow around three wind turbines introducing ECG at the inlet boundary (Storey et al., 2013(Storey et al., , 2014. Besides LES a variety of simulations relying on URANS exist requiring lower resolutions in space and time. For example, the flow around a vertical axis wind turbine undergoing an EOG was computed in 2D by Onol and Yesilyurt (2017). Various EOG were introduced at the inlet and different amplitudes were considered. Similarly, relying on 3D-URANS simulations, Menegozzo et al. (2018) studied the response of a horizontal axis wind turbine undergoing an EOG. As before the gust was introduced at the inlet boundary and freely propagated through the domain. The CFD simulation was coupled with the rotor dynamics, so that the effects induced by the gusts on important parameters of the problem such as the load coefficient were observable.
Besides the challenge of satisfying the global mass conservation mentioned above, the far-field boundary conditions possess two more disadvantages. First, the gust imposed at the inlet has to travel through the entire computational domain before reaching the region of interest, i.e., the building or the wind turbine. Since the grid in the far field is often rather coarse in order to save grid points for the region of interest, the effect of numerical damping of the gust during the approaching phase can be significant. This negative effect can be drastically reduced by a gusttransport mesh, i.e., a secondary fine mesh moving with the gust within the domain using the overset-grid technique (Heinrich, 2014), which is, however, rather costly and not straightforward to implement in general-purpose CFD codes. Second, the entire time interval of this prelude has to be simulated, which in case of a fine grid applied to reduce the damping effect can be very CPU-time consuming.
In order to avoid these drawbacks a low-cost approach called field velocity method (FVM) was proposed in the context of aeronautics. Note that wind gusts are also of interest for the aircraft design, since these events can trigger devastating phenomena such as flutter and stall (Radespiel et al., 2013). The concept of FVM was originally introduced into an unsteady Euler solver by Singh and Baeder (1997) and Parameswaran and Baeder (1997). The objective was to compute the step response of a wing to a sudden change in the angle of attack. The underlying idea of FVM is to alter the convective fluxes by superposition of a pre-defined additional velocity varying in space and time following a gust profile. Since this additional gust velocity is prescribed, the gust can be arranged in the computational domain such that it hits the structure at the first time step of the simulation, avoiding additional CPU costs for transporting the gust from the inlet to the zone of interest. However, the fluid flow and the structure cannot interact with the prescribed wind gust and alter it. Thus, the feedback effect is not taken into account. Despite this important drawback the field velocity method is widely used, see, e.g., Weish€ aupl and Laschka (2001), Ghoreyshi et al. (2018a, b), Heinrich (2014), F€ orster andBreitsamter (2015) and Kelleners and Heinrich (2016), where in the last three references it is denoted disturbance velocity approach (DVA).
URANS results obtained by FVM/DVA were compared with the results obtained by the injection of the gust at the inlet studying 2D (airfoils) and 3D geometries (realistic aircraft shapes) in Ghoreyshi et al. (2018a,b); Heinrich (2014); Kelleners and Heinrich (2016). For wind events with short durations deviations of integral quantities such as the lift force appear, whereas the method delivers similar results as long as the gust duration is long enough.
In order to improve the field velocity method (or DVA) Wales et al. (2014) rearranged the Euler equations similar to the field velocity method, but without the assumption that the airfoil has no feedback effect on the gust. The idea lead to the so-called split velocity method (SVM), which introduces additional gust related source terms into the Euler equations. These gust related source terms correspond to the effect of the body on the gust. The SVM was systematically compared with the FVM and with a simulation injecting the gust at the inlet boundary. The SVM (similar to FVM) offers the advantage of reduced CPU costs in comparison to the injection of the gusts at the inlet. However, compared to FVM the split velocity method leads to more reasonable results for gusts with shorter wavelengths.
The objective of the present study is to develop a methodology to efficiently inject strong (turbulent) wind gusts into the flow domain with the intention to be applied later on to fluid-structure interaction simulations. The intended methodology should generate individual wind gusts or a series of such events by an injection within the computational domain and close to the object of interest. Similar to the formulation used for injecting artificial turbulent fluctuations (Breuer, 2018;De Nayer et al., 2018b;Schmidt and Breuer, 2017;Wood et al., 2016) it is based on a source-term formulation and especially considers the feedback effect of the structure and the resulting turbulent flow on the wind gust itself, which is of major importance for highly turbulent flows and their prediction via LES.
Note that the present paper does not intend to predict realistic wind gusts typical for atmospheric conditions. That is beyond the scope of the contribution focusing on the injection methodology and thus considering solely the injection of gusts into either laminar or turbulent free flows. The application towards atmospheric boundary layer flows is intended for subsequent investigations.
The paper is organized as follows: Section 2 describes the general concept and the derivation of the source-term formulation. In Section 3 different shape functions for the wind gusts are presented and a slightly modified C 2 -1-cosine shape is proposed. The finite-volume flow solver which is used in the present study is briefly described in Section 4. The test cases can be found in Section 5 followed by the results and conclusion sections.
where u is the velocity vector and p the pressure. ρ denotes the fluid density, μ the dynamic viscosity and S the source term not depending on gusts, e.g., the gravity force. V is the volume of the cell. t denotes the physical time.
A new momentum source term S gust (with unit: kg m s À2 ) is added to the previous conservation equation to model the generation of a gust inside the computational domain. The gust induced by this source term is only driven by momentum. No additional mass is injected into the domain. Thus, the mass conservation equation remains unchanged and does not include an additional source term. Accordingly, Eq. (1) is simply extended to: First, this gust momentum source term is formulated in the local basis B 1 ¼ ðO; e ξ ; e η ; e ζ Þ. The direction of the gust is given by the first vector of the local basis, i.e., e ξ . This is a user-defined parameter. e η and e ζ are defined in such a manner that B 1 forms an orthonormal basis. Since convection in a gust plays the dominant role, the gust momentum source term is expressed analog to the convection term of Eq. (2), but with the total velocity of the gust denoted u gust instead of the velocity of the fluid without the gust. That means: The vector of the total gust velocity u gust is defined as follows: u base is the vector corresponding to the velocity of the surrounding fluid without the gust. Using the explicit time-marching scheme described in Section 4, this velocity is given by the local velocity of the previous time step u n . u gust deterministic is the vector, which only contains the velocity profile of the gust without the surrounding fluid. This gust shape is deterministically defined in the local basis B 1 ¼ ðO; e ξ ; e η ; e ζ Þ depicted in Fig. 1. Its definition is based on a user-defined amplitude of the gust A g and userdefined analytic functions representing the shape of the gust in time and space. As mentioned above, the direction of the gust is the first vector of the local basis e ξ . Thus, the velocity of the gust in the local basis reads: The typically used functions for the spatial f ξ ðξÞ, f η ðηÞ, f ζ ðζÞ as well as the temporal f t ðtÞ distributions related to a gust are described in Section 3.
The momentum source term (3) can now be simplified. Since the second and third components of u gust written in B 1 , u gust η j B 1 and u gust ζ j B 1 , are negligible compared with its first component u gust ξ j B 1 (the η and ζ-components are solely composed of the components of u base in B 1 ), only the first component of ðu gust j B 1 Á rÞ remains. Thus, Eq. (3) reads: Furthermore, it is assumed that the gradient in the gust direction of the velocity without gust ∂ u base j B 1 =∂ξ is negligible compared with ∂ u gust deterministic j B 1 =∂ξ. Accordingly, Eq. (6) can be approximated: As usual in a finite-volume discretization, the midpoint rule is applied to approximate the volume integral: Here the expression cc abbreviates the notation cell center. The vector e ξ present under the derivative is constant and can be extracted from this derivative. Combining it with u gust ξ j B1;cc , the velocity u gust j B 1 ;cc appears again and the expression reads: analog to a mass flow rate u gust j B1;cc : Using the Gauss theorem Eq. (9) can be interpreted as a mass flow rate in ξ-direction multiplied by the total gust velocity.
S gust j B 1 has now to be written in the Cartesian basis B 0 ¼ ðO; e x ; e y ; e z Þ. The vector e ξ is a user-defined parameter in B 0 . Since it defines the direction of the gust, it will be denoted n g in the following. The transfer matrix T B 0 →B 1 from basis B 0 to basis B 1 is known: where ξ x , ξ y and ξ z are the x À , yÀ and z À components of the vector e ξ . The same notations are used for the components of the vectors e η and e ζ . Since B 0 and B 1 are orthonormal bases, T B 0 →B 1 is an orthogonal matrix. Thus, the transfer matrix T B 1 →B 0 from basis B 1 to basis B 0 can be easily deduced as the transposed matrix of T B 0 →B 1 : The local coordinates ξ, η and ζ can be also expressed depending on the Cartesian coordinates: It leads to the following expression for the gust velocity: u gust deterministic j B 0 ðt; x; y; zÞ ¼ A g f t ðtÞ f ξ ðt ξ ðx; y; zÞÞ f η ðt η ðx; y; zÞÞ f ζ ðt ζ ðx; y; zÞÞ n g j B0 : Consequently, the source-term formulation in the Cartesian basis B 0 reads: To use the same terminology as in the source-term formulation of Pasquetti and Peres (2015), who injected a micro-jet inside a fluid domain, an expulsion phase and an ingestion phase have to be distinguished during the gust injection. The gust is only formed during the expulsion phase when the velocity rises. The ingestion phase does not beneficially contribute to the generation of the gust. Thus, in order to consider only the expulsion phase in the present gust source term, the idea has to be transformed from the temporal to the spatial coordinate. Taking one of the gust shapes defined below in Section 3 for f ξ , only the left part of these functions (based on the notation of Section 3 that means ξ < ξ g , where ξ g is the center of the gust injection in ξ-direction) leads to a positive mass flow rate as denoted in Eq. (9) and thus to a beneficial source term supporting the gust generation process. In summary, the final source-term formulation reads: Regarding this new source-term formulation the following features have to be clearly mentioned: Contrary to the micro-jet model by Pasquetti and Peres (2015) using a punctual source, the present gust model relies on a volumetric distribution of the source term. Since the underlying physics of micro-jets and gusts significantly differ, the model for the micro-jets relies on a momentum and an additional mass source term, whereas gusts can be generated without a mass source term. Although in lab-scale experiments gusts are often generated by a supply of an additional mass flux (see, e.g., Volpe et al., 2013;Tischmacher and Ruck, 2013;Gromke and Ruck, 2018), in nature that is not the case. Instead the existing fluid of the surrounding is locally accelerated. Consequently, the pure source-term formulation for the momentum equation as suggested here is the more natural choice in case of gusts. Owing to the fact that no additional mass flux is included in the source-term formulation, it is not necessary to relax the divergencefree constraint in the region of the injection as done in the case of the micro-jets (Pasquetti and Peres, 2015).

Gust shapes
In this section the typically used deterministic functions required by Eq. (14) for the spatial and temporal distributions related to a gust are described. Thus, the variable ϕ might be the coordinates (ϕ ¼ fξ;η;ζg) or the time t. The constant ϕ g is the central value of the gust distribution for the corresponding ϕ and L ϕ g its length or time scale. As mentioned in the introduction the Extreme Coherent Gust (ECG) and the Extreme Operating Gust (EOG) are defined in the IEC-Standard (2002). Another option is a Gaussian distribution.
Adapted "1-cosine" shape (Extreme Coherent Gust, ECG): The original "1-cosine" shape found in the literature (IEC-Standard, 2002) is expressed as follows: In order to achieve more control, the original shape is adapted introducing the central value ϕ g (De Nayer et al., 2019): Adapted "Mexican hat" shape (Extreme Operating Gust, EOG): In the same manner the original "Mexican hat" shape (IEC-Standard, 2002) is adapted introducing ϕ g (De Nayer et al., 2019):

Fig. 2.
Typical shapes for the deterministic gust models for a given scale L ϕ g and a given central value ϕ g . Gaussian distribution: This Gaussian distribution is similar to the distribution used in the digital filter method applied to generate synthetic turbulent velocity components (Klein et al., 2003). The coefficient 18 is chosen so that the value of the function is less than 0.01% of the peak value at ϕ ¼ À L ϕ g = 2 and ϕ ¼ L ϕ g =2. Fig. 2 illustrates the three typical distributions for a given scale L ϕ g and a given central value ϕ g .
Since the source-term formulation proposed in Eq. (14) also relies on the derivative ∂f ξ =∂ξ of the shape function, the three corresponding derivatives are given and plotted in Fig. 3.
Adapted "1-cosine" shape (Extreme Coherent Gust, ECG): Adapted "Mexican hat" shape (Extreme Operating Gust, EOG): Gaussian distribution: In the current source-term formulation S gust j B 0 has non-zero values only for ξ < ξ g (see Eq. (14)). Applying this cut on the above presented derivatives leads to a kink in the derivative of each function at ξ ¼ ξ g or more generally at ϕ ¼ ϕ g . In the CFD simulations this discontinuity may produce spurious numerical oscillations during the introduction of the gust. In order to resolve this problem, the function ∂f ϕ =∂ϕ should possess a continuous derivative at ϕ ¼ ϕ g . For this purpose, such a function with a continuous second derivative, denoted ECG-C 2 , is derived from the classic ECG: Fig. 3. Derivatives of the three shape functions for the deterministic gust models for a given scale L ϕ g and a given central value ϕ g . Fig. 4. Proposed C 2 -"1-cosine" gust shape and its derivative for a given scale L ϕ g and a given central value ϕ g . At first, the derivative of the requested function is defined based on the Extreme Coherent Gust distribution: After integration the function ECG-C 2 reads: Fig. 4 shows the ECG-C 2 distribution and its derivative.

Flow solver
As mentioned above the Navier-Stokes equations for an incompressible fluid are solved based on a finite-volume method, i.e., an enhanced version of FASTEST-3D (Breuer et al., 2012;Durst and Sch€ afer, 1996) and corresponding validation studies . The equations are discretized on a curvilinear, block-structured body-fitted grid with a collocated variable arrangement.
The surface and volume integrals are approximated by the midpoint rule. Most flow variables are linearly interpolated to the cell faces leading to a second-order accurate central scheme. The convective fluxes are approximated by the technique of flux blending (Ferziger and Peri c, 2002;Khosla and Rubin, 1974) to stabilize the simulation. For the current case the flux blending includes 3% of a first-order accurate upwind scheme and 97% of a second-order accurate central scheme. This choice is motivated as follows: Based on the standard grid defined in Section 5 and the reference simulation explained in Section 6.2.1 (Table 1) the local P eclet number Pe Δx ¼ ρ u Δx=μ predicted with the grid spacing Δx is much larger than 2. Accordingly, a pure central scheme is particularly susceptible to 2Δx oscillations. The flux blending applied avoids this problem. The subdivision (3%/97%) relies on extensive tests based on the reference case, which on the one hand have shown that below 3% upwind oscillations cannot be avoided. On the other hand, higher portions lead to an increase in numerical diffusion, which can be avoided. Furthermore, the momentum interpolation technique of Rhie and Chow (1983) is applied to couple the pressure and the velocity fields on non-staggered grids.
A semi-implicit predictor-corrector scheme (Breuer et al., 2012) (projection method of second-order accuracy) is used to solve the pressure-velocity coupling problem, which is appropriate for the large-eddy simulation (LES) technique. Here a low-storage multi-stage Runge-Kutta method is applied for time-marching the momentum equations within the predictor-corrector scheme. The subsequent corrector step ensures that mass conservation is achieved in form of a divergence-free velocity field. For this purpose, a Poisson equation for the pressure correction is solved by an incomplete LU decomposition method (Stone, 1968). Overall, that provides second-order accuracy in space and time.
The prediction of turbulent flows relies on LES, where the large scales of the turbulent flow are resolved directly and the small scales are modeled by a subgrid-scale (SGS) model. For this purpose, several models are implemented in FASTEST-3D such as the classical Smagorinsky model (Smagorinsky, 1963) combined with the Van-Driest damping function near rigid walls, the dynamic Smagorinsky model (Germano et al., 1991) and the WALE model (Nicoud and Ducros, 1999). Owing to minor influences of the SGS model considered for the present applications, solely the standard Smagorinsky model with the classical model parameter C s ¼ 0:1 is applied. For LES an appropriate representation of the inflow boundary conditions is often necessary (Wood et al., 2016;Breuer, 2018). In FASTEST-3D an inflow generator based on the digital filter concept of Klein et al. (2003) mimics incoming turbulent perturbations. The original concept was extended to permit the superposition of these turbulent perturbations as source terms inside the domain (Breuer, 2018;De Nayer et al., 2018b;Schmidt and Breuer, 2017).

Test cases
In order to test the proposed source-term formulation for injecting gusts, two test cases are chosen. First, a gust is injected into an undisturbed flow and its evolution in space and time is investigated. Based on this test case the first goal of this study can be evaluated in detail, i.e., the injection of a gust within the computational domain of a simulation. Second, a gust is inserted into a turbulent flow to address the second goal, i.e., to simulate gusts and their propagation within turbulent flows.
Note that the objective here is not to predict realistic wind gusts typical for atmospheric conditions as already motivated in the introduction. Thus, neither an atmospheric boundary layer flow is considered, nor specific characteristics of wind gusts under atmospheric conditions such as their coherence in the vertical direction (Branlard, 2009) or vertical profiles of the streamwise velocity across the gust events (Brasseur, 2001) are taken into account.
The geometry of the free flow considered is a cuboid (length Â height Â depth ¼ 4L Â L Â L) as sketched in Fig. 5. Without loss of generality the west face is assumed to be the inlet, whereas the east face is the outlet. The other faces are periodic in pairs. The grid is block-structured and contains on the standard level 512 Â 128 Â 128 cells, which are equidistant in all directions. In order to assure the quality of the results, a grid independence study is carried out in Section 6.1. For parallelization based on domain decomposition the grid is decomposed into 128 blocks and the simulations run on 128 processors.
Unless otherwise indicated, the following gust shape is assumed as the standard case: In ηand ζ-direction the ECG shape function is chosen.
Since the derivative of the ξ-function with respect to ξ appears in the final source-term formulation, the ECG is replaced by the ECG-C 2 shape in ξ-direction as already discussed in Section 3. Assuming that the ξ-direction and the time are linked by the Taylor hypothesis, the ECG-C 2 shape is also applied for the temporal distribution. Furthermore, the length in ξ-direction L ξ can be computed based on L t : where u conv is the convective velocity of the flow field. The direction of the gust is set to n g j B 0 ¼ e x and its center is located at Fig. 5. In order to carry out parameter studies, the time and length scales of the gust vary in the range 0.25-1. The different values investigated are listed in Tables 2  and 3 in Section 6.2. The background flow without the gust is also directed in x-direction. A constant inlet velocity is set to uj inlet ¼ u ∞ e x with u ∞ ¼ 1 m/s, which leads to a convective velocity of u conv ¼ u ∞ . The time step of the simulation is chosen to Δt ¼ 4 Â 10 À3 s, which guarantees a CFL number less than unity. The density of the fluid is set to ρ ¼ 1 kg m À3 and the dynamic viscosity to μ ¼ 10 À6 kg m À1 s À1 leading to a Reynolds number of Re ¼ 10 6 based on L.
For the turbulent test case the background flow is superimposed by homogeneous isotropic turbulence based on the methodology briefly described in Section 4. For this purpose, the digital filter concept of Klein et al. (2003) requires the following input parameters in order to generate turbulent perturbations. Solely the level of the Reynolds stresses and the definition of one integral time scale (T turb ) and two integral length scales (L turb y , L turb z ) are necessary to generate artificial turbulence with proper autocorrelations in time and two-point correlations in space. Owing to the simplifying assumption of isotropy of the approaching flow, which is the usual way for describing the background turbulence found in experimental facilities such as water or wind tunnels, the Reynolds shear stresses are set to zero and all three normal Reynolds stresses are equal The constant value u 0 u 0 chosen for the entire injection plane located upstream of the gust injection at about x=L % 0:06 (see Fig. 5) depends on the intended turbulence intensity TI ¼ ffiffiffiffiffiffiffi ffi u 0 u 0 p =u ∞ chosen. It is worth to note that in principle the synthetically generated turbulence can also be applied as inflow boundary conditions on the west face of the domain. However, here the source-term formulation for the turbulence injection is favored in order to demonstrate that it can be effectively combined with the source-term formulation of the gust. In the present study two different turbulence intensities are investigated, i.e., TI ¼ 4% and 8%. In dimensionless form the integral time scale is assumed to be given by T turb u ∞ =L ¼ 0:05. With the aforementioned assumption of isotropy and the Taylor hypothesis the two missing integral length scales are L turb y =L ¼ L turb z =L ¼ 0:05.

Grid independence study
In order to assure the quality of the results presented in the subsequent sections, a study on the effect of the grid resolution is carried out first. Starting from the standard grid (denoted grid S) with 512 Â 128 Â 128 cells two coarser and one finer resolutions are taken into account. For this purpose, the number of cells in all directions are halved or  Table 3 Parameter study with varying length scales L η vs. varying gust amplitudes A g (L t u ∞ =L ¼ 0:5, L ζ =L ¼ 0:25): Time-averaged maximal streamwise velocity. doubled, respectively. Thus, the finest grid F consists of 1024Â 256Â 256 cells, whereas the coarse grid C and the very coarse grid VC include 256Â 64 Â 64 and 128 Â 32 Â 32 cells, respectively. The time-step size is kept constant. The evaluation is done for a standard gust configuration, whose parameters are summarized in Table 1. This set of parameters is considered as the reference configuration in the whole study. Since an instantaneous flow field is considered, the evaluation is carried out at two specific instants in time, i.e., shortly after the injection process has finished (t u ∞ = L ¼ 0:8) and after the gust has traveled a certain distance through the domain (t u ∞ =L ¼ 1:6). Fig. 6 depicts the distribution of the streamwise velocity u=u ∞ and the pressure p=ðρ u 2 ∞ Þ along the streamwise coordinate at the center of the gust.
As visible in orange in Fig. 6(a) the evolution of the streamwise velocity and the pressure obtained on the very coarse grid VC differs significantly from the others. This difference increases in form and amplitude with time (see Fig. 6(b)). The coarse grid C delivers better results: The streamwise velocity and pressure distributions are similar to those predicted on the finer meshes. However, they quantitatively differ. Therefore, the grids VC and C are too coarse and will not be used in the following.
The results obtained on the standard grid S and the finest grid F are quasi identical in all areas of the curves just after the injection (see Fig. 6(a)). As time progresses a marginal increase of the deviations appears as shown in the close-ups in Fig. 6(b).
In conclusion, the grid independence study shows that the results converge with increasing grid resolution and that the deviations observed between the standard grid S and the finest grid F are marginal. Consequently, for the subsequent detailed investigations grid S is a good compromise between efficiency and accuracy. Fig. 6. Grid independence study: Dimensionless streamwise velocity and pressure after the injection of a gust into an undisturbed flow at two instants in time (Reference configuration of the gust: ECG-C 2 shape in ξ and t, standard ECG shape in ηand ζ-direction,

Reference simulation
In order to get a deeper insight into the physics of the flow during and after the injection of a gust, a simulation is carried out with the standard parameters summarized in Table 1. Fig. 7 depicts the evolution of the injected gust in time and space within the undisturbed flow. Since the direction of the gust is equal to the x-axis, the streamwise velocity and pressure along the central axis for y= L ¼ z=L ¼ 0 are investigated at certain instants in time. To complete the previous graphs, a 3D view of the wind gust is realized in Fig. 8. The values of the streamwise velocity and the pressure in the xy-plane at z= L ¼ 0 are considered and plotted as the third axis.
During the injection (t u ∞ =L 0:5, which corresponds to the first three instants depicted) the flow is accelerated up to u sim max =u ∞ ¼ 2:24. In this phase the x-distribution of the streamwise velocity is quasi symmetric with respect to x g =L ¼ 0:5. Moreover, comparing the velocity distribution of the third time instant (see Fig. 7(c)) with the theoretical shape of ECG-C 2 (see Fig. 4), the simulated shape looks very similar to the theoretical one. This is a particularly noteworthy result because it clearly Fig. 7. Dimensionless streamwise velocity (left) and pressure (right) of a gust injected into an undisturbed flow at different instants in time and for y= L ¼ 0 and z= L ¼ 0 (Reference configuration of the gust: ECG-C 2 shape in ξ and t, standard ECG shape in ηand ζ-direction, shows that the proposed source-term formulation works as desired. As a reaction to the additional momentum injected into the flow, the pressure takes a negative sinusoidal form. During the injection phase the pressure disturbance Δp=ðρ u 2 ∞ Þ reaches its maximum with about 0.18. In the rest of the simulation the dimensionless pressure varies below AE 0:1. After the injection phase (t u ∞ =L > 0:5) the gust travels through the space with a quasi constant propagation velocity equal to approximatively 1:4 u ∞ . This propagation velocity is evaluated based on the position of the velocity maximum as a function of time. It is larger than the velocity of the surrounding fluid, but much smaller than the maximum velocity of the gust. Indeed, the structure of the gust is composed of fluid elements, which move at highly different velocities: In the middle of the gust the magnitude of their velocities approaches u base þ maxðu gust deterministic Þ =u ∞ , i.e., 2 in the current case. At the periphery this magnitude is nearly equal to ju base j=u ∞ , i.e., 1 in the present case. The fast central region tries to accelerate the slower areas in the surrounding, which decelerates the faster parts. Hence, the resulting (average) propagation velocity is in between ju base j=u ∞ and u base þ During this propagation the shape of the wind gust changes. At first, the velocity distribution becomes asymmetric with respect to the Fig. 8. Three-dimensional view of the dimensionless streamwise velocity (left) and pressure (right) of a gust injected into an undisturbed flow at different instants in time in the xy-plane at z=L ¼ 0 (Reference configuration of the gust: ECG-C 2 shape in ξ and t, standard ECG shape in ηand ζ-direction, maximum velocity (see Fig. 7(d)). The fluid elements with the highest velocities in the gust structure flow faster than the slowest ones. A fast front forms and gradually detaches from the tail of the gust (see Fig. 7(f)) as observed in nature (Branlard, 2009). One important feature which Fig. 7 can not display is clearly visible in Fig. 8. If the gust develops along the x-axis leading to streamwise velocities, which are larger than the base flow, in the region offside the axis the streamwise velocity drops to values below u ∞ . That is obvious for example in Fig. 8(b), which depicts a trough in the streamwise velocity distribution next to the x-axis. Accordingly, the mass conservation is satisfied.
Similar to the streamwise velocity, the pressure distribution evolves in time: Starting from the sinusoidal form during the injection phase, it shows a kind of reversed Mexican hat shape (see Fig. 2(b)) when the gust travels through the domain.
The different physical observations concerning the propagation of the wind gust for the reference simulation are also valid for the simulations carried out in the following investigations.

Effect of the time scale L t on the effective length
During the generation of the gust three length scales and one time scale are required as shown in Section 2. However, because the Taylor hypothesis is assumed in ξ-direction, L ξ and L t are linked by relation (24). Therefore, the number of parameters reduces to three. L ξ , L η and L ζ are parameters of the functions f ξ , f η and f ζ , respectively. Theoretically, for a given dimensionless length scale L ξ =L, the length of the resulting injected gust has to be equal to L ξ =L and so on for the other two directions. In order to verify if the proposed source-term formulation delivers the correct effective lengths, simulations including a standard gust with an amplitude of A g =u ∞ ¼ 1 but with varying time scales L t u ∞ =L ¼ f0:25; 0:5; 1:0g are carried out (see Table 2). Fig. 9 compares the effective lengths of the injected gust in x-and ydirection with their associated theoretical values. The results are depicted along the x-axis with y=L ¼ z=L ¼ 0 for three different time scales. Furthermore, the distributions along the lateral axis are depicted. Since the results in z-direction are identical with those in y-direction, they are not shown here. Fig. 9. Dimensionless streamwise velocity of a gust injected into an undisturbed flow for varying time scales L t u ∞ =L at different instants in time: On the left just after the injection (t ¼ L t ), on the right at t u ∞ =L ¼ 2. The results are taken along the x-axis with y=L ¼ z=L ¼ 0. Additionally, the maximum of each streamwise gust velocity profile is selected and its associated streamwise velocity distribution along the lateral y-axis is depicted. (Configuration of the gust: ECG-C 2 shape in ξ and t, standard ECG shape in ηand ζ-direction, The following observations can be made: Just after the injection of the gust at t ¼ L t the shape of the gust in xdirection (except for the case L t u ∞ =L ¼ 0:25) is already no longer symmetric, but its effective length fits well to the associated theoretical value L ξ =L. Obviously, the asymmetry increases with increasing values of L ξ =L.
After the injection the gust propagates through the flow field. A fast front and a slower tail form as mentioned before. The length of the tail is linked to L ξ =L or L t u ∞ =L, respectively. A small value of L ξ = L leads to a short tail and vice versa. The effective length of the fast front is not strongly influenced by the size of L ξ =L or L t u ∞ =L, respectively. The width of the gust shape in y-direction agrees very well with the theoretical value L η =L directly after the injection and also later. It is not affected by varying the time scale L t u ∞ =L.
6.2.3. Effect of the length scale L η (or L ζ ) on the effective length To investigate the influence of the lateral length scale L η = L on the gust development, this input parameter is now varied within the range f0:25; 0:50; 1:0g (see Table 3). For this purpose, the time scale is fixed to the value L t u ∞ =L ¼ 0:5. Please note that due to the symmetry of the problem the same results can be obtained by varying L ζ instead of L η .
The results are depicted along the central x-axis (y= L ¼ z= L ¼ 0) in Fig. 10. Again, the distributions along the lateral axis are included in the graphs. The following observations can be made: Just after the injection of the gust at t ¼ L t the effective length of the gust shape in x-direction does not depend on the variation of the length scale L η =L. During the propagation phase the effective length of the gust shape in x-direction is affected by the variation of the initially prescribed length scale L η =L. By increasing L η =L, the separation between the fast front and the slower tail is attenuated, impacting the length of the gust. During the propagation phase the effective width of the gust in ydirection is no longer equal to its originally specified length scale L η = L. Indeed, the part with the fast velocity, which is located in the middle of the gust, detaches and moves faster than the rest. The sideward (left and right) edges of the gust are slower. Therefore, the effective width of the gust in y-direction is reduced with increasing time.

Effect of the time scale L t on the strength of the gust
To investigate the influence of the time scale L t on the strength of the generated gust, the maximum dimensionless gust velocity obtained at each time step is time-averaged over the time period between the end of the injection phase (t ¼ L t ) and the given time t u ∞ =L ¼ 2. The resulting value u sim max =u ∞ is compared with the theoretical value u theory max =u ∞ . Table 2 provides an overview of the results obtained for different dimensionless amplitudes of the gust A g =u ∞ varied in the range f0:5; 1:0; 2:0g. Additionally, the dimensionless time scale L t u ∞ =L varies between the values f0:25; 0:50; 1:0g. For this analysis the lengths of the gust in ηand ζ-direction are fixed to L η =L ¼ L ζ =L ¼ 0:25. According to the results summarized in Table 2 the gust velocities obtained for the shortest time scale L t u ∞ =L ¼ 0:25 are too low (about À11% to À15% error) and the ones obtained for the largest time scale L t u ∞ =L ¼ 1:0 are too large (error between þ4% and þ15%). However, the error strongly decreases for L t u ∞ =L ¼ 0:5. Based on these results it can be concluded that a dimensionless time scale approximatively two times larger than the dimensionless length scales L η =L and L ζ =L is required to achieve the best agreement with the intended maximum gust velocity.
6.2.5. Effect of the length scale L η (or L ζ ) on the strength of the gust For the investigation on the influence of the length scale L η =L in the range f0:25;0:50;1:0g, the time scale is fixed to L t u ∞ =L ¼ 0:5. As shown in Table 3 the error made on the time-averaged maximum gust velocity u sim max =u ∞ as defined in the previous subsection is low and between þ4% and À4% for the shortest length scales L η =L ¼ L ζ =L ¼ 0:25. With an increasing length scale L η , the relative error increases. Moreover, the increase of the gust amplitude A g =u ∞ also increases the relative error obtained.

Effect of the different shape functions
The final source-term formulation proposed in Eq. (14) relies on the four shape functions f t , f ξ , f η and f ζ and also on the derivative of f ξ with respect to ξ. Since these shape functions determine the spatial and temporal structure of the gust, the different classic shape functions presented in Section 3 are evaluated and compared with the reference case.
In the present method the ξ-direction and the time t are assumed to be linked by the Taylor hypothesis. Therefore, it is more consistent to use the same shape functions for f ξ and f t as done in the previous investigations. However, to distinguish the effects of each parameter in Eq. (14), the shape functions f ξ and f t will nevertheless be varied independently in this subsection.
First, solely the distribution in time f t is varied. The tested shapes are the ECG, EOG and the Gaussian function. The obtained results are compared in Fig. 11 with the reference simulation based on ECG-C 2 . The other shapes f ξ , f η and f ζ are identical to the reference case (f ξ ¼ECG-C 2 , f η ¼ ECG and f ζ ¼ ECG). The deviations between the velocity profiles in xand y-direction obtained for ECG and ECG-C 2 are minimal as expected. Applying the Gaussian distribution the form of the gust in x-direction remains similar, but its strength is slightly decreased. This observation Fig. 10. Dimensionless streamwise velocity of a gust injected into an undisturbed flow for varying length scales L η =L at different instants in time: On the left just after the injection (t ¼ L t ), on the right at t u ∞ =L ¼ 2. The results are taken along the x-axis with y=L ¼ z=L ¼ 0. Additionally, the maximum of each streamwise gust velocity profile is selected and its associated streamwise velocity distribution along the lateral y-axis is depicted. (Configuration of the gust: ECG-C 2 shape in ξ and t, De Nayer, M. Breuer Journal of Wind Engineering & Industrial Aerodynamics 207 (2020) 104405 can be explained by the fact that the integral value of the Gaussian function over the injection phase is lower than those of ECG or ECG-C 2 .
The use of the Gaussian function for f t does not affect the gust structure in y-direction (or z-direction). However, the results predicted with EOG used for f t differ in form and strength. This can be attributed to the special form of the Mexican hat function (EOG), in which a local minimum occurs before the maximum is reached. Due to the normalization of the Mexican hat function by the difference between the minimum and the maximum value, the maximum itself is also 26% smaller than for the other functions. Since the source term quadratically depends on the gust velocity, the effect on the resulting amplitude of the gust visible in Fig. 11 is even larger.
In the second series of predictions, f t is now fixed to ECG-C 2 as in the reference case and the distribution in ξ-direction f ξ is varied. The tested shapes are the ECG, ECG-C 2 and the Gaussian function. Note that the application of EOG in the main flow direction has been deliberately excluded. The reason for this decision is given by the fact that the special shape of the Mexican hat is lost very quickly during convective transport by the flow, since the fast fluid elements catch up with the slower part very quickly and thereby extinguish it. The shapes in the lateral Fig. 12. Dimensionless streamwise velocity of a gust injected into an undisturbed flow for varying spatial shape functions f ξ at different instants in time: On the left just after the injection (t ¼ L t ), on the right at t u ∞ =L ¼ 2. The results are taken along the x-axis with y=L ¼ z=L ¼ 0. Additionally, the maximum of each streamwise gust velocity profile is selected and its associated streamwise velocity distribution along the lateral y-axis is depicted. (Configuration of the gust: Standard ECG-C 2 shape in time, standard ECG shape in ηand ζ-direction, L t u ∞ =L ¼ 0: Fig. 13. Dimensionless streamwise velocity of a gust injected into an undisturbed flow for varying spatial shape functions f η at different instants in time: On the left just after the injection (t ¼ L t ), on the right at t u ∞ =L ¼ 2. The results are taken along the x-axis with y=L ¼ z=L ¼ 0. Additionally, the maximum of each streamwise gust velocity profile is selected and its associated streamwise velocity distribution along the lateral y-axis is depicted. (Configuration of the gust: Standard ECG-C 2 shape in ξ and t, standard ECG shape in ζ-direction, L t u ∞ =L ¼ 0:  11. Dimensionless streamwise velocity of a gust injected into an undisturbed flow for varying temporal shape functions f t at different instants in time: On the left just after the injection (t ¼ L t ), on the right at t u ∞ =L ¼ 2. The results are taken along the x-axis with y=L ¼ z=L ¼ 0. Additionally, the maximum of each streamwise gust velocity profile is selected and its associated streamwise velocity distribution along the lateral y-axis is depicted. (Configuration of the gust: Standard ECG-C 2 shape in ξ, standard ECG shape in ηand ζ-direction, L t u ∞ =L ¼ 0: directions f η and f ζ do not change (f η ¼ECG and f ζ ¼ ECG). Fig. 12 compares the results with the reference case. By varying f ξ only marginal changes are visible in the velocity distributions along the central axis. Some minor deviations are observed in the strength of the gust. Furthermore, a change in f ξ does not perceptibly affect the structure of the gust in y-or z-direction. Finally, the ECG-C 2 function is used for f t and f ξ and the shape function in η-direction is varied between ECG, EOG and Gaussian. Due to symmetry of the source-term formulation the same results would be obtained by changing f ζ instead of f η . Therefore, only f η is investigated.
The gusts predicted are depicted in Fig. 13. Since the ECG and the Gaussian distributions are comparable (see Fig. 2(a) and (c)), the results obtained for these two shapes are very similar along the y-direction.
Moreover, this change of the η-shape from ECG to Gaussian does not much affect the spatial distribution of the gust in x-direction and also in time.
Interestingly, when f t changes from ECG to Gauss, it has a noticeable impact on the strength. However, changing f η from ECG to Gauss only leads to marginal deviations during the injection phase and small difference after propagation. Again EOG delivers different results due to its special shape: In y-direction the central peak is surrounded by two along the x-axis at y=L ¼ z=L ¼ 0 (Configuration of the gust: ECG-C 2 shape in ξ and t, standard ECG shape in ηand ζ-direction, L t u ∞ =L ¼ 0: areas with a velocity magnitude below the background flow velocity. These two areas are more pronounced than with ECG or Gauss. According to the theoretical shape, the width of the central peak is smaller than with ECG or Gauss and its strength is also lower. This is also visible in x-direction. It is worth to note that the resulting form of the velocity distribution in x-direction after the injection is apparently not affected by EOG, only its maximum. Since the strength of the EOG gust is lower after the injection phase, its propagation in streamwise direction is slower than for the ECG or Gaussian gusts. Additionally, the two lateral areas of lower velocity obtained with EOG propagate slower than the center of the gust and are located at the end of the gust tail. The effective width of the EOG fast front is comparable with those predicted by ECG or Gauss. Only its strength differs.

Propagation of wind gusts in turbulent flow
A major objective of the present work is to demonstrate the possibility to generate wind gusts within turbulent flows. In the previous sections the Reynolds number is high, but the flow is undisturbed. In order to obtain turbulent flows synthetic turbulent perturbations with different turbulence intensities are added to the background flow of the following simulations as explained in Section 5. In these simulations the configuration of the wind gust is identical to the reference case in Section 6.2.1.
Before discussing the results the following should be noted about this test case: The presently used simplifications of generating homogeneous isotropic inflow turbulence (see Section 5) is solely an idealized case for the present study and does not claim to be a description of atmospheric turbulence. Of course, the flow in an atmospheric boundary layer requires inhomogeneous and anisotropic turbulent flow fields. However, two issues have to be taken into account here: First, in the present setup a free flow without the effect of any wall and thus boundary layer is considered with the single objective to demonstrate that the formulation of the source term for the gust and the background turbulence can be superimposed. Second, either with the presently used inflow generator or any other inflow generator (artificial or based on precursor simulations) the present restrictions to homogeneous isotropic turbulence can be easily overcome, which is the topic of ongoing work.
Fig. 14 depicts the evolution of the injected gust in time and space for the case with TI ¼ 4% and 8% along the x-axis at y=L ¼ z=L ¼ 0. To allow a comparison with the undisturbed case the streamwise velocity Fig. 15. Three-dimensional view of the dimensionless streamwise velocity of a gust injected into turbulent flows (left: TI ¼ 4%; right: TI ¼ 8%) at different instants in time in the xy-plane at z=L ¼ 0 (Reference configuration of the gust: ECG-C 2 shape in ξ and t, standard ECG shape in ηand ζ-direction, A g =u ∞ ¼ 1, L t u ∞ = L ¼ 0:5, L η = L ¼ L ζ =L ¼ 0:25). and the pressure are also plotted for the previous case of TI ¼ 0%. Similar to the undisturbed flow depicted in Fig. 8 (left) a 3D view of the streamwise velocity is added in Fig. 15. First, during the injection phase the wind gust shape is only marginally affected by the turbulence perturbations. This observation can be explained by the fact that the amplitude of the perturbations are small compared to the amplitude of the gust. Second, the propagation velocity of the gust is not influenced by the turbulent fluctuations, as visible in Fig. 14(c) and (d). However, the strength of the gust after the injection phase is reduced by the increase of the TI value. The turbulence present in the background augments the destruction of the gust structure. For example, the maximal streamwise velocity of the gust obtained with TI ¼ 8% at t u ∞ =L ¼ 2:0 is only 1:4 u ∞ and the one predicted within an undisturbed flow reached 1:93 u ∞ (see Fig. 14(d)).
Please note that, if the strength decreases with the increase of TI, the effective length in x-direction of the generated gust is not changing significantly. The effective length of the fast front and the tail remains similar. Moreover, no dispersion of the gust in y-and z-direction is noticeable due to the superimposed turbulent fluctuations (see Fig. 15(d)).

Conclusions
Based on the long-term goal of describing the effect of turbulent wind gusts on flexible structures in the context of coupled fluid-structure interaction simulations, an efficient method for introducing the gusts into the computational domain was developed in a first step. The objectives were to avoid CPU-time intensive techniques and to take the feedback effect of the flow field on the gust itself into account. Thus, a source-term formulation for the momentum conservation equation was derived, which allows to place the injection region of the gust close to the region of interest. Since no source term is added to the continuity equation, the local and global mass conservation is not touched and the divergence-free condition is still satisfied for incompressible fluids. Besides applying the classical shape functions for gusts (ECG, EOG) to define the source term, a modified C 2 -"1-cosine" shape ECG-C 2 was derived which guarantees that the source-term distribution has no kink. After a grid independence study a thorough analysis of the new gust injection technique was carried out in undisturbed and turbulent background flows. The main findings are: At the end of the injection phase the predicted shape of the wind gust is very similar to the intended one. Consequently, the proposed source-term formulation works as desired. When the source term is no longer active after the injection phase, the conservation equations of mass and momentum are fully satisfied. Consequently, the generated flow field including the wind gust are a valid solution of the Navier-Stokes equations.
The injected gust travels through the flow field with a nearly constant propagation velocity. In this phase, a fast and steeper front as well as a slower tail of the gust develop as observed in nature. The effective length of the fast front does not strongly depend on the size L ξ =L (or L t u ∞ =L, respectively), whereas the length of the slower tail is linked to this quantity. The prescribed length scale of the gust in the lateral directions significantly influences the separation between the fast front and the slower tail. An assessment of the different shape functions yields the finding that the most significant changes in the form and strength of the generated gust is achieved by the application of the EOG shape function, both in time and in η or ζ-direction.
A parameter study on the effect of the time scale L t has shown that the dimensionless time scale should be about two times larger than the dimensionless length scales L η =L and L ζ =L. The proposed source-term formulation for the wind gusts works in laminar and turbulent flows. It is not restricted to eddy-resolving simulation approaches (DNS or LES) and can also be combined with Reynolds-averaged Navier-Stokes solvers or hybrid LES-URANS techniques. The source-term formulation can be joined with other source-term approaches such as the turbulence generator used in the present study. The background turbulence does not perceptibly affect the propagation velocity of the gust. However, the strength of the gust is significantly attenuated by the turbulent fluctuations and this trend increases with raising turbulence intensities of the underlying flow field.
Next steps are the application of the gust injection technique in practically relevant flow cases (e.g., in atmospheric boundary layers with inhomogeneous and anisotropic turbulent flow fields) and later on the extension towards FSI.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
or in a user function as available in many commercial codes. Note that in contrast to other methods such as the field velocity method or the split velocity method no further modifications of the convective terms or the time integration are required for the present source-term formulation. Thus, its implementation is straightforward and can be done even in commercial codes.
Furthermore, the source-term formulation described above is written down here for the integral form of the Navier-Stokes equations, which is the basis for the finite-volume scheme applied in this study. If the solution method for the governing equations is based on the differential form of the equations (i.e., for a finite-difference scheme) the formulation in Eq. (14) has to be adjusted by canceling the volume of the control volume. Apart from that the source-term formulation is still valid and directly applicable.