Account of ambient turbulence for turbine wakes using a Synthetic-Eddy-Method

The present paper aims at describing the use of a Synthetic-Eddy-Method (SEM), initially proposed by Jarrin et al. [12], in the 3D Lagrangian Vortex method framework. The SEM method is used here in order to generate a far-field incoming flow with a prescribed ambient turbulence intensity. However, for the account of the diffusive term in the Navier-Stokes equations, a classical Particle Strength Exchange model with a LES eddy viscosity is used. Firstly, the general characteristics of the Synthetic-Eddy-Method will be presented together with its integration in the framework of the developed 3D unsteady Lagrangian Vortex software [27]. The capability of the ambient turbulence model to reproduce a perturbed flow that verifies any turbulence intensity I∞ and any anisotropic ratio (σu:σv:σw) will be discussed and validated. Then, the capability of the presented ambient turbulence model to compute turbine wakes will also be presented together with first results. Finally, comparisons will be made between the obtained numerical results against experimental data [22, 23] for two levels of ambient turbulence, namely I∞ = 3% and I∞ = 15%. Although the present study was initially performed in the framework of tidal energy, its application to wind energy is straightforward.


Introduction
The development of renewable energy has been rather intensive for some years, mainly using wind turbines. But tidal turbine are also arising. However, at the scientific level, many studies are still necessary in order to fully understand the behaviour of such arrays. One of these topics is the impact that ambient turbulence has on the behaviour of turbines in a farm. Indeed some studies demonstrated that ambient turbulence highly modifies the behaviour of turbines [26,1,6,5,4,8].
In the previous equation, u = (u, v, w) is the 3D velocity vector,ū stands for the averagevalue of any velocity record u and σ(u) represents the corresponding standard deviation. One of the most noticeable influence can be observed in the wake, downstream of the turbine. Indeed as shown in Ref. [22] for instance, the higher the turbulence intensity I ∞ is, the faster the wake effects decrease. Moreover the interactions between turbines are also highly modified by ambient turbulence as shown by Mycek et al. [23] in their study of the performances and wake of two interacting turbines. The relative spread in in situ ambient turbulence levels inhibits the definition of general guidelines on tidal arrays with respect to turbulence. For wind energy, although the fact that atmospheric boundary layer has a different behaviour than in sea, the situation may be similar. Therefore, the influence of each level of ambient turbulence needs to be characterised, either experimentally or numerically. Consequently numerical simulations have to represent the ambient turbulence, or at least its effects on the performance and wake of the turbines. In the present paper, the objective is not to describe a turbulence model (like a RANS closure, LES or DES models) but how to generate accurate turbulent inflow conditions, that represent a given ambient turbulence. In the literature, some numerical works were carried out in that aim. One of the most famous initialisation of ambient turbulence boundary conditions is TurbSim [13] from NREL. But many others exist. For example, Chatelain et al. [6] used the Mann's algorithm [17] to synthesise a turbulent inflow in the case of wind turbines simulations. Branlard et al. [5,4] used a similar Vortex method approach with frozen and un-frozen approach for turbulence and possibly shear. The present study considers the Synthetic-Eddy-Method (SEM), initially proposed by Jarrin et al. [12,11], because an easier adaptation to the Vortex Method framework is foreseen. Togneri et al. [30,31] already investigated a similar Synthetic-Eddy-Method in order to generate synthetic turbulent inflow conditions for their BEMT-CFD software.
This paper presents the latest numerical developments carried out to integrate a SEM module in the developed 3D Vortex Method software [27]. This software already possesses an LES turbulent model, its turbulent eddy viscosity being based on the work of Mansour et al. [18].
The first part of this paper is devoted to the presentation of the numerical method used to compute the flow around a wind or tidal turbine. Then in a second part, the Synthetic-Eddy-Method is presented as well as its integration in the present software. The last part of this paper focuses on preliminary results obtained using this method to represent the ambient turbulence I ∞ . Eventually, conclusions are drawn on the presented work and an outlook on future development is proposed. The blade is supposed infinitely thin and is presented as cutsection.

