**1. Introduction**

When a surface is submerged in natural water, it is submitted to colonisation by numerous species. The accumulations of biological organisms on this surface is called biofouling. The first step of this phenomenon is the colonisation by micro-organisms like bacteria or micro-algae into a thin layer a few millimetres thick (biofilm) [1]. This thin film serves as a base for larger sessile species such as barnacles and mussels [2], which often organise in habitats allowing the arrival of mobile organisms (shrimps, crabs, . . . ). The speed, the kind and the arrangement of the colonisation depend on many factors [3]: depth, salinity, kind of surface, region, temperature, etc. Larvae are also able to select the substrate suitable for their development depending on the water streaming and chemical properties. Moreover, specific conditions can lead to particularities among the species. For example, in the Alderney Race located between France and United Kingdom, where the tidal stream is extremely energetic, the species evolving in such conditions are mainly smaller, with more developed fixing organs and a smoother surface [4]. Some studies show that biofouling also occurs on moving solids like ship propellers or tidal turbines. The biofouling is the origin of the artificial reef effect created by off-shore artificial constructions such as the wind turbine farms described in [5]. However, despite this positive effect, biofouling is a real challenge for the marine renewable implantation in sea water. Ref. [6] shows its negative impact on boat hulls. The experimental study shows a higher resistance on heterogeneous rough parts of the hull, which is very dependent on the position of the roughnesses [7]. More generally, studies of the effects of small roughnesses on the flow along a flat plate are well-understood in both air [8] and water [9]. Small structures promote

**Citation:** Robin, I.; Bennis, A.-C.; Dauvin, J.-C. 3D Numerical Study of the Impact of Macro-Roughnesses on a Tidal Turbine, on Its Performance and Hydrodynamic Wake. *J. Mar. Sci. Eng.* **2021**, *9*, 1288. https://doi.org/ 10.3390/jmse9111288

Academic Editor: Eugen Rusu

Received: 14 October 2021 Accepted: 9 November 2021 Published: 18 November 2021

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

**Copyright:** © 2021 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/).

the boundary layer transition by making the boundary layer transition move upstream and also increase the rate of turbulence. With the development of renewable marine energies and some special cases in aviation, studies are aiming towards understanding the effects of various roughnesses on wing profiles [10]. Ref. [11] looks into the effect of ice accretion located on the leading edge and finds that the aerodynamic performance of the profile is significantly deteriorated by ice. Ref. [12] experimentally studied the effects of numerous cones distributed homogeneously over the surface of a blade in a wind tunnel. The study is carried out in stationary flow with solid cones ranging from 0.0035 c (chord) to 0.0285 c in height. Results show a strong increase in the drag. Ref. [13] investigates the aerodynamic impact of an isolated barnacle-shaped excrescence on a blade in a wind tunnel. The barnacle is 0.02 c high and is located at 60% of the chord. Three angles for the flow incidence are studied: 5°, 10° and 15° under two different operating modes: a stationary or oscillating blade. On the motionless blade, the barnacle has no impact on lift but considerably increases drag. The pressure coefficient field is very different between the 5° and 15° cases in the motionless case but seems very similar once the blade is set in motion.

Nevertheless, the marine environment at tidal sites is not easy to study in realistic conditions. So, the main parts of these studies are conducted in ideal conditions, easily replicable experimentally. In fact, marine currents are really strong and numerical approach is more suitable to predict biofouling in more realistic conditions. Most papers aim at developing models for small roughnesses. However with the development of CFD, more and more 2D studies are emerging: Ref. [14] investigates the effects of biofouling on a NACA0018 profile using numerous turbulence models (Reynolds Averaged Navier– Stokes (RANS): k-, k-*ω* SST (Shear Stress Transport) and Large Eddy Simulation (LES): Smagorinsky). RANS k-*ω*-SST and LES-Smagorinsky gave close results for the blade performances but RANS models are less efficient than LES models to compute the vortexes in the wake. The results show that if the biofouling species are placed on the first half of the blade chord, their effects are maximised. Conversely, when biofouling is in the second half, the effects are reduced.This phenomenon is explained by the isolated roughnesses that create an early detachment of the boundary layer. Nevertheless, biofouling generates three-dimensional effects in the flow [15,16]. A 2D study is therefore not sufficient to fully understand such effects.

To our knowledge, the only work relating to the impact of biofouling on the performance of tidal turbines applied to a complete turbine is [17], where biofouling is a surface finish applied with a wall function on the whole turbine. In this paper, we propose to investigate the impact of a realistic fouling on an entire turbine. After a short introduction, Section 2 explains the methodology and methods used to carry out the numerical tests. Section 3 is the validation of the model and its comparison with the experimental data of [13]. Section 4 shows our results that are discussed. Section 4 draws the conclusion.

