Embedding of a Blade-Element Analytical Model into the SHYFEM Marine Circulation Code to Predict the Performance of Cross-Flow Turbines

: Our aim was to embed a 2D analytical model of a cross-ﬂow tidal turbine inside the open-source SHYFEM marine circulation code. Other studies on the environmental impact of Tidal Energy Converters use marine circulation codes with simpliﬁed approaches: performance coe ﬃ cients are ﬁxed a priori regardless of the operating conditions and turbine geometrical parameters, and usually, the computational grid is so coarse that the device occupies one or few cells. In this work, a hybrid analytical computational ﬂuid dynamic model based on Blade Element Momentum theory is implemented: since the turbine blades are not present in the grid, the ﬂow is slowed down by means of bottom frictions applied to the seabed corresponding to forces equal and opposite to those that the blades would experience during their rotation. This simpliﬁed approach allowed reproducing the turbine behavior for both mechanical power generation and the turbine e ﬀ ect on the surrounding ﬂow ﬁeld. Moreover, the model was able to predict the interaction between the turbines belonging to a small cluster with hugely shorter calculation time compared to pure Computational Fluid Dynamics.


Introduction
The Earth is called "the blue planet"; in fact, water occupies a large part of the Earth's surface, and moreover, the oceans constitute also immense reserves of energy, which can be exploited only in particular conditions. Close to the coasts, the main resource is represented by tidal currents, which originated from the succession of high and low tides and, for this reason, they are characterized by periodic variations in magnitude and direction that are predictable and quantifiable. The sites of greatest interest in Europe are mainly located in the English Channel [1,2], in the Irish Sea [3], and around Scotland [4], where the proximity to the Atlantic Ocean and the maze of canals, straits, and promontories make the tidal currents particularly intense.
Tidal Energy Converters (TEC) can be classified into seven macro categories, as indicated by the EMEC (European Marine Energy Center) [5]: Horizontal Axis Turbines (HATs); Cross-Flow Turbines (CFTs), also named vertical axis turbines; oscillating hydrofoil, in which a hydrofoil is connected to an oscillating arm; confined turbines in a conduit that concentrates the flow on the device; Archimedes screws that are helical-shaped devices rotated by the water flow; and even tidal kites attached to the study is 2D, and the easiest method to subtract momentum from the flow field-that is, the bottom friction-is adopted. Our objective is to answer these questions: "Is it possible to describe the blade torque evolution during a revolution with grid local refinements that are comparable to the finest ones nowadays adopted in AD marine applications? Is our approach adequate to capture all the turbine main effects on the surrounding flow field, therefore not only the wake characteristics but also the interactions between closely spaced turbines that are responsible for a power increasing in CFTs?"

Methodology
This chapter describes the methodology underlying the work, starting with the validation of the pure CFD model, which is then used to validate the hybrid models instead of experimental data. Subsequently, the foundations of the hybrid BEM-CFD approach are recalled, and then, the grid sensitivity analysis is done with the hybrid model implemented for the software ANSYS-Fluent [39]: this is useful in order to individuate a grid local refinement suitable for marine codes. Finally, the BEM-SHYFEM formulation is presented with a focus on the bottom friction treatment.

Preliminary Validation of the ANSYS CFD Model
The results achieved with the hybrid BEM-CFD models (based on ANSYS or SHYFEM) will be compared with pure CFD simulations in the next sections, since no experimental data are available for the hydrokinetic turbine of interest. Of course, this method requires a preliminary validation of the pure ANSYS-Fluent model with respect to other experimental data from the literature in order to verify the reliability of the setting choices summarized in Table 1 and adopted for both the ANSYS approaches: full CFD and hybrid. The preliminary validation task was done with respect to the data achieved by Bravo et al. [40] in a wind tunnel for a CFT equipped with 3 straight blades, in which the airfoil is the symmetrical NACA0015. These experimental data were chosen due to the relatively high Reynolds (with respect to typical wind tunnel and water flume tests) and low blockage ratio (thanks to the very large section of the wind tunnel). The turbine's performance is usually given by the Power Coefficient, C P , which is defined as: where P is power generated, ρ is the density, U is the free stream flow velocity, and A is the cross-sectional area of the rotor that, in case of CFTs in 2D domains, corresponds to the turbine diameter, d. The curve C P -TSR describes the turbine performance as a function of the Tip Speed Ratio, TSR, which is defined as: where R is the turbine radius and Ω is the rotational speed. It is well known that C P predicted from 2D simulations is overestimated since important 3D effects-first of all, the blade tip losses-cannot be taken into account. In particular, due to the moderate value of the chord-based aspect ratio of this turbine (blade length/chord = 7.5), significant tip losses are expected. Since these 2D CFD simulations cannot directly predict tip losses, a correction factor must be used to reduce C P afterwards in order to allow a meaningful comparison with experimental data. On the basis of the study [41], a plausible value of tip losses for a CFT with this aspect ratio could be ≈25%. Figure 1 depicts the results of CFD (having already subtracted 25% from the 2D original predictions) compared to the experimental results obtained for a flow speed of 10 m/s [40]. It can be seen that the trend, the optimal TSR, and the maximum C P are well matched, with just a small discrepancy in the optimal TSR. In order to match the experimental data at the low TSRs, particular attention was paid to the quality of the calculation grid ( Figure 2); in fact, around each airfoil, many layers of regular quadrilateral cells were created of increasing height with exponential law starting from the blade wall. The following parameters have been kept for all the 2D grids with blades adopted in this study: 710 cells around the blade profile; c × 4 × 10 −6 (where c is the airfoil chord) for the height of the first layer of cells at the blade wall; 1.06 for the cell height growing ratio; 110 layers around the blade profile. It should be observed that the cell smallest height depends on the particular chord value. In this way, around each blade, a very fine and high-quality grid is obtained for a distance of about 1.5 × c from the airfoil/hydrofoil. It could be objected that the experiments adopted in the present work to validate the pure CFD model refer to a turbine [40] with a moderately lower Reynolds and a significantly greater solidity than that of interest. However, we are confident that since the pure CFD model provides realistic results even under more critical conditions for flow separation phenomena (as characterized by a lower Reynolds and higher solidity) than the conditions to which the hybrid models are applied in the present study, it can be considered reliable. where R is the turbine radius and Ω is the rotational speed. It is well known that CP predicted from 2D simulations is overestimated since important 3D effects-first of all, the blade tip losses-cannot be taken into account. In particular, due to the moderate value of the chord-based aspect ratio of this turbine (blade length/chord = 7.5), significant tip losses are expected. Since these 2D CFD simulations cannot directly predict tip losses, a correction factor must be used to reduce CP afterwards in order to allow a meaningful comparison with experimental data. On the basis of the study [41], a plausible value of tip losses for a CFT with this aspect ratio could be ≈25%. Figure 1 depicts the results of CFD (having already subtracted 25% from the 2D original predictions) compared to the experimental results obtained for a flow speed of 10 m/s [40]. It can be seen that the trend, the optimal TSR, and the maximum CP are well matched, with just a small discrepancy in the optimal TSR. In order to match the experimental data at the low TSRs, particular attention was paid to the quality of the calculation grid ( Figure 2); in fact, around each airfoil, many layers of regular quadrilateral cells were created of increasing height with exponential law starting from the blade wall. The following parameters have been kept for all the 2D grids with blades adopted in this study: 710 cells around the blade profile; c × 4 × 10 −6 (where c is the airfoil chord) for the height of the first layer of cells at the blade wall; 1.06 for the cell height growing ratio; 110 layers around the blade profile. It should be observed that the cell smallest height depends on the particular chord value. In this way, around each blade, a very fine and high-quality grid is obtained for a distance of about 1.5 × c from the airfoil/hydrofoil. It could be objected that the experiments adopted in the present work to validate the pure CFD model refer to a turbine [40] with a moderately lower Reynolds and a significantly greater solidity than that of interest. However, we are confident that since the pure CFD model provides realistic results even under more critical conditions for flow separation phenomena (as characterized by a lower Reynolds and higher solidity) than the conditions to which the hybrid models are applied in the present study, it can be considered reliable.

