A Double Multiple Stream Tube (DMST) routine for site assessment to select efficient turbine aspect ratios and solidities in real marine environments

A MATLAB routine, based on a Double Multiple Stream Tube model, developed to quickly predict the performance of cross-flow hydrokinetic turbine, here is presented. The routine evaluate flow data obtained with the open-source marine circulation code SHYFEM. The tool can establish the best locations to place tidal devices taking into account bathymetric constraints and the hydrokinetic potential. Hence, it can be used to decide the best set of geometrical parameters. The geometrical variables of our analysis are turbine frontal area, aspect ratio and solidity. Several sub-models, validated with 3D and 2D CFD simulations, reproduce phenomena such as dynamic stall, fluid dynamic tips losses and the lateral deviation of streamlines approaching the turbine. As a case study, the tool is applied to an area of the northern Adriatic Sea. After having identified some suitable sites to exploit the energy resource, we have compared behaviours of different turbines. The set of geometrical parameters that gives the best performance in terms of power coefficient can vary considering several locations. Conversely, the power production is always greater for turbine with low aspect ratio (for a fixed solidity and area). Indeed, shorter devices benefit from higher hydrokinetic potentials at the top of the water column.


Introduction
Cross Flow Tidal Turbines (CFTTs) have been widely studied in the last decades. They offer many advantages compared to Horizontal Axis Tidal Turbines (HATTs). For example, it is possible to use floating platforms, and set some components above the sea level, such as the generator and the gearbox. In this way, it is easier to reach them for maintaining operations, and at the same time management costs are reduced. Moreover, CFTTs are characterized by hydrodynamic mechanisms that are able to quickly recover wakes [1]. This fact leads to different farm layouts compared with those made up with HATTs: devices can be placed in a closer manner, so it is possible to reach higher values of obtained power per unit area. This is useful in order to limit sea areas where the exploitation of the tidal energy resource occurs. To this end, the effect of adopting counterrotating devices is also beneficial [2][3][4]. Indeed, this expedient leads to favourable hydrodynamic conditions.
The peculiarity of the tidal resource is to be predictable: tidal cycles repeat regularly. The natural phenomenon is characterized by the inversion of the flow direction: inversion does not necessarily occur with a relative angle of 180 degrees between the velocity vectors. The angle between them can also be lower, but CFTTs have the further advantage to work independently from flow direction. Conversely, HATTs must be equipped with yaw systems able to move the device, and follow flow changes. Indeed, it is necessary to set the device so that the flow is perpendicular to the rotor plane to produce energy. This contributes to increase the system complexity and decreases reliability.
CFTTs show lower starting-torque and efficiency compared to HATTs: this flaw can be mitigated by using a blade-pitching mechanism [5][6][7].
In this paper we are going to use a Double Multiple Stream Tube (DMST) model, integrated in a MATLAB routine, to predict turbines behaviours. It is a simplified approach based on Blade Element Momentum (BEM) theory [8]. We couple the DMST model with flow data obtained with the SHYFEM [9] code. It is an open access numerical marine circulation code developed by CNR ISMAR of Venice.
In section two we explain the methodology followed for the analysis, section three is dedicated to show results, and section four summarize the most important evaluations giving some conclusions.

Methodology
The DMST model used in this work has been developed at the University of Pisa [10]. Such a model is based on the Actuator Disc (AD) approach, and uses BEM theory in order to compute forces acting on blades. The model allows to evaluate different situations: indeed, it is possible to vary parameters such as Aspect Ratio (AR), solidity σ and the reference value of the Tip Speed Ratio (TSR) in order to calculate the rotational speed ω. The first two variables are defined as: where L is blades' length, D the turbine diameter, c is the chord and Nb is the number of blades.
The TSR is defined as: where R is the turbine radius and U∞ is the undisturbed flow velocity. We first identify the sites with the greater energy potential and at the same time deep enough to host turbine of various vertical dimensions. Then, we are going to analyse the behaviour of turbine with different geometric characteristics: frontal area, AR, σ.
As will be shown in this paper, this method allows to identify the most efficient turbine geometry for a site with a high energy potential and characterized by bathymetric constraints. Let's consider a generic horizontal plane of a CFTT of height L as shown in figure 1. The incoming flow has an undisturbed velocity equal to U∞ that becomes U1 as the flow reaches the turbine. From a simple geometrical consideration it follows that:

BEM momentum source computation
Re(θ) = c W(θ)/ν (8) where W is the relative velocity to the hydrofoil, α(θ) and Re(θ) are the local values of the angle of attack and Reynolds number respectively. These last two parameters are used in order to calculate the local value of the Lift coefficient CL(α,Re), and the Drag coefficient CD(α,Re). Then, it is possible to evaluate the non-dimensional forces coefficient along tangential and radial direction (CT and CN), and along flow direction and perpendicular to flow direction (CX and CY respectively).

CT = CL sin(α) -CD cos(α)
The extracted power P and the power coefficient CP can be evaluated as: where ρ is water or air density, z is the generic vertical quote along turbine height, and A is the turbine frontal area. The power is integrated on θ between 0 and 2π and on z between 0 and L.

Routine working mode
The basis of the analysis carried out in this work, relies on the DMST model. A DMST model is a simplified approach used to reach a quick performance prediction for a turbine in steady conditions.
The turbine region is represented by dividing it into parallel stream tubes and each of them is divided into two parts: the upstream tube and the downstream tube, as shown in figure 2. The θ value for each half streamtube has a unique value, and it corresponds to the θ value of the centre line of the streamtube. Mass, momentum and energy balance are solved for each half-stream tube using an iterative process.

Fig. 2.
Example of DMST schematic representation [10] Balance equations are the following: where the subscript 1 refers to the upstream tube and the subscript 2 refers to the downstream tube, ṁ is the flow rate. Since A1=A2, the mass balance cannot be respected but the error made is assumed to be small, so it can be neglected.
Introducing the induction factors as a1=U1/U∞ and a2=U2/Ue and combining previous equation we obtain: Knowing U1 and U2 we can calculate the total thrust force using equations (4) to (11) related to the total thrust coefficient CX. Then it is possible to define thrust coefficients from the total thrust force : Using equations (15) to (21) and (23) we obtain: This theoretical result does not match with experimental data for ai <0.6. Therefore, the following empirical correlation is adopted (provided by [11]) : Now it is possible, with an iterative process, to assign the ai value until FX obtained from equations (22) and (23) match. This iterative process is embedded in a MATLAB routine.

Site assessment
As previously mentioned, we evaluated different combination of geometrical parameters: in particular, two values of frontal area (25 m 2 and 50 m 2 ) and, for each of them, four values for the AR, as summarized in table 1. There is a third variable that is σ: we considered two solidities 0,0637 and 0,159, called σ1 and σ2 respectively in the following. Since, solidity influences tips losses [12], and consequently turbine performance, we decided to check different values to assess relative effect.
The geographical area of interest in this study is located in the Northern Adriatic Sea (latitude from 44.5 to 45.5 and longitude from 12 to 13). Available flow data, from a 3D SHYFEM simulation, cover a period from 7 th to 21 th of February 2014, with hourly output (337 hours of simulation, about half lunar cycle).
The choice of turbines location was guided by energy assessments. We first delimited the area of searching by deleting sites not deep enough to host even the shortest turbine. Furthermore, for each turbine we kept a 2 m distance from the free surface of the water column and from the seabed. A proper distance between the blade tips and the free surface is necessary to exploit the beneficial effects of bypassing flow, and to allow the expansion and complete development of the wake [13,14]. Whereas, an adequate distance from the seabed is necessary to avoid an increase in bed shear stress and consequent erosion. This is one of the most impacting effect caused by marine turbine farms [15][16]. Then, we evaluated the averaged available power in time during the entire period. The criterion was to accept only sites that reach at least 50% of the maximum value of available power over the entire domain. A total of 10 sites met our requests, as shown in figure 3. We can see that the most suitable sites are near the three inlets connecting the Venice Lagoon to the Adriatic Sea. The morphology of these sites tends to channel the flow, causing an increase in speed.
To run the routine it is necessary to give as an input a velocity profile along the vertical direction. The purpose is not to evaluate the energy producibility of those sites, but to evaluate turbine performance caused by changing geometrical variables and flow averaged characteristics. We adopted a velocity profile calculated from the averaged power in time for each site to evaluate the behaviour of the device in a single step. Keeping flow characteristics fixed at averaged conditions, facilitates the sensitivity analysis to geometrical parameters.

