# Studying the Wake of a Tidal Turbine with an IBM-LBM Approach Using Realistic Inflow Conditions

Mickael Grondeau, Sylvain Guillou, Jean Charles Poirier, Philippe Mercier, Emmanuel Poizot, Yann Méar

► To cite this version:

Mickael Grondeau, Sylvain Guillou, Jean Charles Poirier, Philippe Mercier, Emmanuel Poizot, et al. Studying the Wake of a Tidal Turbine with an IBM-LBM Approach Using Realistic Inflow Conditions. Energies, 2022, 15 (6), pp.2092. 10.3390/en15062092. hal-03614029

HAL Id: hal-03614029

https://normandie-univ.hal.science/hal-03614029

Submitted on 28 Mar 2022

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
Article

# Studying the Wake of a Tidal Turbine with an IBM-LBM Approach Using Realistic Inflow Conditions

Mickael Grondeau<sup>1,2 ID</sup>, Sylvain S. Guillou<sup>1,* ID</sup>, Jean Charles Poirier<sup>3 ID</sup>, Philippe Mercier<sup>1,* ID</sup>, Emmnuel Poizot<sup>1,4 ID</sup> and Yann Méar<sup>1,4 ID</sup>

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

check for
updates

**Citation:** Grondeau, M.; Guillou, S.S.; Poirier, J.C.; Mercier, P.; Poizot, E.; Méar, Y. Studying the Wake of a Tidal Turbine with an IBM-LBM Approach Using Realistic Inflow Conditions.

*Energies* **2022**, *15*, 2092. https://doi.org/10.3390/en15062092

Academic Editor: Antonio Crespo

Received: 24 February 2022

Accepted: 9 March 2022

Published: 13 March 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

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

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

**Abstract:** The lattice Boltzmann method is used to model a horizontal axis tidal turbine. Because tidal turbines generally operate in highly turbulent flows, a synthetic eddy method is implemented to generate realistic turbulent inflow condition. The approach makes use of the open-source code Palabos. Large eddy simulation is employed. A coupling between an immersed boundary method and a wall model is realized to model the turbine. Calculations are performed at two different turbulence rates. The upstream flow condition is first set up to match with experimental results. Numerical simulations of a tidal turbine with realistic turbulent inflow conditions are then realized with the lattice Boltzmann method. The approach is found to be in good agreement with experimental data. Cases with three different inflow turbulence rates are simulated. An almost linear evolution with the turbulence rate is observed for the axial velocity deficit. An analysis of the propagation of tip-vortices in the close wake is carried out. It is found that turbulence has a great impact on the tip-vortices propagation envelope.

**Keywords:** tidal turbine; wake; turbulence; numerical simulation; lattice Boltzmann method

## 1. Introduction

Worldwide energy consumption keeps increasing and some areas still do not have access to electricity. In order to provide more clean energy, new concepts have been created by academia and industry. Among those concepts are renewable marine energy devices and more precisely tidal turbines. Resulting from gravitational forces, tidal forces create a periodic movement of ocean water level [1]. Induced displacements are not constant and can highly vary from one place to another. In some sites, they create a strong current, like in the Alderney Race where it can go up to 5 ms⁻<fnref n="1" />. Using lift generating devices such as airfoils, tidal turbines extract the kinetic power from those currents. Interactions between tidal flows and their environments (sea bottom, coast, wind) may result in the creation of turbulent structures in the flow [2]. Several in site campaigns have already been performed to evaluate the turbulence rate of tidal sites. Thomson et al. [3] measured a turbulence rate between 8.4% and 11.4% in Puget Sound, USA. Turbulence can also vary along the tidal cycle of a specific site. Milne et al. [4] observed that the stream-wise turbulence rate varies between 7% and more than 20% along a tidal cycle at the sound of Islay, Scotland. The project Turbulence In Marine Environment (TIME) [2] lists several problems associated with marine turbulence. It quotes that turbulence may affect the turbines' integrity, performance and hydrodynamics. To understand turbulence's influence on tidal turbines, experiments in flume tanks have been carried out. Chen et al. [5] realized an experimental study of the

<footnote n="1">¹ Normadie University, UNICAEN, LUSAC, EA4253, 60, Rue Max-Pol Fouchet, CS 20082 Cherbourg, France; m.p.grondeau@soton.ac.uk (M.G.); emmanuel-poizot@lecnam.net (E.P.); yann.mear@lecnam.net (Y.M.)</footnote>

<footnote n="2">² University of Southampton, Boldrewood Innovation Campus, Southampton SO16 7QF, UK</footnote>

<footnote n="3">³ SIREHNA, Naval Group, Technocampus Océan, 5 rue de l'Halbrane, 44340 Bouguenais, France; jean-charles.poirier@sirehna.com</footnote>

<footnote n="4">⁴ Conservatoire National des Arts et Métiers-INTECHMER, Bd de Collignon, 50110 Tourlaville, France</footnote>

<footnote>* Correspondence: sylvain.guillou@unicaen.fr (S.S.G.); philippe.mercier@unicaen.fr (P.M.); Tel.: +33-2-33-01-40-32 (S.S.G.)</footnote>
wake, which is critical for the placement of tidal turbines in a farm. Mycek et al. [6] showed that turbulence greatly affects wake recovery and turbine performances. Blackmore et al. [7] showed that high turbulence rates can generate thrust fluctuation up to 20%, increasing fatigue failure risk. A good understanding of the turbulence influence on a tidal turbine is thus necessary in order to go from prototype to full industrial exploitation. In addition to flume experiments and in situ measurements, computational fluid dynamics (CFD) has also been used to gather such knowledge.

CFD simulations can be divided into two categories, steady simulations and unsteady simulations. Steady models have already proven their capability of predicting the performance of tidal turbines [8]. However, turbulence is chaotic and steady models may miss some important features of turbulent flows. Unsteady simulations can capture those phenomena. Several unsteady approaches have already been carried out. Elie et al. [9] and Grondeau et al. [10] used actuator large eddy simulation (LES) approaches to model tidal turbines. Their approaches do not model the geometry of the blade but rather represent it with a force. Forces coefficients, medium and far wake are correctly computed. Those kinds of methods cost less computational resources than geometry resolved LES and are well suited for far wake propagation and tidal farm studies. However, a priori knowledge of the blade behaviour is required. Geometry resolved LES simulations of tidal turbines have already been performed by Ouro et al. [11] and Ebdon et al. [12]. It was proved in the report that both turbine performances and wake are recovered. Ahmed et al. [13] used a Navier-Stokes LES model to predict the effect of turbulence on the loads applied to a real-size tidal turbine. Ebdon et al. [12] used a DES model to show that the length scale of turbulent structures might have a significant influence on the wake of the tidal turbine. All the previously mentioned studies used Navier-Stokes based models. The lattice Boltzmann method (LBM) is based on the Boltzmann equation. It is an unsteady method and derived from lattice gas automata [14]. LBM is currently not very popular in tidal turbine modelling. However, it is well suited to the detailed description of large areas thanks to its local and explicit formulation [14]. Moreover, it is known to be a low dissipative method [15], which is fundamental in modelling wake propagation. Setting up a reliable LBM model could therefore be of primary importance for tidal farms optimization. Some studies already exist in wind turbine modelling. Deiterding et al. [16] have correctly predicted the performance and wake of a wind turbine with a grid-adaptive LBM code. The computed tip vortices along with their interactions with the mast show that this method is well suited for complex wake modelling. The method was used to obtain a characterization of the turbulence in tidal site [17,18].

The LBM code used for this study works on a non-moving Cartesian mesh. A first comparison of the method with an NS-LES code has been made in [19]. The immersed boundary method (IBM) allows the modelling of complex geometries without any unstructured mesh fitting the geometry. The IBM is used here for modelling both the rotor and stator of the tidal turbine model. Ouro et al. [11] showed the capability of the IBM-NS to model the complex wake of vertical axis tidal turbine. Ouro et al. [20] also used an IBM-NS model to predict the influence of turbulence on the loading of a tidal turbine. The use of a Cartesian mesh can be disadvantageous since it increases the number of nodes in the mesh. A multi-level Cartesian mesh is used [21]. This allows having smaller grid spacing in areas of interest such as in the wake or close to the turbine, thus reducing the total number of nodes.

Turbulence intensity can vary from one site to another but it also varies along the tidal cycle of a specific site. Milne et al. [4] measured that the stream-wise turbulence rate varies between 7% and more than 20% along a cycle at the sound of Islay, Scotland. As the turbulence intensity varies significantly on a tidal site, it is important to know how a turbine reacts to several turbulence intensities. The influence of different turbulence intensities over quantities such as wake velocity and turbulence rate is investigated here. Upstream turbulence is generated with a divergence-free synthetic eddy method (DFSEM), introduced by Poletto et al. [22]. The work presented in this paper aims at developing
an LBM-LES tool to model tidal turbines in turbulent flows. This tool must be cheap enough to model several inflow configurations and turbine designs. This restriction is taken into account during the model setup. The case examined is a horizontal axis tidal turbine (HATT) model that was tested in Ifremer's flume tank [6]. This case provides all the data required to validate the simulations at two turbulent rates, 3% and 15%. Once the model is validated, a turbulent rate of 8% is simulated. Simulations are performed with the open-source code Palabos.

Section 2 describes the numerical and experimental set-up. Section 3 presents comparison of the model results with the experimental data from Mycek et al. In Section 4, the influence of the turbulence rate on the tidal turbine wake is investigated. Finally, Section 5 gives the conclusion and prospects of this study.

## 2. Materials and Methods

In this section, numerical methods and reference experiments are presented. The first part introduces the LBM and models chosen for this study. The experiment is then presented. The boundary conditions used to model the tidal turbine and the boundary conditions which generates inflow turbulence are presented in the third and fourth part. Finally, the numerical setup used to model the turbine is presented.

### 2.1. Lattice Boltzmann Method