BEM Theory Basics
The hybrid model (regardless of whether it is ANSYS-Fluent or SHYFEM) is based on the BEM theory. According to this theory, the "blade element" is defined as the blade section placed at a generic radius r of infinitesimal thickness δr; the forces acting on the B blade elements (where B is the total blades number) are responsible for the axial momentum loss to which the flow passing through the ring of infinitesimal thickness is subjected. As shown in Figure 3, the resulting force can be split in lift and drag, which are directed normal and parallel to the relative velocity W (vector sum of the flow absolute velocity and the blade velocity) approaching the blade element. Important angles that have to be defined are the attack angle α, the pitch angle β, and the inflow angle Φ, which is a sum of α and β. It should be noted that Φ is the angle between the incoming flow and the rotor plane. The values of lift δL and drag δD acting on the blade element are calculated as follows: where CL and CD are the lift and drag coefficients of the particular hydrofoil that is considered. CL and CD depend on the attack angle and Reynolds number, and they need to be measured experimentally in a wind tunnel. Then, the expressions of torque and power are By applying the principle of action and reaction, equal and opposite forces are impressed on the flow: in fact, in hybrid models, the physical blades are not simulated, and therefore, they are not included in the computational domain. The domain region that is swept by the blades is instead

BEM Theory Basics
The hybrid model (regardless of whether it is ANSYS-Fluent or SHYFEM) is based on the BEM theory. According to this theory, the "blade element" is defined as the blade section placed at a generic radius r of infinitesimal thickness δr; the forces acting on the B blade elements (where B is the total blades number) are responsible for the axial momentum loss to which the flow passing through the ring of infinitesimal thickness is subjected. As shown in Figure 3, the resulting force can be split in lift and drag, which are directed normal and parallel to the relative velocity W (vector sum of the flow absolute velocity and the blade velocity) approaching the blade element. Important angles that have to be defined are the attack angle α, the pitch angle β, and the inflow angle Φ, which is a sum of α and β. It should be noted that Φ is the angle between the incoming flow and the rotor plane.

BEM Theory Basics
The hybrid model (regardless of whether it is ANSYS-Fluent or SHYFEM) is based on the BEM theory. According to this theory, the "blade element" is defined as the blade section placed at a generic radius r of infinitesimal thickness δr; the forces acting on the B blade elements (where B is the total blades number) are responsible for the axial momentum loss to which the flow passing through the ring of infinitesimal thickness is subjected. As shown in Figure 3, the resulting force can be split in lift and drag, which are directed normal and parallel to the relative velocity W (vector sum of the flow absolute velocity and the blade velocity) approaching the blade element. Important angles that have to be defined are the attack angle α, the pitch angle β, and the inflow angle Φ, which is a sum of α and β. It should be noted that Φ is the angle between the incoming flow and the rotor plane. The values of lift δL and drag δD acting on the blade element are calculated as follows: where CL and CD are the lift and drag coefficients of the particular hydrofoil that is considered. CL and CD depend on the attack angle and Reynolds number, and they need to be measured experimentally in a wind tunnel. Then, the expressions of torque and power are By applying the principle of action and reaction, equal and opposite forces are impressed on the flow: in fact, in hybrid models, the physical blades are not simulated, and therefore, they are not included in the computational domain. The domain region that is swept by the blades is instead The values of lift δL and drag δD acting on the blade element are calculated as follows: where C L and C D are the lift and drag coefficients of the particular hydrofoil that is considered. C L and C D depend on the attack angle and Reynolds number, and they need to be measured experimentally in a wind tunnel. Then, the expressions of torque δQ and power δP are δQ = (δLcos(Φ) + δDsen(Φ))r δP = δQ × Ω .
By applying the principle of action and reaction, equal and opposite forces are impressed on the flow: in fact, in hybrid models, the physical blades are not simulated, and therefore, they are not included in the computational domain. The domain region that is swept by the blades is instead represented by fluid cells where momentum source terms are applied to slow down the flow. The source terms are in the form of forces per unit of volume: where F i, j is the force applied on the j-th cell of the grid in the i-th direction, ∆θ, ∆h, and ∆r are the azimuthal, height, and radial extension respectively of the considered blade element.