DMST sub-models
The DMST model is integrated with various sub-models to reproduce particular phenomena. One of them is the dynamic stall that characterizes CFTTs. During a blade's rotation, the angle of attack changes significantly and cyclically: the increase of the angle of attack leads to higher values of the lift coefficient compared with those reached in static conditions. The development of the so called Leading Edge Vortex (LEV) [17][18] near the leading edge of the hydrofoil, is responsible for this behaviour. This vortex increases the suction on the hydrofoil preventing and delaying stall. For further increase of the angle of attack, LEV moves toward the trailing edge: this causes the drop of the lift force. If the angle of attack starts decreasing, then a hysteresis occurs in the CL-α curve. The dynamic stall model, present in the DMST routine, is based on the one developed at the University of Pisa, and largely explained in [19]. The latter is based on experimental data obtained for an airfoil that is different from the one considered in this work, that is NACA 0018. For this reason it was necessary to adapt the model through a tuning process.
Other important phenomena requiring specific sub-models are: the streamlines deflection; the blade tip losses. For the implementation of these last two sub-models, as well as for the DMST calibration described in section 3, some 2D and 3D CFD simulations have been performed with the software ANSYS-Fluent v19 [20]. The grids for the 2D CFD simulations are unstructured and made of quad-elements, generated with the software Ansys-Icem by means of the "patch-dependent" technique. For 3D simulations the grids are multi-block structured with the addition of O-grids to thicken the distribution of cells in the areas of greatest interest and at the same time to improve their quality. Two grid levels are used to simulate the blade rotation via the sliding mesh method: a fixed sub-grid with the outer dimensions of the flow domain and a rotating sub-grid including the turbine blades. All around the blades the grid is very fine to make sure that y + at the walls stay below 0.4. To model the turbulence, the k-ω SST (Shear Stress Transport) is adopted [21]; this model is widely used in case of flow characterised by strong adverse pressure as typically happens for wind and tidal cross-flow turbines, especially when operating at low TSR. The algorithm for the velocity-pressure coupling is SIMPLEC. About the spatial discretization scheme, the Least Squares Cell-Based (LSCB) is set for gradient; pressure interpolation, turbulent kinetic energy and specific dissipation rate formulations are based on second order schemes. Temporal discretization is also based on a second order implicit method. The convergence criteria for each time-step is 1×10 −4 for the residuals of continuity, velocity components, turbulence kinetic energy and specific dissipation rate. When using a sliding mesh, to obtain a satisfactory numeric convergence the time-step should not be larger than the time required for advance the mobile interface by a distance corresponding to one cell thus, in order to consider the smallest cell at the interface, a time-step corresponding to 0.5° of revolution was adopted. The overall validation of the 2D CFD model has been described in a previous article [22] in comparison to the experimental data by Bravo et al. [23], whereas the validation of the 3D CFD model has been done versus the wind-tunnel velocity profiles achieved by Vergaerde et al. [24] at several locations downstream a two straight bladed turbine (details of the validation in [25]).
Streamlines deflection can be summarized as follows: although the flow can be approximated as straight, it tends to deviate when approaching the turbine. This because the turbine is seen as an obstacle from the flow, and it tends to escape laterally [26]. We improved this aspect in the DMST routine by applying deviation angle to the flow belonging to each stream tube. The correction depends only on the azimuthal position (identified by the angle θ) and the TSR. A more evident deviation occurs for higher values of TSR (as shown in figure 4). A high TSR value, indeed, makes the turbine less permeable to the flow. Figure 4 was obtained from CFD 2D simulations. Deviation angles, extracted from figure 4, are used to assign the deviation of streamlines by means of linear interpolation for each TSR-θ local values. Fig. 4. Deviation angle for σ1 at TSR 4 and 1,2 respectively a) and b) and for σ2 at TSR 3 and 1 respectively c) and d) The last phenomenon considered is the fluid dynamic losses that occur at blades' tips. Due to the finiteness of the wing, near blade's tip the flow tends to climb over the blade instead of following the hydrofoil profile. This causes a drop in lift force and turbine performance [12]. Instead of using the Prandtl-Glauert or Shen [27] formulation, we adopted a correction based on the power coefficient CP. The correction factor is k, and represents the ratio between the local CP along the blade and the CP of the turbine mid plane. The k factor is extracted by 3D simulations carried out with an undisturbed flow velocity of 1,75 m/s. The trend of the k factor is shown in figure 5.