#### **2. Materials and Methods**

The numerical simulations are divided into two parts. The first part is a 3D model of a single blade with only one barnacle. This aims to validate the numerical model using the experimental results. The second part relates to the modelling of a complete biofouled rotor. That allows us to study the impact of a realistic implantation of sessile species on the performance and the wake of the tidal turbine. Unlike the work of [17], biofouling is explicitly represented in the mesh by integrating the barnacles into the 3D structure and not by using a roughness model applied to the turbine.

#### *2.1. Experimental Setups*

Experiments were carried out by [13,18] and their data were used for comparison with the numerical results.

The results of the **motionless blade experiment** chosen for the validation of the numerical model comes from [13]. Their experimental method is described hereafter. The tests were carried out in the Handley–Page wind tunnel with dimensions of 1.61 m high ×

2.13 m wide × 2.74 m long. The air flow velocity was 45 m·s−<sup>1</sup> with 2.5% turbulence. The blade is made from a NACA 63-619 foil of 0.55 m chord. It passes completely through the height of the wind tunnel. The barnacle is represented by a solid cone with a radius of 20 mm at its base and 10 mm at its top and a height of 11 mm, dimensions nearby to the *Balanus crenatus* found in the Alderney Race [4]. A total of 25 pressure orifices are positioned on and around the barnacle in order to follow the evolution of the pressure field. Three angles of attack are studied: 5°, 10° and 15°. Sensor HDI series gauge sensors measured the dynamic stall while pressure transducers were used to sample the pressure evolution along the blade. The experimental setup also allowed blade oscillation, but unsteady tests were not used here.

The **full rotor simulation** is built according to [18] and uses the IFREMER-LOMC1 horizontal axis turbine. The water velocity was fixed to 0.8 m·s−<sup>1</sup> with 3% of turbulence intensity. The rotor rotation was forced by a motor to fix the Tip Speed Ratio (TSR) to 4. TSR is defined as follows :

$$
\lambda = \frac{|\Omega\_{\mathbf{x}}|R}{\mathcal{U}\_{\infty}},
\tag{1}
$$

where Ω*<sup>x</sup>* is the angular velocity, *R* is the radius of the tidal turbine and *U*<sup>∞</sup> is the inlet flow velocity. The rotor characteristics are presented in Table 1.


**Table 1.** General characteristics of the IFREMER-LOMC turbine

Torque sensors were used to measure the power and drag coefficients of the entire structure, including the hub. Laser Doppler Velocimeters (LDV) were used to monitor the wake and vortexes.

#### *2.2. Governing Equations*

The numerical model is used for validation (in air) and investigation (in water) cases. For the numerical simulations, the following hypothesis are considered:


The 3D Navier–Stokes equations are suitable to solve the fluid motion under the incompressible assumption (with *i* = 1, 2, 3 representing the 3 directions in a Cartesian framework):

$$\frac{\partial \mu\_i}{\partial x\_i} = 0,\tag{2}$$

$$\frac{\partial \mu\_i}{\partial t} + \frac{\partial u\_i u\_j}{\partial x\_j} = -\frac{1}{\rho} \frac{\partial p}{\partial x\_i} + \frac{1}{\rho} f\_i + \nu \frac{\partial^2 u\_i}{\partial x\_j \partial x\_j} \tag{3}$$

where *ui* is the fluid velocity in the *<sup>i</sup>*-direction, *<sup>t</sup>* is the time, *<sup>ρ</sup>* is the fluid density (kg·m<sup>−</sup>3), *p* is the pressure, *fi* represents the volumetric forces in the *i*-direction, and *ν* is the kinematic viscosity.

The RANS k-*ω* SST and LES Smagorinsky turbulence models give close results for the calculation of blade forces. However, the RANS approaches, including the k-*ω* SST model, known to be cheaper, do not give accurate results for the calculation of the wake. On the other hand, LES models like Smagorinsky offer better results for the generation and development of vortices in the wake but the computational cost is high. Thus, the both turbulence models are compared in this work.

#### 2.2.1. Reynolds-Averaged Navier-Stokes Turbulence Model

In RANS models, the velocity *ui* is decomposed into an averaged part (*ui*) and a fluctuating part (*u i* ) such as:

$$
u\_i = \overline{u\_i} + \boldsymbol{u\_i}.\tag{4}$$

The low-frequency component is obtained by applying the Reynolds average (·) to the instantaneous velocity. This average is also applied to the pressure leading to the same decomposition.