Grid Effect Analysis With the ANSYS Hybrid Model
The BEM model was implemented by Rocchio et al. [35] in the ANSYS-Fluent general purpose CFD code. In the present work, for all the sub-models, we kept the same coefficients that were tuned in the original work [35], despite the tidal turbine subject of this study having a higher solidity with respect to the original wind turbine. The reason for this choice was to verify the robustness of the model to differences in turbine geometry. In this section, we characterize the BEM-ANSYS hybrid model behavior in terms of predicting the performance of a tidal CFT and the wake generated by it, and we compare it with pure CFD simulations, varying the grid refinement. The aim was to check if a grid refinement compatible with marine circulation codes, and then relatively low resolution, is still enough to describe the details of the blade torque evolution during the turbine rotation and the wake main characteristics.
The ANSYS hybrid model is implemented in the software by means of a User-Defined Function (UDF). The hydrokinetic turbine here considered is the Kobold turbine with three straight blades [42] where , is the force applied on the j-th cell of the grid in the i-th direction, , ℎ, and are the azimuthal, height, and radial extension respectively of the considered blade element.

Grid Effect Analysis With the ANSYS Hybrid Model
The BEM model was implemented by Rocchio et al. [35] in the ANSYS-Fluent general purpose CFD code. In the present work, for all the sub-models, we kept the same coefficients that were tuned in the original work [35], despite the tidal turbine subject of this study having a higher solidity with respect to the original wind turbine. The reason for this choice was to verify the robustness of the model to differences in turbine geometry. In this section, we characterize the BEM-ANSYS hybrid model behavior in terms of predicting the performance of a tidal CFT and the wake generated by it, and we compare it with pure CFD simulations, varying the grid refinement. The aim was to check if a grid refinement compatible with marine circulation codes, and then relatively low resolution, is still enough to describe the details of the blade torque evolution during the turbine rotation and the wake main characteristics.
The ANSYS hybrid model is implemented in the software by means of a User-Defined Function (UDF). The hydrokinetic turbine here considered is the Kobold turbine with three straight blades [42] ( Figure 4a): the height is 6 m, the radius is 3 m, and the chord is 0.4 m. The solidity σ is 6.37%, and it is defined as  The hybrid grid shows the ring (in red in Figure 4b) of fluid cells where source terms are applied. It can be seen that the kind of grid is "structured multi-blocks" that are only made of quadrilateral cells. The domain dimensions are 64R and 50R, respectively in the x and y direction. The time step is set by dividing the revolution period for the number of ring's cells or for 720, respectively for the hybrid model and pure CFD. In both cases, 60 turbine revolutions have been simulated. The flow inlet velocity is fixed at 1.75 m/s; other setups are the same as already listed in Table 1.
The UDF [35] has several sub-models that are used to simulate phenomena such as dynamic stall, tip losses, and virtual cambering, yet in this work, only the dynamic stall sub-model [43] has been considered for simplicity. A sensitivity analysis to the grid refinement was carried out by increasing progressively the number of cells on the ring from 20 to 148. Figure 5a depicts the definitions commonly adopted to describe the blade path of a CFT during a rotation: the upwind path occurs for θ going from 0° to 180°; the downwind path occurs for θ going from 180° to 360°. The hybrid grid shows the ring (in red in Figure 4b) of fluid cells where source terms are applied. It can be seen that the kind of grid is "structured multi-blocks" that are only made of quadrilateral cells. The domain dimensions are 64R and 50R, respectively in the x and y direction. The time step is set by dividing the revolution period for the number of ring's cells or for 720, respectively for the hybrid model and pure CFD. In both cases, 60 turbine revolutions have been simulated. The flow inlet velocity is fixed at 1.75 m/s; other setups are the same as already listed in Table 1.
The UDF [35] has several sub-models that are used to simulate phenomena such as dynamic stall, tip losses, and virtual cambering, yet in this work, only the dynamic stall sub-model [43] has been considered for simplicity. A sensitivity analysis to the grid refinement was carried out by increasing progressively the number of cells on the ring from 20 to 148. Figure 5a depicts the definitions commonly adopted to describe the blade path of a CFT during a rotation: the upwind path occurs for θ going from 0 • to 180 • ; the downwind path occurs for θ going from 180 • to 360 • . Figure 5b shows the evolution of the instantaneous C P of the single blade during one revolution; the hybrid model is sufficiently able to capture the main detail of the C P evolution at the different grid resolutions except for the coarsest case (i.e., 20 cells on the ring). The anomalous and unphysical drop in the upwind curve (θ ≈ 95 • -110 • ) is due to the virtual camber sub-model. In fact, this sub-model is greatly sensitive to the turbine solidity: for solidity lower than 4-5% it is stable, but it does not adequately work for the solidity of the Kobold turbine, therefore leading to the choice to switch off this sub-model.  Figure 5b shows the evolution of the instantaneous CP of the single blade during one revolution; the hybrid model is sufficiently able to capture the main detail of the CP evolution at the different grid resolutions except for the coarsest case (i.e., 20 cells on the ring). The anomalous and unphysical drop in the upwind curve (θ ≈ 95°-110°) is due to the virtual camber sub-model. In fact, this sub-model is greatly sensitive to the turbine solidity: for solidity lower than 4-5% it is stable, but it does not adequately work for the solidity of the Kobold turbine, therefore leading to the choice to switch off this sub-model. It is worth noting that since our present hybrid model implementation is just 2D, at present, is not possible to reproduce the fast wake recovery that can be experimentally observed, and that is mainly due to the vertical advection induced by the evolution of the tip vortices [44]. Therefore, also for the wake characteristics, we preferred to compare/validate the hybrid models against results from pure 2D CFD, since it undergoes the same limitation. Wake velocity profiles of the velocity x and y components are matched for all the resolutions ( Figure 6), with the exception of the coarsest ones (20 and 28 cells on the ring). Based on these preliminary outcomes, the ring was divided into 44 cells in order to have a not too-fine local resolution, and it would be compatible with marine circulation codes. In dimensionless terms, this refinement corresponds to a cell size of 1/7R.   It is worth noting that since our present hybrid model implementation is just 2D, at present, is not possible to reproduce the fast wake recovery that can be experimentally observed, and that is mainly due to the vertical advection induced by the evolution of the tip vortices [44]. Therefore, also for the wake characteristics, we preferred to compare/validate the hybrid models against results from pure 2D CFD, since it undergoes the same limitation. Wake velocity profiles of the velocity x and y components are matched for all the resolutions ( Figure 6), with the exception of the coarsest ones (20 and 28 cells on the ring). Based on these preliminary outcomes, the ring was divided into 44 cells in order to have a not too-fine local resolution, and it would be compatible with marine circulation codes. In dimensionless terms, this refinement corresponds to a cell size of 1/7R.  Figure 5b shows the evolution of the instantaneous CP of the single blade during one revolution; the hybrid model is sufficiently able to capture the main detail of the CP evolution at the different grid resolutions except for the coarsest case (i.e., 20 cells on the ring). The anomalous and unphysical drop in the upwind curve (θ ≈ 95°-110°) is due to the virtual camber sub-model. In fact, this sub-model is greatly sensitive to the turbine solidity: for solidity lower than 4-5% it is stable, but it does not adequately work for the solidity of the Kobold turbine, therefore leading to the choice to switch off this sub-model. It is worth noting that since our present hybrid model implementation is just 2D, at present, is not possible to reproduce the fast wake recovery that can be experimentally observed, and that is mainly due to the vertical advection induced by the evolution of the tip vortices [44]. Therefore, also for the wake characteristics, we preferred to compare/validate the hybrid models against results from pure 2D CFD, since it undergoes the same limitation. Wake velocity profiles of the velocity x and y components are matched for all the resolutions ( Figure 6), with the exception of the coarsest ones (20 and 28 cells on the ring). Based on these preliminary outcomes, the ring was divided into 44 cells in order to have a not too-fine local resolution, and it would be compatible with marine circulation codes. In dimensionless terms, this refinement corresponds to a cell size of 1/7R.     The qualitative trend is similar for the hybrid and the CFD cases at different TSR values while quantitative discrepancies are noticeable; these differences are due to the fact that the ANSYS hybrid model was validated versus experimental data concerning low solidity wind turbines, as reported in [35]. Yet, it should be noted that the qualitative behavior is captured anyway. For a better quantitative matching of the turbine behavior, it will be necessary to tune some of the sub-model's coefficients versus experimental data available for higher solidity CFTs. However, in this work, we are not interested in improving the hybrid model developed in [35] for a water system: we will assume the model to be acceptable in order to capture the qualitative behavior of a hydrokinetic device, and the ANSYS hybrid model results will be adopted to make a comparison with the SHYFEM hybrid model. What is expected is to record the same behavior for the turbine in both models, because the analytical model underlying the SHYFEM hybrid model is the same as that of the ANSYS hybrid model.
It must be specified that for the hybrid simulations represented in Figure 7, the virtual camber model was disabled: in purely analytical and in hybrid BEM-CFD tools, this model is useful to reproduce the curvature of the flow on the hydrofoil during its rotational motion [45]. In a pure CFD simulation, this effect is naturally reproduced, while in a hybrid model, without a physical blade, this effect is taken into account by giving a virtual camber to the hydrofoil. By doing this, the symmetrical hydrofoil in CFD corresponds to the cambered hydrofoil in hybrid model, as shown in Figure 8a, where a is the axial induction factor. A curved hydrofoil in CFD corresponds to a symmetrical hydrofoil in hybrid simulations, as shown in Figure 8b. In order to keep the virtual camber sub-model in our hybrid simulation switched off, it was necessary to adopt for the pure CFD a hydrofoil with an opposite curvature in order to compensate, and then to cancel, the flow curvature.