The LBM is a numerical method that resolve the Boltzmann equation, Equation (1). It is an instationary weakly compressible method based on a statistical representation of the fluid through a density distribution function $f(x, \xi, t)$. This last represents the probability density for molecules being at a location $x$, which have the velocity $\xi$, at time $t$. The Boltzmann equation reads:

$$ \frac{\partial f(x, \xi, t)}{\partial t} + \xi \cdot \nabla f(x, \xi, t) = \Omega(f(x, \xi, t), f(x, \xi, t)), \quad (1) $$

where $\Omega$ is the collision operator and represents the collision between molecules. In the LBM, velocities $\xi$ are discretized using a predefined set of directions of propagation $c_i$ for molecules. A 2D example is given in Figure 1. Distribution functions are discretized according to these sets, giving the discrete partial distribution functions $f_i$. The method is implemented on a Cartesian lattice with temporal and spatial mesh width $\Delta t$ and $\Delta x$ of size 1.

<figure><img src="image_4.png" /><figcaption>Figure 1. Velocity discretization at a lattice node. The scheme is a 2-dimensions, 9 discrete velocities set (D2Q9).</figcaption></figure>

The discretized Boltzmann equation is solved in two steps called collision and streaming. During the first one, post-collision distribution functions $f'_i$ are computed using the collision operator $\Omega$. The second step is the streaming of post-collision distribution func-
tions onto the neighbouring nodes pointed by the directions of propagation. Those two steps are illustrated in Figure 2 for a D2Q9 velocity scheme.

<figure><img src="image_5.png" /><figcaption>Figure 2. LBM iteration progress. (a) Collide step. (b) Streaming step.</figcaption></figure>

Macroscopic quantities of the fluid can be recovered using distribution functions. Density at a node is calculated with Equation (2) and velocity with Equation (3):

$$ \rho = \sum_{i=0}^{n} f_i, \qquad (2) $$

$$ \rho u = \sum_{i=0}^{n} c_i f_i, \qquad (3) $$

where n is the number of directions of propagation of the velocity scheme. Simulations presented in this study use a 3-dimensional set with 27 directions of propagation (D3Q27).

There are many collision operators. Simulations presented in this article are performed with a recursive and regularized collision operator (RR). The RR procedure was introduced by [23] to improve the stability and accuracy of the LBM. The collision operator is used with two-relaxation times (RR-TRT) and is available in Palabos under the name of *Consistent Smagorinsky Complete Regularized TRT dynamics*. It is used at order 6. Boundary conditions use second-order accurate finite-difference approximations for velocity gradients scheme and can be used to apply velocity/pressure Dirichlet/Neumann conditions (Skordos [24]). This boundary condition is the most stable available in Palabos and is well suited for high Reynolds numbers simulations [25].

Uniform grids are not adapted to model complex geometries at high Reynolds numbers. To overcome this, static mesh refinement is used. The mesh refinement procedure is described in Lagrava et al. [21] and uses a multi-block approach. A convective scaling is set between various levels: $\Delta x_{coarse} = 2\Delta x_{fine}$, $\Delta t_{coarse} = 2\Delta t_{fine}$.

An LES approach is used and the subgrid model is a static Smagorinsky model [26]. It adds a turbulent viscosity calculated with Equation (4):