In this framework, the Navier–Stokes equations are :

$$\frac{\partial \overline{\mu\_i}}{\partial \mathbf{x}\_i} = \mathbf{0},\tag{5}$$

$$\frac{\partial \overline{u\_i}}{\partial t} + \frac{\partial \overline{u\_i} \, \overline{u\_j}}{\partial x\_j} = -\frac{1}{\rho} \frac{\partial \overline{p}}{\partial x\_i} + \frac{1}{\rho} \overline{f\_i} + \nu \frac{\partial^2 \overline{u\_i}}{\partial x\_j \partial x\_j} - \frac{\partial u\_i' u\_j'}{\partial x\_j},\tag{6}$$

where *p* is the mean pressure.

The mean value *ui* is considered as varying slightly in time compared to the variation of fluctuation whereas the mean value of *u <sup>i</sup>* is zero. In order to solve Equations (5) and (6), the knowledge of the Reynolds stress tensor *u i u <sup>j</sup>* is necessary. After approximating this term using a turbulence viscosity depending on both turbulence kinetic energy (*k*) and specific dissipation rate (*ω*), the evolution equations for *k* and *ω* need to be solved to determine these quantities [19]. They are initialised as :

$$k = \frac{3}{2} (I |\mathcal{U}\_{\infty}|)^2,\tag{7}$$

where *I* is the turbulence intensity and *U*∞ is the reference velocity (undisturbed velocity),

$$
\omega = \frac{k^{0.5}}{0.009^{0.25}L} \,\text{\,\,\,}\tag{8}
$$

where *L* is a reference length scale equal to the chord of the profile (c) for the motionless blade simulation and to the rotor radius (R) for the full rotor simulation.

#### 2.2.2. Large Eddy Simulation Turbulence Model

For LES, flow characteristics are separated into two parts according to the turbulent scales by applying a mathematical filter (˜·). *u*˜*<sup>i</sup>* is composed of the large eddies whose size is greater than the size of the filter. *u*∗ carries the smaller eddies with a size inferior to the filter size. In the Smagorinsky turbulence model, the filter size is correlated with the mesh size and *u*˜*<sup>i</sup>* is solved explicitly by solving :

$$\frac{\partial \overline{u}\_{i}}{\partial t} + \frac{\partial \overline{u}\_{i}\overline{u}\_{j}}{\partial x\_{j}} = \frac{1}{\rho} \left( -\frac{\partial \overline{p}}{\partial x\_{i}} + \mathcal{f}\_{i} \right) + (\nu + \nu\_{s\xi s}) \frac{\partial^{2} \overline{u}\_{i}}{\partial x\_{j}x\_{j}} - \frac{\partial \tau\_{ij}}{\partial x\_{i}},\tag{9}$$

where *p*˜ and ˜ *fi* are the filtered pressure and volumetric forces, respectively. *νsgs* is the turbulent eddy viscosity.In the original Smagorinsky model, *νsgs* is computed as :

$$\nu\_{\mathfrak{S}^{\mathfrak{S}}} = \left(\mathbb{C}\_{k}\Delta\right)^{2} \sqrt{2S\_{\vec{i}\vec{j}}S\_{\vec{i}\mathfrak{j}}} \tag{10}$$

where Δ is the width of the filter and *Sij* is the resolved-scale strain rate tensor. However, the Sub Grid Scale Kinetic Energy (SGS TKE) variation of the Smagorinsky model is used here. This choice is made because of the ability of the SGS TKE to evaluate the forces and model the flows in the near walls. The classical Smagorisky model does not allow modelling of the viscous sub-layer of the boundary layer without drastically increasing the number of cells [20]. Thus, *νsgs* is written:

$$\upsilon\_{s\\$^\ast} = \mathbb{C}\_c \Delta k\_1^{0.5} \, , \tag{11}$$

where *Ce* is a constant of the SGS Kinetic Energy model constant and *k*<sup>1</sup> is the turbulent kinetic energy computed according to:

$$\frac{\partial k\_1}{\partial t} + \tilde{u}\_l \frac{k\_1}{\partial x\_i} = c\_\varepsilon \Delta\_x \sqrt{k\_1} \tilde{S}\_{i\bar{j}} \tilde{S}\_{i\bar{j}} - \mathbb{C}\_\varepsilon \frac{k\_1^\frac{3}{2}}{\Delta\_x} - \frac{1}{\partial x\_{\bar{j}}} \left[ \left( \nu + \frac{c\_\varepsilon}{\sqrt{k\_1}} \Delta \sqrt{k\_1} \right) \frac{\partial k\_1}{\partial x\_i} \right],\tag{12}$$