SHYFEM Grid
To reproduce the same resolution in the turbine zone, the SHYFEM finite-elements grid has a ring with 56 triangular elements ( Figure 9) with an average size corresponding to a little more than 1/9R. The dimensions of the computing domain are the same as those used for the ANSYS grid, and in this case, the resolution is not uniform all over the grid: an unstructured mesh easily allows high refinement in the area of interest and progressively lowers the resolution, up to 4/3R at the domain boundary. At 2R downstream from the turbine axis, the element's averaged size is about 1/3R, while at 8R, it is about 2/3R.

Friction's Formulation in the SHYFEM Hybrid Model
Differently from ANSYS Fluent, SHYFEM simulates the effect of the turbine's presence not as a momentum source but as a variation on bottom friction. In the present formulation, the model solves the shallow water equations, reproducing the barotropic transport and reflecting the horizontal flow patterns around the turbines: where H is the height of the water column given by the sum of ζ, the level of the free surface of the water, which is variable due to, e.g., waves and tides, and h, the undisturbed water column depth; F is a coefficient that describes the time scale of friction; g is the gravitational acceleration; and X and

SHYFEM Grid
To reproduce the same resolution in the turbine zone, the SHYFEM finite-elements grid has a ring with 56 triangular elements ( Figure 9) with an average size corresponding to a little more than 1/9R.

SHYFEM Grid
To reproduce the same resolution in the turbine zone, the SHYFEM finite-elements grid has a ring with 56 triangular elements ( Figure 9) with an average size corresponding to a little more than 1/9R. The dimensions of the computing domain are the same as those used for the ANSYS grid, and in this case, the resolution is not uniform all over the grid: an unstructured mesh easily allows high refinement in the area of interest and progressively lowers the resolution, up to 4/3R at the domain boundary. At 2R downstream from the turbine axis, the element's averaged size is about 1/3R, while at 8R, it is about 2/3R.

Friction's Formulation in the SHYFEM Hybrid Model
Differently from ANSYS Fluent, SHYFEM simulates the effect of the turbine's presence not as a momentum source but as a variation on bottom friction. In the present formulation, the model solves the shallow water equations, reproducing the barotropic transport and reflecting the horizontal flow patterns around the turbines: where H is the height of the water column given by the sum of ζ, the level of the free surface of the water, which is variable due to, e.g., waves and tides, and h, the undisturbed water column depth; F is a coefficient that describes the time scale of friction; g is the gravitational acceleration; and X and The dimensions of the computing domain are the same as those used for the ANSYS grid, and in this case, the resolution is not uniform all over the grid: an unstructured mesh easily allows high refinement in the area of interest and progressively lowers the resolution, up to 4/3R at the domain boundary. At 2R downstream from the turbine axis, the element's averaged size is about 1/3R, while at 8R, it is about 2/3R.