$$ v_{sgs}(\mathbf{x}, t) = (C_s' \Delta x)^2 (2|\bar{S}(\mathbf{x}, t)|^2)^{1/2}, \qquad (4) $$
where $\bar{S}$ is the filtered strain rate tensor. The implementation of the static Smagorinsky model in the LBM is described in Malaspinas and Sagaut [27]. The coefficient $C'_s$ is linked to the Smagorinsky constant $C_s$ by using the Van Driest damping function, described in Equation (5), in order to cancel the turbulent viscosity close to the turbine's surface. The distance normal to the surface is $y$, $v$ is the molecular viscosity of the fluid and $u_τ$ is the friction velocity. The value $C_s = 0.18$ has been chosen for this study.

$$C'_s = C_s(1 - \exp(-\frac{y u_τ}{26 \nu})). \qquad (5)$$

## 2.2. Reference Experiment

The experiment as reference is the one of Mycek et al. [6]. He studied the wake and performances of a three-bladed tidal turbine model with two turbulent inflow conditions. The experiment has been carried out at Ifremer's flume tank. The tidal turbine was developed by Ifremer and is illustrated in Figure 3. It is placed in the middle of the tank and the blockage ratio is approximately 4.8%.

<figure><img src="image_6.png" /><figcaption>Figure 3. Tidal turbine model, Mycek et al. [6].</figcaption></figure>

Two turbulent inflow conditions were tested. The first one with $I = 3\%$, where flow straighteners were used. The second one with $I = 15\%$, where no flow straighteners were used. Several tip speed ratios (TSR) were also tested. The TSR of the numerical simulations is set close to the operating point at TSR = 3.67. Characteristics chosen for the LBM-LES study are listed in Table 1.

<table><caption>Table 1. Experimental configurations from Mycek et al. [6] studied with the LBM-LES.</caption><thead><tr><th>Quantity</th><th>Value</th><th>Unit</th></tr></thead><tbody><tr><td>Blade radius R = D/2</td><td>0.35</td><td>m</td></tr><tr><td>Flow velocity &lt; U<sub>x∞</sub> &gt;</td><td>0.8</td><td>m.s<sup>-1</sup></td></tr><tr><td>Tip speed ratio TSR = ωR/ &lt; U<sub>x∞</sub> &gt;</td><td>3.67</td><td>~</td></tr><tr><td>Reynolds number Re = &lt; U<sub>x∞</sub> &gt; R/v</td><td>280 × 10<sup>3</sup></td><td>~</td></tr><tr><td>Turbulence intensity I</td><td>[3, 15]</td><td>%</td></tr></tbody></table>

Performances of the turbine are evaluated with the thrust coefficient $C_T$ and power coefficient $C_P$, defined in Equations (6) and (7). $F_x$ is the axial force, $M_x$ the axial torque and $\omega$ the rotational speed of the rotor.

$$C_T = \frac{F_x}{0.5 \rho \Pi R^2 < U_{x\infty} >^2}. \qquad (6)$$

$$C_P = \frac{M_x \omega}{0.5 \rho \Pi R^2 < U_{x\infty} >^3}. \qquad (7)$$

The axial force is measured with a load cell placed at the fixing point of the turbine mast, so $F_x$ is the force applied to both stator and rotor. The torque is measured with a
torque sensor located on the rotor axis. Wake measurements are made with laser Doppler velocimetry (LDV) and the frame origin is the rotor centre. The flow velocity $\mathbf{U}$ and the location of the measurements $x$ are presented in dimensionless units with $\mathbf{U}^* = \mathbf{U} / \langle U_{x\infty} \rangle$ and $x^* = x / R$. The experimental set-up is presented in Figure 4.

<figure><img src="image_7.png" /><figcaption>Figure 4. Experimental set-up of Mycek et al. [6]. The of the coordinate system’s origin is the rotor centre.</figcaption></figure>

### 2.3. Modelling the Turbine

Boundary conditions for straight walls introduced in Section 2.1 are not adequate for modelling complex geometries like tidal turbines. Using such conditions would end with a staircase approximation of the surface and an inaccurate wake. The technique applied here is the immersed boundary method (IBM), which is well suited for modelling complex geometries [11]. The IBM was first introduced by Peskin [28] for simulating blood flow in the heart. The first use of IBM in an LBM simulation was done by Feng and Michaelides [29]. In Palabos IBM-LBM, the solid surface is discretized into a set of Lagrangian points. There should be enough points to correctly describe the solid and the spacing between points should not be superior to $\Delta x$. The IBM uses interpolations to compute the fluid velocity and density at each Lagrangian points $P$ (Figure 5). Quantities are interpolated from lattice macroscopic quantities computed from pre-collision distribution functions. The solid boundary is modelled through a force acting on the fluid at the Lagrangian point’s location. The forcing method used in Palabos is based on [30]. The momentum at a Lagrangian point is written:

$$ \rho_P \mathbf{U}_{\text{boundary}} = \rho_P \mathbf{U}_P + F_P, \qquad (8) $$

where $\mathbf{U}_{\text{boundary}}$ is the velocity of the boundary and $F_P$ is the forcing term. The force $F_P$ can be written: $F_P = \rho_P(\mathbf{U}_{\text{boundary}} - \mathbf{U}_P)$. The no-slip condition is then enforced by redistributing the force $F_P$ to the fluid momentum (Equation (9)):

$$ \rho(x) \mathbf{U}(x) = \rho(x) \mathbf{U}(x) + \sum_P \tau S_P F_P D(x - x_P), \qquad (9) $$

where $S$ is the surface portion of the boundary represented by $P$ and $D(x)$ is the Dirac distribution function. The relaxation frequency $\tau$ comes from the collision operator. This procedure is applied to every Lagrangian points of the boundary. The collision step is then executed with a modified momentum for the lattice nodes surrounding the Lagrangian points. More advanced IBM methods such as [31] could be explored in the future.
<figure><img src="image_8.png" /><figcaption>Figure 5. Solid surface represented with Lagrangian points on a lattice.</figcaption></figure>

Simulating high Reynolds numbers flow with LBM-LES can be very resource consum-
ing. Indeed, LES mesh size constraints at the wall are almost the same as direct numerical
simulation (DNS) ones [32]. Using grid refinement helps but Cartesian meshes have many
nodes compared to structured or unstructured meshes. To reduce the computational cost
of simulations, a wall model is used. It is implemented using a pre-existing wall model of a
beta version of Palabos. It is based on the model described in Ruffin et al. [33] and uses
Spalding’s formulation of the boundary layer [34]:

$$
y^+ = u^+ + e^{-\kappa B} \left[ e^{\kappa u^+} - 1 - \kappa u^+ - \frac{(\kappa u^+)^2}{2} - \frac{(\kappa u^+)^3}{6} \right], \quad (10)
$$

where $y^+$ is the normalized distance to the solid surface: $y^+ = u_τy/v$. The velocity $u^+$ is calculated with: $u^+ = U/u_τ$. The variable $κ$ is the von Kármán constant and $B$ is a constant. The Newton–Raphson method is used for computing $u_τ$. The turbulent viscosity is computed like in Ruffin et al. [33] and applied to the lattice by modifying the relaxation time like in Malaspinas and Sagaut [35]. Spalding’s law of the wall needs two parameters: $κ$ and $B$. For the simulations presented thereafter, $κ = 0.42$. The parameter $B$ can slightly vary depending on the mesh size or the upstream turbulence, it is typically equal to $B = 2.6 \pm 0.2$.

2.4. *Generating Realistic Inflow Conditions*

There are several ways of implementing a realistic turbulent inflow condition. In Gron-
deau et al. [36], the synthetic eddy method (SEM) developed by Jarin et al. [37] has been
implemented in Palabos. The SEM generates a random turbulent inflow condition that
matches a specified average Reynolds tensor and average axial velocity. It gave good results
at low Reynolds numbers. As the Reynolds number increases, simulations can become un-
stable due to density fluctuations generated by the SEM boundary condition. To overcome
this issue, the divergence-free SEM (DFSEM) of Poletto et al. [22] has been implemented in
Palabos and is used here. The DFSEM is closed to the SEM except for the calculation of
velocity fluctuations. The algorithm can be summarized by the following steps:

1. Definition of the DFSEM parameters: average stream-wise velocity < U >, eddy length L, Reynolds tensor, number of eddies N;

2. Creation of the eddy box around the inlet of the LBM domain, see Figure 6;

3. Creation of N eddies at a random location and with random intensities;

4. Displacement of eddies in the stream-wise direction;

5. Computation of eddies contributions to the fluctuations of the average stream-wise velocity;

6. Repeat 4 and 5.

When an eddy is leaving the DFSEM domain it is deleted. A new eddy is created at
$x = 0$ in the DFSEM reference frame with random $y, z$ positions and intensities. Each eddy
adds velocity fluctuations to the surrounding inlet nodes. The radius of an eddy influence
zone is defined by the eddy length $L$. Fluctuations are computed using Equation (11).

$$
u'_{\text{node}} = \sqrt{\frac{1}{N}} \sum_{k=1}^{N} \frac{q_L(|\mathbf{r}_k|)}{|\mathbf{r}_k|^3} \mathbf{r}_k \times \alpha_k, \quad (11)
$$
where $q_L$ is the shape function and $r_k = \frac{x-x_k}{L}$, with $x$ the position of the node and $x_k$ the position of an eddy. The shape function used is defined in Equation (12).

$$q_L(a) = \sqrt{\frac{16V_B}{15\pi L^3}} \sin(\pi a)^2 a, \qquad (12)$$

where $V_B$ is the volume of the DFSEM domain. Polleto et al. (2013) suggested to use an eddy density $d$ around 1 in order to have a good compromise between computational cost and accuracy. The eddy density is defined in Equation (13).

$$d = \frac{L^3 N}{V_B}. \qquad (13)$$

The DFSEM boundary condition is applied through the Dirichlet boundary condition of Palabos. The objective is to generate synthetic turbulence as close as possible to the turbulence of Mycek et al.'s [6] experiment. In addition to the data provided by Mycek et al. [6], an analysis of the upstream and downstream turbulence was made by Medina et al. [38]. The general definition of the turbulence intensity, noted $I$ or $I3D(\%)$, is based on the three components of the velocity fluctuations, and is given by Equation (14), where $<u_i^2>$ are the time-averaged velocity fluctuations and $<U_i>$ are the time-averaged velocity components. In addition, the stream-wise turbulence intensity $I1D(\%)$ is based on the velocity component in the current direction. Similarly, the two-dimensional turbulence intensity $I2D(\%)$ is based on the velocity fluctuation in the flow direction and in the transverse direction. The turbulence rates $I_{1D}$ and $I_{2D}$ are defined in equation Equations (15) and (16). Measured turbulence rate $I$, $I_{1D}$ and $I_{2D}$ were also provided by Ifremer. Depending on the source, the turbulence rate $I$ is found between 12.5% and 15% for the case where no flow straighteners are used. According to those sources, the turbulence rate is constant in the stream-wise direction. When using a synthetic turbulent inflow condition, one should verify that the turbulence intensity is conserved over a long enough distance. Indeed, for a high turbulence rate, the artificial turbulent structures rapidly lose their intensity and the turbulence rate decreases according to the stream-wise direction. Accordingly, we choose to generate a flow with a turbulence rate of $I = 12.5\%$. This case is denominated "$I = 12.5\%$" from now on. When flow straighteners are used, the turbulence rate is close to $I = 3\%$. This value is kept for this study and this case is denominated "$I = 3\%$".

$$I = I_{3D} = 100 \sqrt{\frac{1/3(<u_x^2> + <u_y^2> + <u_z^2>)}{<U_x>^2 + <U_y>^2 + <U_z>^2}} \qquad (14)$$

$$I_{1D} = 100 \sqrt{\frac{<u_x^2>}{<U_x>^2}} \qquad (15)$$

$$I_{2D} = 100 \sqrt{\frac{0.5(<u_x^2> + <u_y^2>)}{(<U_x>^2 + <U_y>^2)}} \qquad (16)$$

The study made by Medina et al. [39] has brought to light that the turbulence of the flume tank is fully developed. The length scale of turbulent structures was estimated to be between 0.25 m and 0.6 m for $I = 12.5\%$ and between 0.65 m and 0.9 m for $I = 3\%$. The turbulence is isotropic for $I = 12.5\%$, whereas for $I = 3\%$, the stream-wise velocity fluctuation is stronger than the other two components. The implemented DFSEM does not allow strong anisotropy ratios and the velocity fluctuations are smoothed over the three directions.

The DFSEM is used to generate a realistic inflow boundary condition. The implemented DFSEM needs four parameters to work: the convection velocity of the DFSEM eddies, the number of DFSEM eddies, the size of DFSEM eddies and the imposed Reynolds tensor. The velocity of the DFSEM eddies is equal to the average upstream stream-wise velocity
$<U_{x_\infty}>$. To define the number of eddies, Polleto et al. (2013) suggested using an eddy density $d$ around 1 in order to have a good compromise between computational cost and accuracy. The eddy density is defined in Equation (13). It is chosen to $d = 1$ for the case with $I = 3\%$. For the case with $I = 12.5\%$, the eddy density is chosen to be $d = 2$ to avoid spatial disparity. The size of DFSEM eddies and the imposed Reynolds tensor are chosen carefully after dedicated studies, presented in Section 3.1. Those studies are carried out without the turbine in order to save computational resources. It allows for longer simulated durations, enabling advanced statistics like spatial autocorrelation. For those studies, the numerical domain is 4 m high, 4 m wide and 12 m long. The mesh is uniform and, unless stated otherwise, $\Delta x = 0.035$ m.

<figure><img src="image_9.png" /><figcaption>**Figure 6.** Working principle of the DFSEM in an LBM simulation.</figcaption></figure>

## 2.5. Setup of the Turbine Modelling

This subsection introduces the numerical domains used for modelling the turbine presented in Section 2.2. In order to have elements regarding the mesh sensitivity of the approach, three mesh configurations are studied. The static mesh refinement of Palabos is used here. The three meshes are summarized in Table 2. The boundary layer model used is only valid in the logarithmic layer, viscous sub-layer and buffer layer [33]. Therefore, the coarsest mesh is chosen in order to have a $\Delta x_{max}^+ \approx 300$ for TSR = 3.67. The finest level of the mesh has the highest number and the coarsest level has level number 0.

<table><caption>**Table 2.** Meshes selected for the LBM mesh sensitivity study. The finest level is the numerical domain surrounding the rotor of the turbine. h.CPU is the wall-clock time multiplied by the number of CPU used.</caption><thead><tr><th>Mesh</th><th>Finest Level</th><th>Δx at Finest Level</th><th>Number of Nodes<br>Revolution in h.CPU</th><th>Computational Time Per</th></tr></thead><tbody><tr><td>No. 1</td><td>6</td><td>1.39 × 10<sup>-3</sup> m</td><td>83 × 10<sup>6</sup></td><td>450 h.CPU</td></tr><tr><td>No. 2</td><td>6</td><td>1.04 × 10<sup>-3</sup> m</td><td>60 × 10<sup>6</sup></td><td>1300 h.CPU</td></tr><tr><td>No. 3</td><td>5</td><td>7.32 × 10<sup>-4</sup> m</td><td>124 × 10<sup>6</sup></td><td>1910 h.CPU</td></tr></tbody></table>

Meshes around the turbine are shown in Figure 7. Only the area swept by the rotor is refined at the finest level. Indeed, the flow around the blades is the area where the Reynolds number is the highest and where the boundary layer is the thinnest. The finest mesh volume for mesh No. 2 and 3 is much smaller than for mesh No. 1. After running simulations with mesh No. 1, it seemed impractical to use the same finest mesh volume for mesh No. 2 and 3. In order to lower the computational cost, considerable work has been done to lower as much as possible the number of nodes in the finest domain. Since convective scaling is used, all three cases are run with a constant Mach number of $Ma = 0.38$.
<figure><img src="image_10.png" /><figcaption>**Figure 7.** Mesh around the tidal turbine [6] used for LBM-LES simulations. The domain is centred at the rotor centre. (a): mesh No. 1 is represented in cubes of 14 × 14 × 14 nodes. (b,c): mesh No. 2 and 3 are represented in cubes of 16 × 16 × 16 nodes.</figcaption></figure>

Full meshes are visible in Figure 8. The turbine is located more than 5D downstream of the inlet. Considering the small blockage ratio, the flume tank is not modelled and the turbine is 6D away from lateral boundaries. The wake is divided into two areas, the close wake up to 4D and the far wake from 4D and up to 11.5D. Because of how the mesh refinement works in Palabos, the mesh size in the wake of the three cases is slightly different. Table 3 summarizes the mesh size in the wake of the turbine. Liu and Hu [40] studied the turbine of Mycek et al. [6] with an ALM unsteady-RANS model. A converged solution have been obtained with a mesh size in the wake of 8.75 × 10⁻³ m. This is slightly larger than the wake mesh size of Mesh No 2. A three-dimensional view of the domain is shown in Figure 9.

The Dirichlet boundary condition modified with DFSEM is used on the inlet and with $U = (< U_{x∞} >, 0, 0)$ on sides boundaries. The outlet is a pressure-imposed boundary condition. On every boundary, except for the inlet, a pressure-wave absorbing function is applied to distribution functions. This function is described in Xu and Sagaut [41]. Simulations are carried out at the *Centre Régional Informatique et d'Application Numérique de Normandie* (CRIANN, France) on Broadwell CPU cores (2.4 GHz).
<figure><img src="image_11.png" /><figcaption>(a) Mesh No. 1—xy view</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>(b) Mesh No. 2—xy view</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>(c) Mesh No. 3—xy view</figcaption></figure>

**Figure 8.** Mesh used in LBM-LES simulations of the tidal turbine [6]. Domain is centred at the rotor centre. (a): mesh No. 1 is represented in cubes of 14 × 14 × 14 nodes. (b,c): mesh No. 2 and 3 are represented in cubes of 16 × 16 × 16 nodes.

<table><caption>**Table 3.** Mesh sizes in the wake of the turbine.</caption><thead><tr><td>Mesh</td><td>Close Wake Level</td><td>Δx Close Wake</td><td>Far Wake Level</td><td>Δx Far Wake</td></tr></thead><tbody><tr><td>No. 1</td><td>3</td><td>5.58 × 10<sup>-3</sup> m</td><td>2</td><td>11.16 × 10<sup>-3</sup> m</td></tr><tr><td>No. 2</td><td>3</td><td>8.30 × 10<sup>-3</sup> m</td><td>2</td><td>16.60 × 10<sup>-3</sup> m</td></tr><tr><td>No. 3</td><td>3</td><td>5.58 × 10<sup>-3</sup> m</td><td>2</td><td>11.16 × 10<sup>-3</sup> m</td></tr></tbody></table>
<figure><img src="image_14.png" /><figcaption>Figure 9. Numerical domain of LBM-LES simulations of the tidal turbine [6]. Mesh No. 1 is displayed on the boundaries of the domain, it is represented in cubes of 14 × 14 × 14 nodes. The domain is centred at the rotor centre and dimensions are scaled by the turbine diameter D.</figcaption></figure>

### 3. Inlet Calibration and Validation

This section first presents the results of the study dedicated to the choice of DFSEM parameters. The chosen sets of DFSEM parameters are then used for the LBM-LES modeling of the turbine. Results of those simulations are presented afterwards.

#### 3.1. Upstream Turbulence for the Tidal Turbine

The choice of DFSEM parameters is first investigated for the case with $I = 3\%$ and then for the case with $I = 12.5\%$. Stream-wise evolutions of quantities such as turbulence rate, power spectral density or integral length scale are investigated.

##### 3.1.1. DFSEM, Low Upstream Turbulence Rate Case ($I = 3\%$)

For this case, the choice of parameters is straightforward. The size of eddies is $L = 1.0$ m, which is close to the maximum of the measured integral length scale $L_{max} = 0.9$ m. $L$ is equal to the diameter of eddies since they are sphere shaped. The Reynolds tensor is chosen as a diagonal matrix and all diagonal terms are equal: $u_x'^2 = u_y'^2 = u_z'^2 = 8.928 \times 10^{-4} \text{ m}^2 \cdot \text{s}^{-2}$.

Figure 10 shows the evolution of the turbulence rate $I_{1D}$ according to the stream-wise direction. The turbulence rate is almost conserved along the flume tank. An autocorrelation function is used to compute the integral length scale $L$ of turbulent structures. This function is applied in the plane of normal $z$ at $z = 2$ m and according to the direction $y$. It is defined in Equation (17):

$$AC(x) = \frac{3}{N_{AC}(\langle u_x' \rangle + \langle u_y' \rangle + \langle u_z' \rangle)} \sum_{i=0}^{N_{AC}} u(x_i) \cdot U(x_0), \quad (17)$$

where $N_{AC}$ is the number of iterations during which the autocorrelation is applied. The point $x_0$ is located at the centre of the domain. Figure 11 shows the evolution of the estimated integral length scale according to the stream-wise direction.
<figure><img src="image_15.png" /><figcaption>**Figure 10.** Turbulence rate $I_{1D}$ at the centre of the LBM domain using DFSEM with $L = 1.0$ m and $u_i'^2 = 8.928 \times 10^{-4} \text{ m}^2\cdot\text{s}^{-2}$. $x$ is the distance from the inlet of the LBM domain.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Figure 11.** Results of the autocorrelation function applied to the velocity field computed from a LBM-LES simulation using DFSEM with $L = 1.0$ m and $u_i'^2 = 8.928 \times 10^{-4} \text{ m}^2\cdot\text{s}^{-2}$. $x$ is the distance from the inlet of the LBM domain. $y$ is the distance from the centre of the domain in the span-wise direction.</figcaption></figure>

The integral length scale of turbulent structures generated by the DFSEM boundary condition is conserved along the stream-wise direction and is close to the prescribed eddy length $L$.

Turbulent flows are complex phenomena where coherent structures of several scales evolve. Large structures are the most energetic and energy goes from big structures to small structures. This energy cascade is observable with a power spectral density (PSD) calculated from a sampling of the velocity. The PSD from a developed turbulent flow has a slope of -5/3, characteristic of the energy cascade. The stream-wise velocity component is recorded at the centre of the flume at $f = 112$ Hz in $x = 2.0, 4.0$ and $9.0$ m. The PSD calculated from this recording is plotted in Figure 12.

<figure><img src="image_17.png" /><figcaption>**Figure 12.** Power spectral density calculated from the stream-wise velocity component computed from a LBM-LES simulation using DFSEM with $L = 1.0$ m and $u_i'^2 = 8.928 \times 10^{-4} \text{ m}^2\cdot\text{s}^{-2}$. $x$ is the distance from the inlet of the LBM domain.</figcaption></figure>
A drop in the PSD for frequencies higher than 1 Hz is observed. This drop has two origins: only big structures are generated by the DFSEM boundary condition and the relatively coarse mesh does not allow for small structures to be conserved. The PSD would certainly be closer to the -5/3 slope with a finer mesh but the volume upstream of the turbine is too big for a finer mesh size. For the frequencies inferior to 1 Hz, the PSD is close to the -5/3 slope.

To conclude, the turbulent flow generated with the DFSEM is close to the experimental flow upstream of the turbine for I = 3%. The chosen set of parameters is: $u_i'^2 = 8.928 \times 10^{-4} \text{ m}^2 \cdot \text{s}^{-2}$ and L = 1.0 m. Those parameters are used to model the tidal turbine with the LBM-LES with I = 3% in Section 3.2.1.

### 3.1.2. DFSEM, High Upstream Turbulence Rate Case (I = 12.5%)

The turbulent flow generated by the DFSEM loses its intensity as it travels in the flume. This effect increases with the turbulence rate. Consequently, the choice of parameters is less straightforward for high inflow turbulence rates and several sets of parameters must be tested. Three sizes of eddies are tested: L = 0.25, 0.66 and 1.2 m. Two Reynolds tensors are tested. They are both diagonal and diagonal terms are equal to: $u_x'^2 = u_y'^2 = u_z'^2 = 0.01638 \text{ m}^2 \cdot \text{s}^{-2}$. Two mesh sizes are tested: Δx = 0.035 and 0.025 m. The evolution of the turbulence rate $I_{1D}$ according to the stream-wise direction, for several sets of parameters, is plotted in Figure 13.

<figure><img src="image_18.png" /><figcaption>**Figure 13.** Turbulence rate $I_{1D}$ at the centre of the LBM domain for different parameters of the DFSEM. The targeted turbulence rate is 12.5%. x is the distance from the inlet of the LBM domain.</figcaption></figure>

In cases where L is equal to 0.25 m and 0.6 m both, have a strong decay of the turbulence rate. The eddy size seems to have a great influence on the decay of the turbulence rate. Cases with $u_i'^2 = 0.01638 \text{ m}^2 \cdot \text{s}^{-2}$ and $u_j'^2 = 0.04000 \text{ m}^2 \cdot \text{s}^{-2}$ converge toward the same turbulence rate. Increasing the velocity fluctuations at the inlet is not a solution. A finer mesh size does not have a significant influence on the decay of the turbulence rate. It appears that the best solution to conserve a high enough turbulence rate is to increase the size of eddies.

The autocorrelation function used previously is applied to the case with L = 1.2 m. Results are plotted in Figure 14. The estimated integral length scale is conserved and close to the prescribed size of eddies. Following Blackmore et al. [42] and Ebdon et al. [12] the turbulence integral length scale has a great influence on the behaviour of the turbine. Because the prescribed size L is conserved in the canal, a compromise must be chosen. The set L = 1.2 m and $u_i'^2 = 0.01638 \text{ m}^2 \cdot \text{s}^{-2}$ allows conserving the turbulence rate above 10%. This set is chosen for the case I = 12.5%. Since L is twice as large as L, this set of parameters may induce some discrepancies. Table 4 presents the evolution of the turbulence $I_{1D}$, $I_{2D}$, and $I_{3D}$ defined in Equations (14)–(16). The slope linear regression through the data points employed in [43,44] was used to estimate the rate of change in the distance from the inlet of the LBM domain with turbulence rates. It was observed that as the distance from the inlet
of the LBM domain increases, I1D (%) decreases at the rate of −0.547, I2D (%) decreases at the rate of −0.533, while I3D (%) decreases at the rate of −0.543. Rates at which the 1D, 2D and 3D turbulence intensities decreases are extremely close, which tends to indicate that the turbulence generated by the DFSEM is isotropic.

<figure><img src="image_19.png" /><figcaption>**Figure 14.** Results of the autocorrelation function applied to the velocity field computed from a LBM-LES simulation using DFSEM with $L = 1.2$ m, $u_i^{\prime 2} = 0.01638 \text{ m}^2 \cdot \text{s}^{-2}$ and $\Delta x = 0.035$ m. $x$ is the distance from the inlet of the LBM domain. $y$ is the distance from the centre of the domain in the span-wise direction.</figcaption></figure>

<table><caption>**Table 4.** Turbulence rates from a LBM-LES simulation using DFSEM with $L = 1.2$ m, $u_i^{\prime 2} = 0.01638 \text{ m}^2 \cdot \text{s}^{-2}$ and $\Delta x = 0.025$ m. $x$ is the distance from the inlet of the LBM domain. Turbulence rates are calculated at the centre of the LBM domain.</caption><thead><tr><th>x (m)</th><th>I<sub>1D</sub> (%)</th><th>I<sub>2D</sub> (%)</th><th>I<sub>3D</sub> (%)</th></tr></thead><tbody><tr><td>0.1</td><td>13.00</td><td>12.76</td><td>12.88</td></tr><tr><td>2.0</td><td>12.65</td><td>12.28</td><td>12.25</td></tr><tr><td>6.0</td><td>10.34</td><td>10.06</td><td>9.98</td></tr><tr><td>9.0</td><td>8.25</td><td>8.11</td><td>8.13</td></tr><tr><td>Slope</td><td>-0.547</td><td>-0.533</td><td>-0.543</td></tr></tbody></table>

The PSD of the stream-wise velocity is plotted in Figure 15. Close to the inlet, a drop is observed at the frequency of DFSEM eddies. This drop quickly disappears and the slope of the PSD converges toward the −5/3 slope around $x = 4.0$ m.

<figure><img src="image_20.png" /><figcaption>**Figure 15.** Power spectral density calculated from the stream-wise velocity component computed from a LBM-LES simulation using DFSEM with $L = 1.2$ m, $u_i^{\prime 2} = 0.01638 \text{ m}^2 \cdot \text{s}^{-2}$ and $\Delta x = 0.035$ m. $x$ is the distance from the inlet of the LBM domain.</figcaption></figure>
To conclude, a compromise has been made in order to conserve the turbulence rate over a long enough distance. The chosen set of parameters is: $u_i^2 = 0.01638 \text{ m}^2.\text{s}^{-2}$ and $L = 1.2 \text{ m}$. A distance of convergence of approximately 4 m has been observed. The turbine location, presented in Section 2.5, respects this conclusion.

## 3.2. Tidal Turbine Modelling

Results from the LBM-LES simulations of the turbine are presented here, at first with a low upstream turbulence rate, then with a high upstream turbulence rate. Simulations results are compared with experimental results, investigated wake quantities are the average axial velocity and the average turbulence rate. Force coefficients from simulations are also compared to experimental ones.

### 3.2.1. Turbine, Low Upstream Turbulence Rate Case ($I = 3\%$)

The simulation of the tidal turbine, described in Section 2.2, with an upstream turbulence rate of $I = 3\%$ is performed with an LBM-LES code. Quantities in the wake and force coefficients are averaged over 12 s at a sampling rate of 398 Hz. A convergence duration of 16 s is respected before the calculation of average quantities.

Table 5 presents force coefficients from both the experiment and simulations. It appears that mesh No. 2 and 3 are converged in thrust coefficient $C_T$. Mesh No. 1 gives a relatively good $C_T$ prediction with an error of less than 10%. The difference in power coefficient $C_p$ keeps decreasing with the mesh size. Nevertheless, the error in $C_p$ is always lower than 20%, which is relatively good considering the fact that the boundary layer is not fully resolved. Mesh No. 2 is a compromise between accuracy and computational cost. For the rest of the study, simulations are carried out with mesh Nos. 1 and 2. Mesh No. 1 is kept because it cost three times less than mesh No. 2 and may give good wake prediction.

<table><caption>**Table 5.** Turbine force coefficients from [6] and from LBM-LES simulations at TSR = 3.67 and $I = 3\%$. Differences are calculated between the experiment and simulations.</caption><thead><tr><th>Mesh</th><th>C<sub>T</sub></th><th>Difference</th><th>C<sub>P</sub></th><th>Difference</th></tr></thead><tbody><tr><td>Experiment</td><td>0.805</td><td>\</td><td>0.420</td><td>\</td></tr><tr><td>No. 1</td><td>0.74</td><td>-7.8%</td><td>0.50</td><td>19.7%</td></tr><tr><td>No. 2</td><td>0.79</td><td>-1.2%</td><td>0.48</td><td>13.5%</td></tr><tr><td>No. 3</td><td>0.82</td><td>1.4%</td><td>0.46</td><td>9.1%</td></tr></tbody></table>

Average axial velocity profiles and turbulence intensity profiles are plotted in Figures 16 and 17. At 1.2D downstream of the turbine, in the close wake, the difference between LBM-LES simulations with mesh No. 1 and mesh No. 2 is non-negligible. The average axial velocity is better calculated with mesh No. 2. The interaction between the mast and the wake creates an asymmetry of the wake and is visible on the experimental turbulence profile at $x = 1.2D$. This asymmetry is not on the simulations profiles. This interaction is probably not properly computed by the boundary layer model. The turbulence peak around $y = \pm 0.5D$ is correctly predicted by both LBM-LES simulations, it is most probably coming from tip-vortices. Tip-vortices can be observed in Figure 18. After 2D, both the axial velocity and the turbulence intensity are correctly calculated by LBM-LES simulations with mesh Nos. 1 and 2. On the turbulence profile obtained with mesh No. 1, a small irregularity is observed 5D downstream of the turbine. Considering the shape of this irregularity, it may come from the tip-vortices. This could indicate that either the tip-vortices have slightly different behaviour in this simulation, or the quantities are not averaged over a duration long enough.
<figure><img src="image_21.png" /><figcaption>**Figure 16.** Average axial velocity profiles in the wake of the tidal turbine at TSR = 3.67 and I = 3%. Results are from experiment [6] and from LBM-LES simulations with mesh Nos. 1 and 2. The complete set of profiles is in Appendix A.</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>**Figure 17.** Average turbulence intensity profiles I<sub>2D</sub> in the wake of the tidal turbine at TSR = 3.67 and I = 3%. Results are from experiment [6] and from LBM-LES simulations with mesh Nos. 1 and 2. The complete set of profiles is in Appendix A.</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>**Figure 18.** Three-dimensional view of the turbine with an upstream turbulence rate of I = 3%, instantaneous U<sub>z</sub> velocity. Tip-vortices are observable with iso-contour of the lambda2 criterion. LBM-LES simulation with mesh No. 1.</figcaption></figure>
### 3.2.2. Turbine, High Upstream Turbulence Rate Case (I = 12.5%)

The simulation of the tidal turbine described in Section 2.2 with an upstream turbulence rate of *I* = 12.5% is performed with an LBM-LES code. As stated in Section 3.2.1, LBM-LES simulations are only carried out with mesh No. 1 and 2. Quantities in the wake and force coefficients are averaged over 16 s at a sampling rate of 398 Hz. It is slightly more than for *I* = 3% because a higher turbulence rate requires a longer sampling duration. Convergence lasts 16 s.

Table 6 shows force coefficients from the LBM-LES simulations and experiment. Mesh Nos. 1 and 2 give good thrust predictions. However, they both overestimate the power coefficient and gives relatively close power predictions. It has been pointed out in Section 3.1 that the chosen size of eddies could induce discrepancies. Indeed, Blackmore et al. [42] have observed that increasing the integral length scale could increase force coefficients. This could explain the observed differences in $C_p$.

<table><caption>**Table 6.** Turbine force coefficients from [6] at TSR = 3.67 and *I* = 12.5% and from LBM-LES simulations at TSR = 3.67 with *I* = 12.5%.</caption><thead><tr><th>Mesh</th><th>C<sub>T</sub></th><th>Difference</th><th>C<sub>P</sub></th><th>Difference</th></tr></thead><tbody><tr><td>Experiment</td><td>0.73</td><td>\</td><td>0.37</td><td>\</td></tr><tr><td>No. 1</td><td>0.70</td><td>-4.4%</td><td>0.462</td><td>23.0%</td></tr><tr><td>No. 2</td><td>0.75</td><td>3.4%</td><td>0.48</td><td>28.6%</td></tr></tbody></table>

Average axial velocity profiles and turbulence intensity $I_{2D}$ profiles are plotted in Figures 19 and 20. Differences are small between the two meshes, except at 1.2D where the average axial velocity is a bit more accurate with mesh No. 2. At 1.2D downstream, a slight shift of the wake in the *y* direction is observed on the experimental profile. Data provided by Ifremer show that the upstream flow is not perfectly uni-axial when no flow straighteners are used. A small *y* component exists and may affect the direction of propagation of the wake. It could explain the difference. Except for this difference, the axial velocity profiles from the simulations are close to the experimental ones.

<figure><img src="image_24.png" /><figcaption>**Figure 19.** Average axial velocity profiles in the wake of the tidal turbine at TSR = 3.67 and *I* = 12.5%. Results are from experiment [6] and from LBM-LES simulations with mesh Nos. 1 and 2. The complete set of profiles is in Appendix B.</figcaption></figure>
<figure><img src="image_25.png" /><figcaption>**Figure 20.** Average turbulence intensity profiles *I*<sub>2D</sub> in the wake of the tidal turbine at *TSR* = 3.67 and *I* = 12.5%. Results are from experiment [6] and from LBM-LES simulations with mesh Nos. 1 and 2. The complete set of profiles is in Appendix B.</figcaption></figure>

The maximum turbulence intensity predicted by the LBM-LES simulation is identical to the experimental one. There is a lack of turbulence elsewhere in the close wake. For instance, in the $-1 < y^* < 1$ interval at $x = 1.2D$, the average experimental $I_{2D}$ intensity is 22%, whereas it is 16% in the simulation with mesh No. 1. There is no obvious explanation for this difference but several sources may be suggested. The boundary layer model may react badly to the high rate of upstream turbulence or the integral length scale may influence the production of turbulence. This difference decreases in the far wake. Despite not being revealed through a located peak on averaged profiles, tip-vortices are observed in Figure 21. The ambient turbulence's influence is seen to be notable as tip-vortices rapidly lose consistency.

<figure><img src="image_26.png" /><figcaption>**Figure 21.** Three-dimensional view of the turbine with an upstream turbulence rate of *I* = 12.5%, instantaneous *U*<sub>z</sub> velocity. Tip-vortices are observable with iso-contour of the *lambda2* criterion. LBM-LES simulation with mesh No. 1.</figcaption></figure>

The LBM-LES simulations have proven to be reliable to predict the wake of tidal turbines with a high upstream turbulence intensity. A more detailed study of the wake is proposed in the next section.

**4. Wake Analysis**

Velocity deficit and increase of turbulence have been observed in turbines' wake. Thus, previous results have shown that the upstream turbulence has a great influence over the wake of a tidal turbine. In this section, a cross-comparison between simulations at *I* = 3, 8
and 12.5% is carried out. Then a spectral analysis is proposed. Finally, the propagation of tip-vortices is studied. A third LBM-LES simulation with an upstream intensity of $I = 8\%$ is realized. This turbulence intensity has been observed in different areas favourable to tidal stream facilities [3]. Except for the upstream turbulence, this simulation is identical to the simulation with an upstream turbulence intensity of $I = 12.5\%$.

## 4.1. Cross-Comparison

A comparison based on the averaged axial velocity, the average turbulence rate and the average shear stress is proposed here. Quantities from the case at $I = 8\%$ are obtained with the same procedure as for the case at $I = 12.5\%$.

The axial velocity in the wake has a different evolution according to the upstream turbulence rate (Figure 22). The turbine footprint is still visible on the average axial velocity profile of the simulation at $I = 8\%$ at $x = 10D$. The axial velocity $<U_x>$ for the simulation at $I = 8\%$ is always included between the axial velocity of the simulations at $I = 3\%$ and $I = 12.5\%$. The turbulence rate $I_{2D_{bis}}$ in the wake is plotted in Figure 23 and is defined in Equation (18):

$$I_{2D_{bis}} = 100\sqrt{\frac{0.5(<u_x'^2> + <u_y'^2>)}{<U_{x\infty}>^2}} \quad (18)$$

The wake, observed with the quantity $I_{2D_{bis}}$, can be divided into three zones. A first zone from $x = 0D$ to $x = 5D$ where the turbulence rate $I_{2D_{bis}}$ is different for the three cases and increases with the upstream turbulence rate. A second zone between $x = 5D$ and $x = 9D$ is observed. In this zone, the turbulence rates $I_{2D_{bis}}$ of cases at $I = 3\%$ and $I = 8\%$ are equal. The turbulence rate $I_{2D_{bis}}$ from the case $I = 12.5\%$ is greater. A third zone at $x = 10D$ can be identified. In this zone, the turbulence rate $I_{2D_{bis}}$ is the same for all three cases. The shear stress $<u_x'u_y>$ has a different behaviour (Figure 24). From $x = 6D$ to $x = 10D$, the shear stress is very close for all three cases. However, the lack of samples makes observations difficult.

The diameter of the wake can be estimated from the axial velocity profiles (Figure 22). It is calculated assuming that the averaged axial velocity profiles are close to a Gaussian distribution. The diameter of the wake can be estimated as $D_w = (6/1.177)y_{1/2}$, where $y_{1/2}$ is the half-width of the profile at $<U_x> = (<U_{x\infty}> + <U_{xmin}>)/2$. The calculated $D_w$ are summarized in Table 7. The diameter of the wake increases according to the $x$ direction. It is wider for high upstream turbulence rates. Those observations are in agreement with the observations made by [6].

Iso-surfaces of the $\lambda_2$ criterion are shown in Figure 25. The $\lambda_2$ criterion was introduced by Jeong and Hussain [45] to identify the centre of turbulent structures in turbulent flows. More details about this criterion and others can be found in Kolár [46]. Only iso-surfaces included between $R = 0.15$ m and $R = 0.42$ m are shown in Figure 25. The iso-surface of $\lambda_2$ criterion is shown with no upstream turbulence in Figure 25a and with upstream turbulence $I = 3\%$ in Figure 25b. The tip-vortices are visible on both figures. Tip-vortices are coherent turbulent structures with high energy. They are one of the main sources of turbulence in the wake. Studying their propagation is important for wake prediction. The upstream turbulence has an influence on the propagation of those structures. It seems to affect their envelope and the duration during which they are coherent.

Tip-vortices with an upstream turbulence of $I = 8\%$ can be observed in Figure 25c and with $I = 12.5\%$ in Figure 25d. They are highly disrupted by the upstream turbulence. Snapshots of the wake can hardly provide information about the propagation of turbulent structures such as tip-vortices.

Even if only average quantities of the wake are considered, it is difficult to predict a priori the evolution of the wake as a function of the upstream turbulence. The average velocity deficit on the rotation axis seems to have a linear behaviour whereas it is hard to perceive a pattern for average fluctuations and tip-vortices. The next two analyses bring
some elements in order to have a better understanding of the phenomena that take place in
the wake.

<figure><img src="image_27.png" /><figcaption>**Figure 22.** Average axial velocity profiles in the wake of the tidal turbine simulated in LBM-LES at TSR = 3.67 with upstream turbulence intensities of *I* = [3, 8, 12.5]%.</figcaption></figure>
<figure><img src="image_28.png" /><figcaption>**Figure 23.** Average turbulence intensity profiles $I_{2D_{bis}}$ (calculated with $U_{x\infty}$) in the wake of the tidal turbine simulated in LBM-LES at TSR = 3.67 with upstream turbulence intensities of $I = [3, 8, 12.5] \%$.</figcaption></figure>
<figure><img src="image_29.png" /><figcaption>**Figure 24.** Average shear stress profiles $<u_x'u_y>$ in the wake of the tidal turbine simulated in LBM-LES at TSR = 3.67 with upstream turbulence intensities of $I = [3, 8, 12.5]$%.</figcaption></figure>

<table><caption>**Table 7.** Diameter of the wake of tidal turbine at TSR = 3.67 calculated from averaged axial velocity profiles from LBM-LES simulations with upstream turbulence intensities of $I = [3, 8, 12.5]$%.</caption><thead><tr><th>Upstream Turbulence</th><th>Diameter at x = 4D</th><th>Difference</th><th>Diameter at x = 5D</th><th>Difference</th><th>Diameter at x = 7D</th><th>Difference</th></tr></thead><tbody><tr><td>3%</td><td>1.23 m</td><td>\</td><td>1.33 m</td><td>\</td><td>1.50 m</td><td>\</td></tr><tr><td>8%</td><td>1.57 m</td><td>+28%</td><td>1.61 m</td><td>+21%</td><td>2.07 m</td><td>+38%</td></tr><tr><td>12.5%</td><td>1.54 m</td><td>-2%</td><td>2.13 m</td><td>+32%</td><td>2.41 m</td><td>+16%</td></tr></tbody></table>
<figure><img src="image_30.png" /><figcaption>**Figure 25.** Iso-surface of $\lambda_2$ criterion showing tip-vortices of tidal turbine simulated in LBM-LES at TSR = 3.67 with upstream turbulence intensities of $I = [0, 3, 8, 12.5]\%$. (a) $I = 0\%$. (b) $I = 3\%$. (c) $I = 8\%$. (d) $I = 12.5\%$.</figcaption></figure>

## 4.2. *Spectral Analysis*

Experimental spectral data in the wake of the tidal turbine are available in Medina et al. [38]. They used the same turbine as the one simulated here but at TSR = 4.0.

Figure 26a compares LBM results with the experiment made by [38]. Power spectral density (PSD) computed from the LBM simulation with an upstream turbulence rate of 12.5% is close to the experimental one up to 10 Hz frequency. After 10 Hz, the energy of the PSD from the LBM simulation drops below the experimental one. Numerical dissipation and mesh size, $\Delta x = 5.58 \times 10^{-3}$ m, are the main cause of this difference. PSD from the simulation with $I = 3\%$ is further away from the experiment. The PSD from the LBM simulation with $I = 8\%$ is included between the $I = 3\%$ PSD and $I = 12.5\%$ PSD, which could be expected. It can be observed that it is closer to the $I = 3\%$ PSD. A similar observation was made with $I_{2Dbis}$ in the previous section.

Far from the turbine at $x = 8.0D$ (Figure 26b), PSD of the three cases are quite similar. This observation and the observations made in the cross-comparison suggest that the turbulence in the wake converge toward a state independent from the upstream turbulence intensity. It is to be noted that none of the chosen upstream turbulence rates exceed the turbulent rates calculated at $x = 10D$ downstream of the turbine.
Figure 26c shows the PSD calculated from the axial velocity in the wake at $y = R$. A peak is observed on the PSD from the case $I = 3\%$ at the frequency $3f_0$. The frequency $f_0 = 1.335$ Hz is the rotating frequency of the rotor. This peak is not present on the PSD from the cases $I = [8, 12.5]\%$. There are two explanations. Either the tip-vortices have already lost their consistency or they are too much disrupted to be observed at a given location.

The PSD from the simulations has been found to be quite reliable compared to the experimental PSD. Some observations have been made regarding the evolution of the turbulence in the wake. It has been suggested that the turbulence converges toward a state independent from the upstream turbulence. The chosen spectral approach is not suited to observe the evolution of some specific structures like tip-vortices. A spatial approach is proposed hereafter to analyse the propagation of tip-vortices.

<figure><img src="image_31.png" /><figcaption>**Figure 26.** LBM-LES simulation of the tidal turbine at TSR = 3.67. PSD calculated from the axial velocity $U_x$ for $I = [3, 8, 12.5]\%$. Experimental data are from [38] at TSR = 4.0. (a) $x = 4.0D$, $y = z = 0$ m. (b) $x = 8.0D$, $y = z = 0$ m. (c) $x = 2D$, $y = R$ and $z = 0$ m.</figcaption></figure>
4.3. *Spatial Analysis*

Spatial analysis is based on instantaneous pressure maps of the wake. Those maps are
in the xz and xy plans and contain the axis of rotation. There is 13 moments (26 maps) for
the cases I = [0,3]% and 30 moments (60 maps) for the cases I = [8, 12.5]%. Intersections
between plans and tip-vortices are identified with pressure minima. Locations of intersec-
tions are then recorded regardless of the moments or maps they are from. The sampling
rate of the moments is 2 Hz.

The wake is divided according to the stream-wise direction into segments of 0.5 m spaced 0.25 m apart. The number of intersections per segment is $N_i$. If $x_0$ is the lower limit of a segment, the first segment is at $x_0 = 0.5$ m and the last at $x_0 = 2.5$ m. The number of intersections in the first segment is $N_0$ and the ratio $N_t = N_i/N_0$ is the intersections rate. It gives an estimation of the loss of consistency of tip-vortices. Figure 27 shows the evolution of $N_t$ according to the stream-wise direction. Without upstream turbulence, the tip-vortices are still consistent after $x = 3D$. It confirms what is observed in Figure 25a. Even a low turbulence intensity ($I = 3\%$) has an influence over tip-vortices consistency. This was also observed in Figure 25b and in Figure 23. On that last one, the peak of turbulence intensity $I_{2D_{bis}}$ generated by tip-vortices is hardly observable at $x = 3D$. With an upstream turbulence intensity of $I = [8, 12.5]\%$, tip-vortices rapidly lose their consistency (Figure 25d). An analysis based on the height of tip-vortices is now carried out to have more detailed knowledge of the influence of turbulence over tip-vortices propagation.

<figure><img src="image_32.png" /><figcaption>**Figure 27.** LBM-LES simulation of the tidal turbine at TSR = 3.67. Intersections rate $N_t$ at the segments center for $I = [0, 3, 8, 12.5]\%$.</figcaption></figure>

Figure 28 shows the distance $y^*$ from the intersections to the axis $y = z = 0$ m. The height $y^*$ is averaged over 0.3 m segments from $x_0 = 0$ m until $N_t < 30\%$. The height $y^*$ decreases faster as the turbulence increases. A higher upstream turbulence intensity thus increases the trend of tip-vortices to move toward the axis $y = z = 0$ m. The low turbulence rate $I = 3\%$ as a significant impact over $y^*$ and the evolution of $y^*$ for the cases with $I = [8, 12.5]$ are close.

An analysis of tip-vortices propagation has been carried out. Several trends that are
consistent with observations made in the previous section have been highlighted. Since one
of the main sources of turbulence in the close wake is tip-vortices, additional comments
about the results of the previous section can be made. In Figure 20, a difference between
the simulation and the experiment is observed at $x = 1.2D$. The turbulence intensity from
the simulation is too low compared to the experiment. By considering observations made
on tip-vortices, a third source of error can be suggested. Tip-vortices may move too quickly
toward the axis $y = z = 0$ m, modifying the shape of turbulence intensity profiles. It is
known that the integral length scale for the case $I = 12.5\%$ is larger in the simulation than in
the experiment. The integral length scale of the turbulence may influence the propagation of tip-vortices, and modify the turbulence intensity profiles.

<figure><img src="image_33.png" /><figcaption>**Figure 28.** LBM-LES simulation of the tidal turbine at TSR = 3.67. Average height of tip-vortices intersections for $I = [0, 3, 8, 12.5]\%$.</figcaption></figure>

**5. Conclusions**

1. The modelling of a tidal turbine in a turbulence flow with an open-source LBM solver has been performed. To achieve this, an IBM boundary condition coupled to a wall model and a synthetic turbulence generation algorithm were implemented in the code.

2. A realistic upstream turbulent condition has been validated for an LBM-LES tidal turbine modelling. For low turbulence rates, it has been observed that the choice of parameter is straightforward. The characteristics of the synthetic turbulence match experimental ones. It has been found that for a high turbulence rate, compromises have to be made in order to keep the turbulence rate high enough. As a consequence, the integral length scale of the synthetic turbulence is greater than the experimental one.

3. Blade resolved IBM-LBM-LES simulations of a three-bladed HATT have been carried out for two turbulent rates. Three meshes have been tested to estimate the mesh sensitivity of the approach. It has been observed that the mesh convergence in the thrust coefficient is quickly reached. Differences in power coefficient are still observed for the finest mesh. For the high turbulence rate case, those differences may come from the upstream synthetic turbulence.

4. Average wake quantities from LBM-LES simulations have been compared with experimental ones. It has been observed that refining the mesh has little influence over the average wake. It slightly improves the predicted average axial velocity in the close wake. Elsewhere in the wake, average quantities computed with the LBM-LES are in accordance with experimental data. The influence of turbulence on the turbine behaviour is well predicted by the model.

5. A third case has been simulated with an upstream turbulence rate of $I = 8\%$. A cross-comparison has been carried out between wakes of the three different cases $I = [3, 8, 12.5]\%$. The evolution of the averaged velocity deficit is quite intuitive. This is not the case for other quantities. Turbulence intensity and shear stress from cases $I = 3\%$ and $I = 8\%$ are close whereas for tip-vortices evolution, cases $I = 8\%$ and $I = 12.5\%$ are close. Based on turbulence intensity, shear stress and PSD, far wakes of the three cases converge toward similar states.

6. It has been observed with iso-surfaces of $\lambda_2$ criterion that typical wake structures like tip-vortices are greatly influenced by ambient turbulence. An analysis of the evolution of tip-vortices in the wake has been proposed. The envelope of tip-vortices is shortened in the x direction and shrunk toward the axis with the increase in turbulence rate. The close wake is thus greatly influenced by the upstream turbulence.
7. Improvements to the model or further investigations can be suggested. Studying the influence of the integral length scale on the wake is an interesting lead. It would confirm or invalidate some of the hypotheses made. Finally, improvements of the boundary layer model are one of the top priorities. It would allow for proper force computation. If forces from the simulations are reliable enough, it is then possible to estimate the turbulence-induced fatigue. Such information is substantial for turbine optimization.

**Author Contributions:** Software, M.G., P.M. and E.P.; validation, M.G.; writing—original draft preparation, M.G. and S.S.G.; writing—review and editing, M.G., S.S.G., P.M., J.C.P., E.P. and Y.M.; supervision, S.S.G.; project administration, S.S.G.; funding acquisition, S.S.G. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported by Agence Nationale de la Recherche et de la Technologie (ANRT) CIFRE contract 2015/1194, with the collaboration of SIREHNA.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Acknowledgments:** This work was supported by ANRT CIFRE contract 2015/1194, with the collaboration of SIREHNA. S.S.G. acknowledges the financial support of the Tidal Stream Industry Energiser (TIGER) project, co-financed by the European Regional Development Fund through the Interreg France (Channel) England Programme. Part of the computational means was founded by Manche County Council. Most of this work was performed using the computing resources of CRIANN (Normandy, France). The authors would like to thank IFREMER for the provided data and blade geometry. The authors would also like to thank the people that carried out the Palabos project.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results

**Abbreviations**

The following abbreviations are used in this manuscript:

<table><tr><td>CFD</td><td>Computational Fluid Dynamics</td></tr><tr><td>DFSEM</td><td>Divergence Free Synthetic Eddy Method</td></tr><tr><td>DNS</td><td>Direct Numerical Simulation</td></tr><tr><td>HATT</td><td>Horizontal Axis Tidal Turbine</td></tr><tr><td>IBM</td><td>Immersed Boundary Method</td></tr><tr><td>LBM</td><td>Lattice Boltzmann Method</td></tr><tr><td>LES</td><td>Large Eddy Simulation</td></tr><tr><td>PSD</td><td>Power Spectral Density</td></tr><tr><td>SEM</td><td>Synthetic Eddy Method</td></tr><tr><td>TSR</td><td>Tip Speed Ratio</td></tr></table>
<figure><img src="image_34.png" /><figcaption>Appendix A. All Profiles for I = 3%</figcaption></figure>

**Figure A1.** Average axial velocity profiles in the wake of the tidal turbine at TSR = 3.67 and I = 3%. Results are from experiment [6] and from LBM-LES simulations with mesh Nos. 1 and 2. Complete profile set.
<figure><img src="image_35.png" /><figcaption>**Figure A2.** Average turbulence intensity profiles $I_{2D}$ in the wake of the tidal turbine at TSR = 3.67 and $I = 3\%$. Results are from experiment [6] and from LBM-LES simulations with mesh Nos. 1 and 2. Complete profile set.</figcaption></figure>
<figure><img src="image_36.png" /><figcaption>Appendix B. All Profiles for $I = 12.5\%$</figcaption></figure>

**Figure A3.** Average axial velocity profiles in the wake of the tidal turbine at TSR = 3.67 and $I = 12.5\%$. Results are from experiment [6] and from LBM-LES simulations with mesh Nos. 1 and 2. Complete profile set.
<figure><img src="image_37.png" /><figcaption>**Figure A4.** Average turbulence intensity profiles *I*<sub>2D</sub> in the wake of the tidal turbine at *TSR* = 3.67 and *I* = 12.5%. Results are from experiment [6] and from LBM-LES simulations with mesh Nos. 1 and 2. Complete profile set.</figcaption></figure>

References

1. Hardisty, J. *The Analysis of Tidal Stream Power*; John Wiley & Sons Ltd.: New York, USA, 2009.

2. Clark, T. Turbulence in Marine Environments (TiME): A framework for understanding turbulence and its effects on tidal devices. In Proceedings of the 11th European Wave and Tidal Energy Conference, Nantes, France, 6–11 September 2015.

3. Thomson, J.; Polagye, B.; Durgesh, V.; Richmond, M.C. Measurements of turbulence at two tidal energy sites in Puget Sound, WA. *IEEE J. Ocean. Eng.* **2012**, *37*, 363–374. [CrossRef]

4. Milne, I.A.; Sharma, R.N.; Flay, R.G.J.; Bickerton, S. Characteristics of the turbulence in the flow at a tidal stream power site. *R. Soc.* **2013**, *371*, 196. [CrossRef]
5. Chen, Y.; Lin, B.; Lin, J.; Wang, S. Experimental study of wake structure behind a horizontal axis tidal stream turbine. *Appl. Energy* **2017**, *196*, 82–96. [CrossRef]

6. Mycek, P.; Gaurier, B.; Germain, G.; Pinon, G.; Rivaolen, E. Experimental study of the turbulence intensity effects on marine current turbines behaviour. Part I: One single turbine. *Renew. Energy* **2014**, *66*, 729–746. [CrossRef]

7. Blackmore, T.; Gaurier, B.; Myers, L.; Germain, G.; Bahaj, A.S. The effect of freestream turbulence on tidal turbines. In Proceedings of the 11th European Wave and Tidal Energy Conference, Nantes, France, 6–11 September 2015.

8. Batten, W.M.J.; Harrison, M.E.; Bahaj, A.S. Accuracy of the actuator disc-RANS approach for predicting the performance and wake of tidal turbines. *Philos. Trans. R. Soc. A* **2013**, *371*, 20120293. [CrossRef] [PubMed]

9. Elie, B.; Oger, G.; P.-E., G.; Alessandrini, B. Simulation of horizontal axis tidal turbine wake using a Weakly-Compressible Cartesian Hydrodynamic solver with local mesh refinement. *Renew. Energy* **2017**, *108*, 336–354. [CrossRef]

10. Grondeau, M.; Guillou, S.; Mercier, P.; Poizot, E. Wake of a ducted vertical axis tidal turbine in turbulent flows, LBM actuator-line approach. *Energies* **2019**, *12*, 4273. [CrossRef]

11. Ouro, P.; Stoesser, T. Wake generated downstream of a vertical axis tidal turbine. In Proceedings of the 12th European Wave and Tidal Energy Conference, Cork, Ireland, 27 August–1 September 2017.

12. Ebdon, T.; O'Doherty, D.M.; O'Doherty, T.; Mason-Jones, A. Modelling the effect of turbulence length scale on tidal turbine wakes using advanced turbulence models. In Proceedings of the 12th European Wave and Tidal Energy Conference, Cork, Ireland, 27 August –1 September 2017.

13. Ahmed, U.; Apsley, D.D.; Afgan, I.; Stallard, T.; Stansby, P.K. Fluctuating loads on a tidal turbine due to velocity shear and turbulence: comparison of CFD with fields data. *Renew. Energy* **2017**, *112*, 235–246. [CrossRef]

14. Guo, Z.; Shu, C. *Lattice Boltzmann Method and Its Application in Engineering*; World Scientific: Singapore, 2013.

15. Marié, S.; Ricot, D.; sagaut, P. Comparison between lattice Boltzmann Method and Navier—Stokes high order schemes for computational aeroacoustics. *J. Comput. Phys.* **2009**, *228*, 1056–1070. [CrossRef]

16. Deiterding, R.; Wood, S. Predictive wind turbine simulation with an adaptive lattice Boltzmann method for moving boundaries. *J. Phys. Conf. Ser.* **2016**, 753, 082005. [CrossRef]

17. Mercier, P.; Grondeau, M.; Guillou, S.; Thiebot, J.; Poizot, E. Numerical study of the turbulent eddies generated by the seabed roughness. Case study at a tidal power site. *Appl. Ocean Res. 2020*, *97*, 102082. [CrossRef]

18. Mercier, P.; Guillou, S. The impact of the seabed morphology on turbulence generation in a strong tidal. *Phys. Fluids 2020*, *33*, 055125. [CrossRef]

19. Grondeau, M.; Poirier, J.C.; Guillou, S.; Méar, Y.; Mercier, P.; Poizot, E. Modelling the wake of a tidal turbine with upstream turbulence: LBM-LES versus Navier-Stokes LES. *Int. Mar. Energy J.* **2020**, *3*, 83–89. [CrossRef]

20. Ouro, P.; Harrold, M.; Stoesser, T.; Bromley, P. Hydrodynamic loadings on a horizontal axis tidal turbine prototype. *J. Fluids Struct.* **2017**, *71*, 78–95. [CrossRef]

21. Lagrava, S.; Daniel, W.; Malaspinas, O.; Latt, J.; Chopard, B. Advances in multi-domain lattice Boltzmann grid refinement. *J. Comput. Phys.* **2012**, *231*, 4808–4822. [CrossRef]

22. Polleto, R.; Craft, T.; Revell, A. A new divergence free synthetic eddy method for the reproduction of inlet flow conditions for LES. *Flow Turbul. Combust.* **2013**, *91*, 519–539. [CrossRef]

23. Malaspinas, O. Increasing stability and accuracy of the lattice Boltzmann scheme: recursivity and regularization. *arXiv 2015*, arxiv:1505.06900v1.

24. Skordos, P.A. Initial and boundary conditions for the lattice Boltzmann method. *Phys. Rev. E* **1993**, *48*, 4823–4842. [CrossRef]

25. Latt, J.; Chopard, B.; Malaspinas, O.; Deville, M.; Michler, A. Straight velocity boundaries in the lattice Boltzmann method. *Phys. Rev. E* **2008**, *77*, 056703. [CrossRef]

26. Smagorinsky, J. General circulation experiments with the primitive equations. *Mon. Weather Rev.* **1963**, *91*, 99–163. [CrossRef]

27. Malaspinas, O.; Sagaut, P. Consistent subgrid scale modelling for lattice Boltzmann methods. *J. Fluid Mech.* **2012**, *100*, 514–542. [CrossRef]

28. Peskin, C.S. Numerical Analysis of Blood Flow in the Heart. *J. Comput. Phys.* **1977**, *25*, 220–252. [CrossRef]

29. Feng, Z.G.; Michaelides, E.E. The immersed boundary-lattice Boltzmann method for solving fluid-particles interaction problems. *J. Comput. Phys.* **2003**, *195*, 602–628. [CrossRef]

30. Guo, Z.; Zheng, C.; Shi, B. Non-equilibrium extrapolation method for velocity and pressure boundary conditions in the lattice Boltzmann method. *Chin. Phys.* **2002**, *11*, 366–374.

31. Yu, Z.; Shao, X. A Direct-Forcing Fictitious Domain Method for Particulate Flows. *J. Comput. Phys.* **2007**, *227*, 292–314. [CrossRef]

32. Pope, S.B. *Turbulent Flows*; Cambridge University Press: Cambridge, UK, 2000.

33. Ruffin, S.M.; Lee, J.A. Adaptation of a k-epsilon model to a cartesian grid based methodology. *Int. J. Math. Models Methods Appl. Sci.* **2009**, *3*, 238–245.

34. Spalding, D. A single formula for the law of the wall. *J. Appl. Mech.* **1961**, *28*, 455–458. [CrossRef]

35. Malaspinas, O.; Sagaut, P. Wall model for large-eddy simulation based on the lattice Boltzmann method. *J. Comput. Phys.* **2014**, *275*, 25–40. [CrossRef]

36. Grondeau, M.; Mercier, P.; Guillou, S.; Poirier, J.; Poizot, E.; Mear, Y. Numerical simulation of a tidal turbine model in a turbulent flow with the Lattice Boltzmann Method. In Proceedings of the 12th European Wave and Tidal Energy Conference, Cork, Ireland, 27 August –1 September 2017.
37. Jarrin, N.; Prosser, R.; Uribe, J.C.; Benhamadouche, S.; Laurence, D. Reconstruction of turbulent fluctuations for hybrid RANS/LES simulations using a Synthetic-Eddy-Method. *Int. J. Heat Fluid Flow* **2009**, *30*, 435–442. [CrossRef]

38. Medina, O.D.; Schmitt, F.G.; Calif, R.; Germain, G.; Gaurier, B. Turbulence analysis and multiscale correlations between synchronised flow velocity and marine turbine power production. *Renew. Energy* **2017**, *112*, 314–327. [CrossRef]

39. Medina, O.D.; Schmitt, F.G.; Germain, G.; Gaurier, B. Correlation between synchronised power and flow measurements, a way to characterize turbulence effects on a marine current turbine. In Proceedings of the European Wave and Tidal Energy Conference (EWTEC), Nantes, France, 6–11 September 2015.

40. Liu, C.; Hu, C. An actuator line-immersed boundary method for simulation of multiple tidal turbines. *Renew. Energy* **2019**, *136*, 473–490. [CrossRef]

41. Xu, H.; Sagaut, P. Analysis of the absorbing layers for the weakly-compressible lattice Boltzmann schemes. *arXiv* **2012**, arxiv:1203.6350.

42. Blackmore, T.; Myers, L.E.; Bahaj, A.S. Effects of turbulence on tidal turbines: Implications to performance, blade loads, and condition monitoring. *Int. J. Mar. Energy* **2016**, *14*, 1–26. [CrossRef]

43. Shah, N.; Animasaun, I.; Ibraheem, R.; Babatunde, H.; Sandeep, N.; Pop, I. Scrutinization of the effects of Grashof number on the flow of different fluids driven by convection over various surfaces. *J. Mol. Liq.* **2018**, *249*, 980–990. [CrossRef]

44. Animasaun, I.; Shah, N.; Wakif, A.; Mahanthesh, B.; Sivaraj, R.; Koriko, O. *Ratio of Momentum Diffusivity to Thermal Diffusivity: Introduction, Meta-analysis, and Scrutinization.*; Taylor & Francis Ltd.: London, UK, 2022.

45. Jeong, J.; Hussain, F. On the identification of vortex. *J. Fluid Mech.* **1995**, *285*, 69–94. [CrossRef]

46. Kolár, V. Vortex identification: New requirements and limitations. *Int. J. Heat Fluid Flow* **2007**, *28*, 638–652. [CrossRef]