Numerical method
The developed software to simulate turbines and wakes is based on the Vortex method [29,14,9]. This method is an unsteady Lagrangian method, where the flow is discretised using vorticity carrying particles (see Fig. 1) and the turbines are represented using a panel method [3]. In this numerical method, the Navier-Stokes equations for an unsteady and incompressible flow are used in their velocity/vorticity (u, ω) formulation: One can refer to [27] for more detailed information of the implementation. But the real novelty of the present paper is the use of the Synthetic Eddy Method (SEM) to generate far-field ambient turbulence as it will be discussed in Section 3. The Helmholtz decomposition of the velocity field (eq. (4)) is used in order to decompose the velocity vector: Thus the velocity field is divided into three parts: • a potential component u φ representing the influence of the turbines blades on the flow (see Fig. 1 and 2), a panel method (boundary integrals) is used to solve the Laplace equation for φ. • a rotational component u ψ representing the wake of the turbines (see Fig. 1

and 2), a 3D
Vortex method is used to solve a Poisson equation for ψ. • a velocity vector u ∞ representing the upstream tidal current.  Here, the ambient turbulence intensity is set to I ∞ = 0%. The image is reproduced from Mycek et al. [24].
In order to make the link between the rotational and the potential parts of the velocity, particles are emitted at the blade's trailing edge using a Kutta condition [14]. The emitted particles are then advected in the flow as shown in Fig. 2. The rotational part of the velocity field is calculated as the sum of the influence of these vorticity carrying particles. The expression of the component u ψ is obtained in every point of the fluid domain by solving the Biot-Savart law [32]. The calculation of this u ψ -component for each particles represents the major part of the computational cost. For that reason, a Treecode algorithm inspired from Lindsay and Krasny [15] is used to reduce the computational time. Then, the results are post-treated in order to obtain a map of the velocity field as presented in Fig. 3 is reproduced from Mycek et al. [24], where a judicious pre-conditioner is presented for solving the potential component u φ in case of several turbines in interaction. In that respect, a preconditioned Bi-CGSTAB solver is presented in Ref. [24] that really speeds up the computation of this velocity component. For more details on the computations, including Biot-Savart equation desingularisation, remeshing or time-stepping, etc., one can refer to [27].

Synthetic-Eddy-Method
The Synthetic-Eddy-Method was initially developed to generate an inflow velocity boundary condition for a given turbulence intensity I ∞ and a given anisotropic ratio (σ u :σ v :σ w ) by adding a perturbation term, which is added to the far-field mean velocity of the flow The present paper aims at describing its implementation in the Lagrangian Vortex framework, and not only as inflow velocity condition but applied to the whole flow domain.

General overview of the method
Similarly to a Reynolds decomposition of the velocity, a fluctuating term u σ is added to the upstream mean velocity u ∞ and the upstream velocity u ∞ then becomes: On the contrary to the already presented method [12,11], this perturbation term is not only added at the inlet of the fluid domain but everywhere in the studied space. This aspect will be treated later in sub-sec. 3.3. The perturbation term u σ is calculated as the influence of N generated "turbulent structures", also called "eddies" in [12,11], randomly placed in the studied space. Then the perturbation term u σ induced by those N turbulent structures is simply defined as the sum of the influences of each structure k: with x being any point of the flow and u σ,k the velocity induced by a single turbulent structure k. This last term can be expressed as: where V σ is the volume of the three-dimensional space containing all the N "turbulent structures". At this point, the turbulent structure k and more particularly its intensity c k needs to be defined together with the function F λ .

Intensity and shape function
Each turbulent structure k is defined by its position x k in the flow and its intensity c k . The intensity c k is defined as: The different k i,j terms that intervene in equation (8) are independent random values that follow a normal or Gaussian distribution centred around 0 with a variance of 1. This k i,j term is supposed to represent the random aspect of turbulence. The elements a i,j represents any  (9)).
With these equations (8) and (9), the link between the intensities c k of the turbulent structures and the Reynolds Stress Tensor R ensures the generation of a velocity field that statistically replicates any given turbulence intensity I ∞ and any given anisotropic ratio (σ u :σ v :σ w ) [16,12]. Indeed the three components of the anisotropic ratio (σ u :σ v :σ w ) are just the square root of the diagonal components over the component R 11 of the Reynolds Stress Tensor R. Moreover the turbulence intensity I ∞ expression can be rewritten in function of the trace of the Reynolds Stress Tensor R as shown in equation (10): This last equation (10) ensures that the generated turbulence will have the desired turbulence intensity I ∞ and the desired anisotropic ratio (σ u :σ v :σ w ). The second part of equation (7) is F λ also called the shape function and it can be expressed as: λ defines the size of the zone of influence of each turbulent structure k. This size is probably linked to a turbulent length scale (Taylor or Kolmogorov length scales) but this will need further investigation. As it is defined by a vector λ, and presented in equation (11), the zone of influence can ideally have different sizes regarding each coordinate (λ i ). And additionally, each turbulent structure could have its own size vector λ but this will be treated in the following. The subfunctions f λ that intervene in the evaluation of the shape function F λ have to meet some conditions to ensure that the ambient turbulence has the chosen characteristics: argmax y (f λ (y)) = 0, (12a) In the preliminary results presented in section 4, the basic shape function F λ (eq. (11)) indicated in Jarrin et al. [12,11] is used. Indeed its sub-function f λ simply is a tent function, represented by a triangular function centred at zero with a base of 2λ. Figure 4 depicts tent functions for different values of λ. It is anticipated that this basic tent function (eq. 13) will not meet all the requirement of smoothness, especially if ones want to computed the velocity gradient, etc. However, in order to start with, this basic function issued from the bibliography will be used but other shape function were already tested and validated. The tent function for different sizes λ of turbulent structure. Figure 5. General representation scheme of the developed ambient turbulence model.