Friction's Formulation in the SHYFEM Hybrid Model
Differently from ANSYS Fluent, SHYFEM simulates the effect of the turbine's presence not as a momentum source but as a variation on bottom friction. In the present formulation, the model solves the shallow water equations, reproducing the barotropic transport and reflecting the horizontal flow patterns around the turbines: where H is the height of the water column given by the sum of ζ, the level of the free surface of the water, which is variable due to, e.g., waves and tides, and h, the undisturbed water column depth; F is a coefficient that describes the time scale of friction; g is the gravitational acceleration; and X and Y are explicit terms forcing the hydrodynamics, such as wind and atmospheric pressure. Finally, U and V are transports.
The shallow water approximation assumes that the horizontal scale of the process is larger than the vertical scale, and the water column is reproduced by a unique layer. SHYFEM, in its 3D formulation, assumes the Shallow Water approximation in order to solve the flow field in each layer: then, vertical advection and viscosity terms induce the vertical dynamics, allowing connection among layers on the water column; this means that an exchange of mass, temperature, and other variables occur between different layers. In the present implementation, as a first substantial simplification, the focus is on the reproduction of the horizontal current field around the turbine; therefore, a 2D approach is preferred, and the effects of the turbine's presence in water in terms of slowdown on barotropic flows are induced by variation in bottom friction. Since this is the first application of a hybrid model with a BEM approach in a marine circulation model, it was decided to adopt the friction method to allow a more stable calculation. This kind of approach is suitable for 2D application, while it would not be suitable for a 3D simulation. In coastal and transitional environment studies, the different friction formulations are used to reproduce different bottom environments (e.g., rocky, muddy, sandy) and concur in allowing correct sediment transport dynamic, deposition-resuspension, or wave pattern reproduction [46,47]. The friction is introduced in the Shallow Water equation thanks to the F term defined as follows: where λ represent the drag coefficient (a-dimensional) that can be set constant. The formulation of λ in the grid area where the turbine is located, according to other works found in the literature [1,48], is given by: where λ is proportional to the force module that the blade would have experimented. Lift and Drag (L and D) forces are given for a single blade and are expressed per unit area of the grid element to which they are applied. Forces are applied to a ring's cell as if a single blade is present in that cell for the entire revolution period. Since this is similar to having the blade omnipresent in the whole ring at any time, it is necessary to divide the forces for the time frame in which the blade is effectively present in the cell, dividing the forces for the number of the ring cells icount: the effect is to have a blade that rests on every element just for an icount-th of the revolution time. Furthermore, since there are B blades, it is necessary to multiply forces for the number of blade n B . Developing the equation, we obtain: Applying friction only on the ring's elements, we can see fluctuations in the downwind performance curve, as shown in Figure 10. This simulation has been made at the optimal Tip Speed Ratio of 2.7.
Applying friction only on the ring's elements, we can see fluctuations in the downwind performance curve, as shown in Figure 10. This simulation has been made at the optimal Tip Speed Ratio of 2.7. The dynamic stall model in downwind leads to unrealistic downwind fluctuations; therefore, it was decided to disable it. This approximation is acceptable, because the downwind flux is slowed down by the power generation occurred in upwind: it is unlikely of a having high angle of attack that can cause downwind stall.
On the other hand, applying the friction only on the ring's elements (with the dynamic stall model disabled in downwind), the performances are not similar to those obtained by pure CFD simulations: the upwind production is reduced, and the downwind production is excessive ( Figure  10). Therefore, a different approach was needed to improve the hybrid model results: the turbinedependent friction has been applied not only on the ring's elements but also on the inside elements. For the inside elements, the θ value is given in a different way: for the upwind elements, we divide the turbine in 18 virtual stripes and assign the θ value for all the inside elements as the one of the ring's element on the same stripe ( Figure 11a); for the inside downwind elements, the θ value is assigned as the geometrical angle of the element (Figure 11b).  Applying the friction on the entire disc, the flow that passes through the central stripe of the turbine is excessively slowed down, since there are many elements with turbine-dependent friction. On the contrary, the flow that passes through lateral stripes of the turbine is not slowed down enough. Making the turbine virtually square prevents these inconveniences, multiplying the number of turbine elements by a 1.2732 factor that represents the ratio between the surface of a square with The dynamic stall model in downwind leads to unrealistic downwind fluctuations; therefore, it was decided to disable it. This approximation is acceptable, because the downwind flux is slowed down by the power generation occurred in upwind: it is unlikely of a having high angle of attack that can cause downwind stall.
On the other hand, applying the friction only on the ring's elements (with the dynamic stall model disabled in downwind), the performances are not similar to those obtained by pure CFD simulations: the upwind production is reduced, and the downwind production is excessive ( Figure 10). Therefore, a different approach was needed to improve the hybrid model results: the turbine-dependent friction has been applied not only on the ring's elements but also on the inside elements. For the inside elements, the θ value is given in a different way: for the upwind elements, we divide the turbine in 18 virtual stripes and assign the θ value for all the inside elements as the one of the ring's element on the same stripe ( Figure 11a); for the inside downwind elements, the θ value is assigned as the geometrical angle of the element (Figure 11b). Applying friction only on the ring's elements, we can see fluctuations in the downwind performance curve, as shown in Figure 10. This simulation has been made at the optimal Tip Speed Ratio of 2.7. The dynamic stall model in downwind leads to unrealistic downwind fluctuations; therefore, it was decided to disable it. This approximation is acceptable, because the downwind flux is slowed down by the power generation occurred in upwind: it is unlikely of a having high angle of attack that can cause downwind stall.
On the other hand, applying the friction only on the ring's elements (with the dynamic stall model disabled in downwind), the performances are not similar to those obtained by pure CFD simulations: the upwind production is reduced, and the downwind production is excessive ( Figure  10). Therefore, a different approach was needed to improve the hybrid model results: the turbinedependent friction has been applied not only on the ring's elements but also on the inside elements. For the inside elements, the θ value is given in a different way: for the upwind elements, we divide the turbine in 18 virtual stripes and assign the θ value for all the inside elements as the one of the ring's element on the same stripe ( Figure 11a); for the inside downwind elements, the θ value is assigned as the geometrical angle of the element (Figure 11b).  Applying the friction on the entire disc, the flow that passes through the central stripe of the turbine is excessively slowed down, since there are many elements with turbine-dependent friction. On the contrary, the flow that passes through lateral stripes of the turbine is not slowed down enough. Making the turbine virtually square prevents these inconveniences, multiplying the number of turbine elements by a 1.2732 factor that represents the ratio between the surface of a square with Applying the friction on the entire disc, the flow that passes through the central stripe of the turbine is excessively slowed down, since there are many elements with turbine-dependent friction. On the contrary, the flow that passes through lateral stripes of the turbine is not slowed down enough. Making the turbine virtually square prevents these inconveniences, multiplying the number of turbine elements by a 1.2732 factor that represents the ratio between the surface of a square with edge 2R and a disc with radius R, and then adopting a 2R/chord factor to correct the friction in each element (see Figure 12). The chord is the geometrical chord of the circle. edge 2R and a disc with radius R, and then adopting a 2R/chord factor to correct the friction in each element (see Figure 12). The chord is the geometrical chord of the circle. This method allows having an appropriate imposition of turbine-dependent friction to a wider surface without having problems, such as premature slowdown of the flow, to reproduce the blade's power production along a revolution.