Analysis guideline
The described methodology identified 10 high energy potential sites. These sites, being characterized by several depths, do not allow to compare all the 16 turbine geometries that we have hypothesized. It is worth noting, that the methodology has a general validity, and it can be used for different aims. For example, if the aim is to find the turbine geometry that maximizes annual energy production, the DMST must be applied to each of the 10 sites. In this case, we will be able to compare the best solutions from an energetical point of view, and choose between them taking into consideration investment costs and possibly environmental issues. However, this is not the purpose of this paper. We want to show the utility of comparing performance coming from different turbines. That is the reason why, in the following, we will consider only three sites which are 5, 6 and 10, chosen both for the energetical potential and the velocity profile shape along the vertical direction. In figure 6 we can see different velocity profiles: the red rectangle delimits the maximum depth occupied by the turbine. Site 10 shows the biggest velocity variation along blades' length. Fig. 6. Velocity profile along the vertical direction for some sites of interest Moreover, in this work we have fixed the optimal TSR of the turbine in the upper part of turbine blades' to privilege blades' sections invested by higher flow velocities. The optimal TSR is fixed at a reference height computed as follows: the height, where the power centre of gravity along z is located, is established as height weighted on power available at that height. Then, the result of the weighted average is enhanced arbitrarily of 1/3 (L-Hweighted), obtaining the definitive reference height. 1/3 (L-Hweighted) is a compromise between exploiting better the highest velocities in the upper part of the water column, and maintaining the optimal section far enough from tips effects. Indeed, if we consider a section too close to tips it is possible that hydrodynamic tips losses frustrate our intentions. This is particularly true for σ2 where tips losses are relevant, as can be seen from figure 5b. The rotational speed is calculated by considering the optimal TSR value and the undisturbed flow velocity at the reference height.

Model calibration
For what concern the dynamic stall sub-model, it was necessary to carry out a tuning process. As already mentioned indeed, the original sub-model was adapted for a different airfoil. Two tuning processes were made separately for each solidity, and compared to the relative 2D CFD simulations, carried out with ANSYS Fluent and an undisturbed flow speed equal to 1,75 m/s. Fig. 7. CP-TSR curve that compares the CFD trend with results of the tuning process for σ1 a) and for σ2 b) Figure 7a shows the trend for σ1 of the power coefficient at different TSR: the model overestimates the optimal TSR with respect to the CFD results, and a quick drop in performance can be seen at TSR lower than the optimal one. Whereas, the σ2 tuning process has brought to a curve that has a similar CFD trend. In this case, the optimal TSR value is well captured (figure 7b).
It is worth noting that, for both solidities, the model is sensitive to changes in the Reynolds number. Varying turbine geometries, or considering significantly different undisturbed flow velocity, can lead to changes in the CP-TSR curve. This is particularly true for σ1: figure 8a shows that the model leads to higher performance when Reynolds number increases, as expected. However, the optimal TSR value can vary with geometry. In figure  8a, we can see that the optimal value for the turbine "Area 50 AR 0,67" is 2,7, whereas for turbine "Area 50 AR 2" it is 2,8. This behaviour seems to affect low solidity turbines: indeed, high solidity turbines have always the optimal TSR at 1,4. This must be considered, and it is particularly important seen the aim of this work. We want to evaluate the behaviour of different turbine geometries by setting the optimal working TSR at the reference height. Figure 8 shows the 2D trends of performance for each turbine, obtained by running the DMST routine at an undisturbed flow velocity equal to that at the reference height for each case. We have considered site 5 velocity profile.