Integration of the Synthetic-Eddy-Method into the Vortex method
The numerical method presented in Section 2 is an unsteady Lagrangian method defined for an unbounded domain. Indeed, the fluid domain is supposed to be infinite but the zone where particles are present, grows with the particles emission and their advection in the flow. In fact, there is no domain boundaries as in Eulerian formulation. Unfortunately the Synthetic-Eddy-Method need the definition of a precise volume V σ (see eq. (7)) containing all the turbulent structures, which is not required for classical Vortex method. Thus in order to combine the Synthetic-Eddy-Method with the Vortex method, a studied space E S needed to be defined. The desired velocity fluctuations induced by the Synthetic-Eddy-Method will be applied inside this E S domain. Once the studied space E S is defined, a turbulent space E σ of volume V σ can be created. This turbulent space E σ contains all the turbulent structures used in the Synthetic-Eddy-Method. If one wants to have a statistically correct velocity fluctuations everywhere inside the studied space E S , some conditions on the turbulent space E σ are needed, especially at the limit. Therefore, the turbulent space E σ needs to be a little larger than the studied space E S : with x a position vector and λ = λ i , (i = 1, 3) the size vector of the turbulent structures. Those conditions simply ensure that all the study space is influenced by the turbulent structures. Now that the turbulent space E σ is defined, the generation of the turbulent structures can be performed. E σ is initialised with the N turbulent structures k and their intensities c k are calculated using equation (8) with k i,j taking the value 1 with a probability of 1/2 and the value −1 with an equal probability. Then each turbulent structure is randomly placed in the turbulent space E σ .
During the simulation, the turbulent structures will be advected by the flow. Here, different strategies can be applied: either the structures can be advected by u ∞ , which is somehow too simplistic, or the structures can be advected by u as defined in equation (4). Additionally, these structures could also be used in order to compute the stretching term of the Navier-Stokes equation (first term of the right end side of equation (3) but not detailed here). In this preliminary study, the first option was chosen that is: the structures are advected by u ∞ without influencing the stretching term. These structures work somehow as an under-layer, which means that they will influence the velocity distribution on the turbines (and hence the particles emission process) and also the particles in the wake. But they will not be influenced reciprocally.  Regarding the second option (the structures are advected by the actual fluid velocity u), several questions arise. In fact, the validity of the method relies on the randomness of the turbulent structures locations; will this property still be valid when advected by u? Thus, with the second approach, the Synthetic-Eddy-Method may not replicate the prescribed Reynolds stress tensor. Some more work needs to be performed on that prior to choose this second option.
Finally during the advection of the turbulent structures, each structure leaving the turbulent space E σ will be suppressed and replaced by a newly generated one at the inlet of the turbulent space E σ .

Numerical results
The numerical results presented in this section are issued from preliminary computations and will only be used as proof of concept. At first in sub-section 4.1, the ambient turbulence module will be run without any turbine to ensure that the model correctly reproduces the different input parameters. Then, computations will be presented in sub-section 4.2 with the presence of a single turbine for a moderate and a higher turbulence intensity I ∞ levels.