Results of the BEM-SHYFEM Hybrid Model
The chapter shows the predictive capability of the SHYFEM hybrid model for a single turbine in comparison to a pure CFD and ANSYS hybrid model and for two clusters of three turbines in comparison to pure CFD. Some considerations about grid resolution effects are presented.

Behavior of the Single Turbine at Different TSR
Several simulations have been performed at different TSR, which were obtained by varying the rotational speed of the turbine and maintaining the free stream velocity fixed to 1.75 m/s. Other important setups are listed in Table 2. The results obtained with SHYFEM are very close to those obtained with ANSYS's UDF, as shown in Figure 13. The behavior of the turbine is reported for different TSRs. It is noteworthy that not only the turbine overall CP sufficiently matches the respective CFD result for all the TSR, but also the behavior of blade instantaneous CP satisfies the typical azimuthal evolution in CFTs, showing, for instance, dynamic stall in upwind at very low TSR and power output decreasing in downwind as the TSR increases. What emerges is a quantitative difference between the two hybrid models and the pure CFD: as explained before, this is due to the fact that the model developed in [35] is tuned for an airfoil with lower solidity. We report also the ANSYS hybrid model in graphics just to make a comparison between the two hybrid models. Since both are based on the same analytical model, the same behavior is expected. As shown in Figure 13, expectations are satisfied for all TSR. By comparing the values of the source terms in ANSYS and SHYFEM simulations (Table 3), it is evident that the total source terms are similar: differences are limited to 10% (the only exception is TSR 3.2). This method allows having an appropriate imposition of turbine-dependent friction to a wider surface without having problems, such as premature slowdown of the flow, to reproduce the blade's power production along a revolution.

Results of the BEM-SHYFEM Hybrid Model
The chapter shows the predictive capability of the SHYFEM hybrid model for a single turbine in comparison to a pure CFD and ANSYS hybrid model and for two clusters of three turbines in comparison to pure CFD. Some considerations about grid resolution effects are presented.

Behavior of the Single Turbine at Different TSR
Several simulations have been performed at different TSR, which were obtained by varying the rotational speed of the turbine and maintaining the free stream velocity fixed to 1.75 m/s. Other important setups are listed in Table 2. The results obtained with SHYFEM are very close to those obtained with ANSYS's UDF, as shown in Figure 13. The behavior of the turbine is reported for different TSRs. It is noteworthy that not only the turbine overall C P sufficiently matches the respective CFD result for all the TSR, but also the behavior of blade instantaneous C P satisfies the typical azimuthal evolution in CFTs, showing, for instance, dynamic stall in upwind at very low TSR and power output decreasing in downwind as the TSR increases. What emerges is a quantitative difference between the two hybrid models and the pure CFD: as explained before, this is due to the fact that the model developed in [35] is tuned for an airfoil with lower solidity. We report also the ANSYS hybrid model in graphics just to make a comparison between the two hybrid models. Since both are based on the same analytical model, the same behavior is expected. As shown in Figure 13, expectations are satisfied for all TSR. By comparing the values of the source terms in ANSYS and SHYFEM simulations (Table 3), it is evident that the total source terms are similar: differences are limited to 10% (the only exception is TSR 3.2).  This guarantees that the energy balance is respected, although the friction's ubication has been dislocated with respect to power production, as described in Section 2.5. However, it should be noted that the relative amounts of the x and y sources are significantly different in ANSYS and in SHYFEM. This is because SHYFEM uses a unique formulation for the friction term F, while in the ANSYS's  This guarantees that the energy balance is respected, although the friction's ubication has been dislocated with respect to power production, as described in Section 2.5. However, it should be noted that the relative amounts of the x and y sources are significantly different in ANSYS and in SHYFEM. This is because SHYFEM uses a unique formulation for the friction term F, while in the ANSYS's UDF, the x and y force components to be adopted in the momentum equations are calculated by correctly projecting the lift and drag along the x and y directions.
By a qualitative analysis of the flow field velocity magnitude (Figure 14), it can be seen that SHYFEM results are similar to those obtained with ANSYS's UDF. However, in the SHYFEM hybrid model, the grid resolution around and downstream from the turbine notably affects the wake's length and the magnitude of the flow lateral accelerations. More details will be provided in Section 3.2.
J. Mar. Sci. Eng. 2020, 8, 1010 14 of 21 UDF, the x and y force components to be adopted in the momentum equations are calculated by correctly projecting the lift and drag along the x and y directions. By a qualitative analysis of the flow field velocity magnitude (Figure 14), it can be seen that SHYFEM results are similar to those obtained with ANSYS's UDF. However, in the SHYFEM hybrid model, the grid resolution around and downstream from the turbine notably affects the wake's length and the magnitude of the flow lateral accelerations. More details will be provided in Section 3.2.

Application to a Small Turbine's Cluster
Then, the SHYFEM model was applied to a set of three turbines in two different spatial configurations: in the first case, turbines are placed aligned with a distance between axes of 3R, while in the second case, the devices are placed at the vertices of an equilateral triangle, with distance between axes of 9.2R, as sketched in Figure 15. In this analysis, two different kinds of grid were used: one with uniform resolution all over the domain fixed at 1/6R and the other with variable resolution from about 1/6R in the turbine zone to 4/3R at the boundaries. The extension of these grids was increased to 100R along x and 80R along y (300 m and 240 m long x and y respectively) to avoid boundary blockage effects.  Figure 16 shows the flow field obtained with SHYFEM: it can be noticed, as mentioned before, how grid resolution affects the flow field. Higher resolution grids allow better visualizing velocity gradients concentrated in limited areas; furthermore, the wake length increases. This happens for both spatial configurations.
For the aligned triad, acceleration corridors appear between turbines, since the flow cannot expand freely; therefore, there is a decrease in the y component of speed and consequently also an