CP trend
In figure 9 we can see the power coefficient (calculated taking into account tip losses) for the 3 sites chosen and for all the geometrical turbine characteristics. Site 5 is the one with the best performance in each case, because of the higher energy content. For turbine with 25 m 2 frontal area (figure 9a), it can be observed that for both solidities an increase in AR, and therefore in turbine height, lead to a trend that is not monotonously increasing. Here we have three competing phenomena: the increase in AR leads to lower tip losses, but also lower Reynolds numbers, and turbines affected by a flow with a higher Δv (Δv is intend to be the difference in velocity module between the top and bottom of the turbine). A higher Δv means higher ΔTSR along the blade, and this causes a drop in performance.
A longer blade has a larger portion of the section working at TSR more and more distant from the optimal one. In figure 9a we can see that switching from AR 0,67 to AR 1,11 causes a reduction in performance. This is because a small increase in the AR leads to a decrease in tip losses that is negligible compared to the decline in performance, due to larger range variation in TSR along the blade and lower Reynolds numbers.  . "Net of tip losses" is the ratio between power output with tip losses and without tip losses. "Net of other losses" is the ratio between CP without tip losses and the CP without tip losses of the turbine with the same frontal area but AR 0,67. "Overall net" is the product of the previous two terms.
For the frontal area equal to 50 m 2 (figure 9b), we note that only for site 10 there is a decreasing trend of CP with increasing AR. This is due to the peculiarity of the velocity profile shape in site 10. Indeed, in this case the ΔTSR plays the most important role. For instance, turbine "Area 50 AR 1,55" with σ2, in site 5 has a relative ΔTSR that is the 17,6% with respect to the optimal TSR value 1,4, whereas in site 10 the same turbine has a relative ΔTSR of 42,5%. Site 10 is not deep enough for turbine "Area 50 AR 2", which is why it is not shown in figure 9b.
In figure 10 we can see the competition between the tip losses effect and the TSR and Reynolds number effect. For site 10 we report here the dimensionless power generated after tip losses curtailment (red columns) that enhance with the increase of AR. This is more evident for σ2 because tip losses are greater. We can also see that the tip losses have no linear trend with the increase in AR, but are more noticeable the more the AR increases. The blue columns represent the ratio between the CP and the CP of the turbine with same frontal area but AR 0,67 (both CP without tip losses). This allows evaluating the influence of other fluid dynamic losses due to TRS and Reynolds number changes. The green columns is the product of previous terms and is a sort of ratio between the CP and the CP without tip losses of the case with same frontal area and AR 0,67 (that became a kind of comparing term for all the other cases). Hence, green columns represent the ratio between CP and the maximum CP available for that site. Green columns have almost the same height for σ1, meaning that the competing phenomena balance out. Whereas, for σ2 green columns are higher where the relative CP in figure 9 is higher. Figure 9 shows also the relevant difference in power coefficient between the two solidities. σ1 always has the greatest power production, keeping fixed all the other conditions. The power coefficient at the optimal TSR from the tuning curve ( figure 7) shows a difference that amounts about 13% between the two solidities. However, in figure 9 we can see differences in performance for the two solidities that are about 20% and more. We must consider that turbines with σ2 are affected by higher tip losses (as can be observed also from figure 10), but this is not sufficient to clarify the behaviour. Indeed, if we consider performance coefficient obtained without applying the tip losses correction, we can see comparable percentage differences between performance for the two solidities. The explanation cannot be found in the quick drop in performance that characterize high solidity turbine. If we consider a fixed AR for the turbine, and the same site, turbines with different solidities are affected by the same percentage ΔTSR along blades (percentage ΔTSR relative to the optimal TSR for each solidity). Consequently, we have seen that a similar percentage ΔCP, relative to the optimal CP along blade can be observed. Therefore, the only explanation must be the influence of the Reynolds number. Indeed, all other conditions being equal, the change in solidity affects the Reynolds number. We already saw the importance of this parameter, as show in figure 8. If we consider the tuning curve in figure 7 and the curve turbine "Area 50 AR 0,67" in figure 8, we have a reduction in Re that is about 48% for each solidity. This because we have a high reduction in the undisturbed flow velocity, that goes from 1,75 m/s for the tuning process to about 0,82 m/s, used in figure 8. This causes a drop in performance more relevant for σ2 than σ1. So, σ2 seems more sensitive to variation in Re.
Therefore, from figure 9 emerges that for site 5 and 6 higher performance are recorded for turbine with frontal area 50 m 2 and AR 2 for each solidity, whereas for site 10 maximum performance results for area 50 m 2 and AR 0,67 for each solidity.