Results on the ambient turbulence model without turbine
The ambient turbulence module was firstly studied without turbine to ensure that the model correctly reproduces the different input parameters such as the turbulence intensity I ∞ and the anisotropic ratio (σ u :σ v :σ w ). This verification is quite straight forward as I ∞ , σ u , σ v and σ w can be deduced from the Reynolds Stress Tensor R (see eq. (10)).  σ(λ) = 0% σ(λ) = 25% σ(λ) = 50% σ(λ) = 75% σ(λ) = 100% Figure 7. Influence of enabling a variation of the turbulent structure sizes (represented by σ(λ)) on the power spectral density with a Gaussian function, a mean λ = (0.5, 0.5, 0.5), a filling ratio of R f = 2 and a given turbulence intensity of I ∞ .
To this aim, 3D turbulent velocity fields were generated on a fluid domain E S of dimensionless size 6×6×6. In order to validate the model, the generated flow-field characteristics are analysed. Table 1 gives examples of comparison between the Reynolds Stress Tensor R prescribed to the model and the recalculated ones from the generated flow. These examples are for diagonal Reynolds Stress Tensors R, which correspond to an ambient turbulent I ∞ of 15%. But the validation was also performed for other intensity values. The anisotropic ratio was arbitrarily set to the one observed in Milne et al. [21]. The recomputed values from the generated flow show an excellent agreement with the prescribed ones and thus regardless of the number of used turbulent structures, their size, functions, etc. This good agreement between the prescribed and recomputed values of the Reynolds Stress Tensor R proves the model capacity to generate velocity fluctuations replicating a given turbulence intensity I ∞ and anisotropic ratio (σ u :σ v :σ w ). Figure 6 displays the results obtained with the Synthetic-Eddy-Method for I ∞ = 15%. It shows an example of "turbulent" flow-field generated with the model for the same anisotropic ratio (σ u :σ v :σ w )=(1.0,0.75,0.56). One can observe that the generated flow seems to be very realistic, although the SEM model only is a mathematical tool that produces velocity fluctuations, without solving the fluid dynamics equations. Additionally, even though the turbulent structures sizes were equally set to λ = (1.0, 1.0, 1.0), this length cannot be clearly observed. At the present time, work is under progress to identify what is the link between this size vector λ and the turbulent length scales of the generated flow. This vector λ revealed to be an important parameter for this Synthetic-Eddy-Method and its influence needs to be further investigated. Figure 7 displays the power spectral density function obtained with the use of a Gaussian shape function, with a mean value for λ = (0.5, 0.5, 0.5) and increasing the spread in the λ values, characterised by σ(λ). σ(λ) represents the standard deviation in the size of the turbulent structures sizes λ. For instance, a σ(λ) = 50% for a mean λ = (1.0, 1.0, 1.0) means that each λ i value may vary between 0 and 2 and that the spread is characterised by a standard deviation of +/ − 0.5 around 1.0. The differences between those curves are quite striking. Indeed, when a single size of structure is used, the obtained power spectral density is spiked for the higher frequencies. This spiking behaviour of single sized-structures was analysed theoretically and Further to the previous study of turbulent flow without any turbine, few computations were run with a turbine. These runs were performed for the two turbulent intensities I ∞ of 3% and 15% , namely the lowest and highest turbulent intensity values for which experimental results are accessible in the IFREMER wave and current flume tank. Numerical parameters were those from Ref. [27] and the computation were scaled on the Reynolds number of the experiments, that is to say Re = u ∞ R ν = 0.28 × 10 6 . The obtained results are very encouraging and a qualitative analysis can already be performed. Figures 8 and 9 depict the numerically computed wakes for I ∞ = 3% and 15% and compared to the experiments performed by Mycek et al. [22]. The numerical velocity map presented on Figure 8 seems to indicate that the velocity deficit is underestimated for an ambient turbulence set to I ∞ = 3%. For the turbulent intensity of I ∞ = 15% an over-estimation seems to be observed. This discrepancies could partially be attributed to the fact that, since now, only rather coarse computations were tested. With finer discretisations, as those already used in [27], more accurate velocity deficits are expected. As already mentioned, these are only preliminary results, which need to be further validated with more computations (finer discretisations, a mesh independence analysis, other parameters validation, etc.). Hopefully, in a near future, more intense turbulent conditions are scheduled (several turbulence intensities and other anisotropic cases) in order to cover a larger range of ambient turbulence intensities.

Conclusions and Further works
The important role played by the ambient turbulence intensity on wind or tidal turbines' behaviour was addressed using a numerical investigation. A new module based on this Synthetic-Eddy-Method initially proposed by Jarrin et al. [12,11] was implemented in the numerical code in order to represent the ambient turbulence around turbines. The general characteristics of the Synthetic-Eddy-Method were presented together with its integration in the framework of a 3D unsteady Lagrangian Vortex method. The capability of the ambient turbulence model to reproduce a perturbed flow that verifies any turbulence intensity I ∞ and any anisotropic ratio (σ u :σ v :σ w ) was validated. Moreover the compatibility of the presented ambient turbulence model with the use of a turbine on the software was proven and first results in terms of turbine wake were presented. As a matter of short-term perspective, a better qualification of the Synthetic-Eddy-Method is needed in order to better understand the input parameters of the model. For instance, the role of the shape function need to be identified and improvements are already planned.
Another development on the ambient turbulence model could be the modification of the Synthetic-Eddy-Method in order to generate a divergence-free flow. In fact, divergence-free flow is a compulsory hypothesis inherent to the Vortex method in order to apply the Biot & Savart law. Unfortunately, the present implementation is for sure not divergence-free. For low to moderate turbulence intensities, this should not introduce a tremendous error in the flow. But increasing the ambient turbulence intensity might become an issue. Such kind of development were already initiated by Poletto et al. [28] but need a lot of work in order to be introduced in the framework of Vortex method.