Application to a Small Turbine's Cluster
Then, the SHYFEM model was applied to a set of three turbines in two different spatial configurations: in the first case, turbines are placed aligned with a distance between axes of 3R, while in the second case, the devices are placed at the vertices of an equilateral triangle, with distance between axes of 9.2R, as sketched in Figure 15. In this analysis, two different kinds of grid were used: one with uniform resolution all over the domain fixed at 1/6R and the other with variable resolution from about 1/6R in the turbine zone to 4/3R at the boundaries. The extension of these grids was increased to 100R along x and 80R along y (300 m and 240 m long x and y respectively) to avoid boundary blockage effects. UDF, the x and y force components to be adopted in the momentum equations are calculated by correctly projecting the lift and drag along the x and y directions. By a qualitative analysis of the flow field velocity magnitude (Figure 14), it can be seen that SHYFEM results are similar to those obtained with ANSYS's UDF. However, in the SHYFEM hybrid model, the grid resolution around and downstream from the turbine notably affects the wake's length and the magnitude of the flow lateral accelerations. More details will be provided in Section 3.2.

Application to a Small Turbine's Cluster
Then, the SHYFEM model was applied to a set of three turbines in two different spatial configurations: in the first case, turbines are placed aligned with a distance between axes of 3R, while in the second case, the devices are placed at the vertices of an equilateral triangle, with distance between axes of 9.2R, as sketched in Figure 15. In this analysis, two different kinds of grid were used: one with uniform resolution all over the domain fixed at 1/6R and the other with variable resolution from about 1/6R in the turbine zone to 4/3R at the boundaries. The extension of these grids was increased to 100R along x and 80R along y (300 m and 240 m long x and y respectively) to avoid boundary blockage effects.  Figure 16 shows the flow field obtained with SHYFEM: it can be noticed, as mentioned before, how grid resolution affects the flow field. Higher resolution grids allow better visualizing velocity gradients concentrated in limited areas; furthermore, the wake length increases. This happens for both spatial configurations.
For the aligned triad, acceleration corridors appear between turbines, since the flow cannot  Figure 16 shows the flow field obtained with SHYFEM: it can be noticed, as mentioned before, how grid resolution affects the flow field. Higher resolution grids allow better visualizing velocity gradients concentrated in limited areas; furthermore, the wake length increases. This happens for both spatial configurations.
respectively at the beginning and at the end of upwind path. In fact, in Figure 17, it is shown that turbine 2 anticipates the production, while turbine 3 postpones production with respect to the isolated turbine. The behavior of the three turbines is in accordance with pure CFD results. Turbines 2 and 3 have clockwise rotational verse: in this case, the θ reference system is not the same as that proposed in Figure 5a but it is opposite with respect to the y axis. The direction of increase of the angle θ is concordant with the direction of rotation.
Instead, the higher production recorded in downwind for all three turbines is due to the fact that wakes are not free to expand; therefore, speeds in wake are higher compared to the isolated turbine case. In fact, if the wake is constrained by the presence of other devices' wakes, the section of the stream tube is reduced, and consequently, higher velocity is recorded. In this way, the mass balance is respected. For the triad in a triangular configuration, the incoming flow speed and then the power output are reduced for turbine 1 (Figure 18) due to the greater fluid dynamic resistance offered to the flow by the triad compared to the single device. Turbines 4 and 5, on the other hand, increase power generation in upwind, since the accelerated water flowing laterally to turbine 1 is projected onto the side turbines, which therefore benefit from higher flow velocities. The triangular triad in pure CFD simulations shows the same behavior. For the aligned triad, acceleration corridors appear between turbines, since the flow cannot expand freely; therefore, there is a decrease in the y component of speed and consequently also an increase in the angle of attack. The higher velocities together with the greater attack angles lead to significantly greater performance, where the acceleration corridors are located with respect to the blade path: considering turbine 1, acceleration corridors are located on both sides; this means at the beginning and at the end of upwind path for this turbine, this leads to a larger upwind C P curve than in the isolated case, as confirmed in Figure 17. For turbines 2 and 3, acceleration corridors are present respectively at the beginning and at the end of upwind path. In fact, in Figure 17, it is shown that turbine 2 anticipates the production, while turbine 3 postpones production with respect to the isolated turbine. The behavior of the three turbines is in accordance with pure CFD results. Turbines 2 and 3 have clockwise rotational verse: in this case, the θ reference system is not the same as that proposed in Figure 5a but it is opposite with respect to the y axis. The direction of increase of the angle θ is concordant with the direction of rotation.
Instead, the higher production recorded in downwind for all three turbines is due to the fact that wakes are not free to expand; therefore, speeds in wake are higher compared to the isolated turbine case. In fact, if the wake is constrained by the presence of other devices' wakes, the section of the stream tube is reduced, and consequently, higher velocity is recorded. In this way, the mass balance is respected.  In Table 4, the percentage variations in performance for each single device making part of the triad, compared to the same device working alone and the performance changes for the entire cluster are presented. The results show a qualitative trend similar to the pure CFD behavior. The major quantitative differences are recorded for the side-by-side triad devices.  For the triad in a triangular configuration, the incoming flow speed and then the power output are reduced for turbine 1 (Figure 18) due to the greater fluid dynamic resistance offered to the flow by the triad compared to the single device. Turbines 4 and 5, on the other hand, increase power generation in upwind, since the accelerated water flowing laterally to turbine 1 is projected onto the side turbines, which therefore benefit from higher flow velocities. The triangular triad in pure CFD simulations shows the same behavior.  In Table 4, the percentage variations in performance for each single device making part of the triad, compared to the same device working alone and the performance changes for the entire cluster are presented. The results show a qualitative trend similar to the pure CFD behavior. The major quantitative differences are recorded for the side-by-side triad devices.  In Table 4, the percentage variations in performance for each single device making part of the triad, compared to the same device working alone and the performance changes for the entire cluster are presented. The results show a qualitative trend similar to the pure CFD behavior. The major quantitative differences are recorded for the side-by-side triad devices. Wake velocity profiles are well reproduced by the SHYFEM hybrid model compared to pure CFD simulations, as shown in Figure 19. The influence of the resolution is more evident at a distance of 6R behind the turbine for both spatial configurations of turbines: this is because at 3R, the two grids still have comparable resolution, while at 6R, the difference is greater. At 3R in SHYFEM, the flow is slowed down too much, which is in accordance to the higher entity of x source terms (as shown in Table 3). The resolution affects the results at a distance of 6R where the higher resolution grid has a slower wake recovery and exhibits more distant behavior than pure CFD. It is useful to remark that these are 2D simulations: in this case, wakes are much longer compared to 3D wakes because 2D simulations do not allow vertical advection, which supports the re-energization of the wake [44]. Since 3D wakes will be shorter [13], the area where it will be necessary to have high resolution behind the turbine will be limited. Resolution has more importance in 2D than in 3D. Therefore, for all the above considerations, it can be concluded that it is not necessary to have high resolution all over the domain in order to simulate the turbine effect on the surrounding flow field. This allows enormous time savings in computing time. Wake velocity profiles are well reproduced by the SHYFEM hybrid model compared to pure CFD simulations, as shown in Figure 19. The influence of the resolution is more evident at a distance of 6R behind the turbine for both spatial configurations of turbines: this is because at 3R, the two grids still have comparable resolution, while at 6R, the difference is greater. At 3R in SHYFEM, the flow is slowed down too much, which is in accordance to the higher entity of x source terms (as shown in Table 3). The resolution affects the results at a distance of 6R where the higher resolution grid has a slower wake recovery and exhibits more distant behavior than pure CFD. It is useful to remark that these are 2D simulations: in this case, wakes are much longer compared to 3D wakes because 2D simulations do not allow vertical advection, which supports the re-energization of the wake [44]. Since 3D wakes will be shorter [13], the area where it will be necessary to have high resolution behind the turbine will be limited. Resolution has more importance in 2D than in 3D. Therefore, for all the above considerations, it can be concluded that it is not necessary to have high resolution all over the domain in order to simulate the turbine effect on the surrounding flow field. This allows enormous time savings in computing time. It should be remarked that our current 2D SHYFEM hybrid model is similar to having a single layer to analyze, in which the depth is equal to the turbine height; consequently, none of the following 3D phenomena can be considered: accelerations of the flow above and below the turbine due to vertical blockage, vortices at the blade tips, mixing with re-energization of the wake by vertical It should be remarked that our current 2D SHYFEM hybrid model is similar to having a single layer to analyze, in which the depth is equal to the turbine height; consequently, none of the following 3D phenomena can be considered: accelerations of the flow above and below the turbine due to vertical blockage, vortices at the blade tips, mixing with re-energization of the wake by vertical convection. Although this limitation allows matching the results from pure 2D CFD simulations (since in those cases, it is similar to carrying out the calculations on the mid plane of a turbine with blades so long as the effects of fluid dynamic tip losses are avoided), it prevents the application of the hybrid model to identify the optimal layouts for turbine clusters. As future work, we plan to extend the embedding of the BEM model to 3D. In realistic 3D simulations, in which the water column extends above and below the turbine, the friction method should be replaced with a methodology capable of reproducing the accelerations above and below the turbine (vertical blockage). To this end, it will be necessary to evaluate whether to switch to the "source terms" method [28]. Some studies in fact show that for a horizontal axis turbine, there are no major differences in using a 3D with "source terms" or a 2D with "Bed friction" if preliminary 3D simulations are made in order to calculate the integrated average speed [49].