Power trend
It is useful to highlight that the highest power coefficients do not correspond to the maximum power production. Indeed, the power coefficient is a measurement of the device's performance relative to the power available. Since, as the turbine height increases the flow velocity decreases, this causes a reduction in the available power. In our study we kept the area value fixed (considering two values 25 m 2 and 50 m 2 ), and an increase in AR leads to taller but thinner turbine. In this way, we reduce the portion of area that works at the higher velocity and replace it with a portion that works with more and more low velocity as qualitatively shown in figure 11. The advantages coming from the increase in AR, that most of all concerns the reduction in fluid dynamics tips losses, are negligible compared to the loss in available power.
The trend of the extracted power per unit area (E) and the power coefficient are summarized in table 2 and 3 for σ1 and σ2 respectively. Probably, results can change if we consider geographical sites with more turbulent flow. In this case, we can imagine a velocity profile that is flatter compared to those in figure 6: it quickly reaches the undisturbed flow speed. In this way different turbine can be reached by a velocity profile that has no relevant gradient. In any case, this analysis highlighted the importance of considering devices in the operating condition. It reminds us to read with critical spirit literature results obtained with software that are not able to reproduce realistic operating conditions.

Conclusions
A MATLAB routine, based on DMST model, was run giving as input flow data obtained from 3D SHYFEM simulations. The routine was used for site assessment purposes, taking into consideration bathymetric constraints and the energetic potential. A sub-set of locations was selected for the energy content, for the depth of the site, and the peculiar velocity profile shape. We were interested in evaluating the turbine performance by changing some geometrical variables The results show that CP -AR trends (considering fixed all other geometrical variables) are influenced by competing phenomena: increase in AR means lower tip losses, but also lower Reynolds number, and greater TSR variation along blade. Depending also on site flow characteristics, we have found that site 5 and 6 reach higher performance with "Area 50 AR 2" for both solidities, whereas site 10 reaches higher performance with "Area 50 AR 0,67" for both solidities. The CP is a measurement of the turbine performance related to the available power. So it is not a parameter representative of the maximum power production. The shorter turbines produce the maximum power per unit area for both solidities, this means turbine "Area 25 AR 0,67". Moreover, to obtain the same overall frontal area of 50 m 2 , we can consider two paired counter rotating turbines. In this way we can exploit also other fluid dynamic mechanisms that furtherly enhance power production.
The turbine rotational speed was calculated in this work, by using the optimal TSR fixed at the reference height and the undisturbed flow velocity at that height. The so obtained rotational speed is not necessarily the one which maximizes performance. In future work, an iterative process can be adopt to predict the optimal rotational speed.
To conclude, the DMST allows to quickly identify suitable sites for tidal energy exploitation and to choose the best turbine geometry in terms of performance or maximum obtainable power. Hence, it is a useful tool for preliminary analyses, before further computationally burdensome full 3D CFD simulations.
Affiliations of authors should be typed in 9-point Times. They should be preceded by a numerical superscript corresponding to the same superscript after the name of the author concerned. Please ensure that affiliations are as full and complete as possible and include the country.