Conclusions
The Blade Element turbine analytical model for cross-flow turbines developed for SHYFEM allows obtaining results close to the pure CFD, and to the hybrid BEM-CFD model previously implemented in the ANSYS-Fluent software, yet with much lower computational effort and therefore with considerable time savings. To have an order of magnitude SHYFEM, the simulation time for a single device is about 10 min, having a set simulation end time of 1500 s, while for a small cluster of turbines, the simulation time is about 8 h with the same simulation end time. A pure CFD simulation with three rotors and the same domain extension requires about 7 days, so calculation of the time savings is relevant. It is possible to reach these results thanks to grid local refinements that are of the same order of magnitude of the finest ones nowadays adopted by other researchers using the Actuator Disc approach. Yet, in our model, turbine geometry and rotational speed are considered.
In particular, it was verified that just about 50 triangular elements on the annular ring swept by the blades are enough to reproduce even the evolution of the instantaneous C P of the blade during a turbine revolution. This correspond to cell sizes of about 1/17 of the turbine diameter. The SHYFEM hybrid model matches the results of the pure CFD and ANSYS hybrid model also in a wide range of TSR values.
Furthermore, considering the total amount of the source terms, they differ slightly (about 10%) from the results obtained with the ANSYS hybrid model at various TSR: this ensures that although the power generation area has been decoupled from the application of the friction in the turbine-dependent friction approach, the energy balance is in any case respected.
The sensitivity analysis to the grid refinement highlights the necessity of having an adequate resolution in the turbine area and in the wake region, since C P and the turbine effects on the surrounding and on the downstream flows depends on grid spacing.
The hybrid BEM-SHYFEM model was applied to a cluster of three turbines in two different arrangements: the results are qualitatively close to those obtained with pure CFD simulations, in terms not only of C P increasing/decreasing with respect to the isolated turbine, but also regarding the local mutual interactions between the devices. This is not trivial, given the highly asymmetrical behavior of CFTs with respect to the plain passing through the turbine axis and parallel to the current.
After this first preliminary study, the hybrid BEM-SHYFEM model for cross-flow turbines will be adapted for 3D simulations in order to also describe crucial 3D mechanisms such as blade tip losses, vertical blockage, vertical advection, and wake recovery. The aim is to implement a tool that is useful to plan optimal layouts of tidal turbines, with the dual purpose of maximizing the mechanical energy production and minimizing the environmental effects in terms of hydrodynamic impact and turbulence generation. The difference from current studies is that it will be possible to embed turbines in real-life applications and the flow field thanks to the ability of SHYFEM to describe bathymetry, currents, tides, and surface waves as they occur in actual operating conditions. Turbines will be surrounded by a realistic flow field as if they were placed in a real area of interest; by means of the hybrid model, which is able to capture not only the average power generation of the turbine but the exact behavior of a single blade during its rotation, it will be possible to make accurate studies about energy harvested from these devices and their impact on the flow field.