*Article* **Incorporation of a Non-Constant Thrust Force Coefficient to Assess Tidal-Stream Energy**

**Lilia Flores Mateos 1,2,3,\*,† and Michael Hartnett 1,2,3,†**


Received: 16 September 2019; Accepted: 25 October 2019; Published: 31 October 2019

**Abstract:** A novel method for modelling tidal-stream energy capture at the regional scale is used to evaluate the performance of two marine turbine arrays configured as a fence and a partial fence. These configurations were used to study bounded and unbounded flow scenarios, respectively. The method implemented uses turbine operating conditions (TOC) and the parametrisation of changes produced by power extraction within the turbine near-field to compute a non-constant thrust coefficient, and it is referred to as a momentum sink TOC. Additionally, the effects of using a shock-capture capability to evaluate the resource are studied by comparing the performance of a gradually varying flow (GVF) and a rapidly varying flow (RVF) solver. Tidal-stream energy assessment of bounded flow scenarios through a full fence configuration is better performed using a GVF solver, because the head drop is more accurately simulated; however, the solver underestimates velocity reductions due to power extraction. On the other hand, assessment of unbounded flow scenarios through a partial fence was better performed by the RVF solver. This scheme approximated the head drop and velocity reduction more accurately, thus suggesting that resource assessment with realistic turbine configurations requires the correct solution of the discontinuities produced in the tidal-stream by power extraction.

**Keywords:** tidal-stream energy; thrust force coefficient; momentum sink; unbounded flow; open channel flows; shock-capturing capability

#### **1. Introduction**

Resource characterisation is the initial step of any tidal energy project. However, conventional methodologies to assess tidal-stream power are based on infinite extent flow and the representation of a turbine as actuator disc [1–6]. This approach excludes the influence of boundary conditions such as free-surface and seabed in the power extraction analysis, which leads to the Lanchester-Betz–Joukowsky limit and indicates that in optimal conditions turbines can convert up to 59% (or 16/27) of kinetic energy into mechanical energy [7]. To better represent the aquatic medium of tidal turbines analytical approaches introduced features of a finite flow; firstly by introducing a rigid lid surface and lately by considering a deformable surface and a turbine downstream region where turbine wake mixing occurs. This analytical model is known as the linear momentum actuator disc in open channel flow (LMAD-OCH) [8]. Analysis of an actuator disc within a finite flow unveiled the importance of turbine bypass flow, blockage ratio (*B*), and the turbine downstream mixing region; where *B* indicates the fraction of a channel cross-sectional area occupied by the turbine.

A large number of potential regions for tidal power extraction are located in coastal areas; in these regions the use of finite flow is more appropriated as tidal streams are strongly influenced by the seabed and free-surface conditions. The consideration of these boundaries is relevant because they can significantly increase the amount of power extracted from the flow [9–11]. This is because in shallow waters the seabed and free surface constraint effect in the turbine is more significant and produces a stronger blockage effect, which in turn increases the maximum power extractable by the turbine. On the other hand, a recent methodology, which undertakes natural boundaries in the power extraction analysis, uses a rather expensive computational technique and constrains the head drop produced by power extraction. This methodology uses a rapidly varying flow solver and shock fitting techniques to implement a line sink of momentum [11,12]. The line sink of momentum is based on the LMAD-OCH analytical model, which enable to determine a relation between upstream and downstream depth-average flow velocities and depths as a function of the turbine operating conditions. This relation is given by the relative change of head drop across the turbine. Line sink of momentum approach uses an RVF solver to compute the rapid changes produced by the tidal turbine power extraction in the water depth and flow velocities. The authors of [11,12] used a Godunov-type scheme, which implements a shock-fitting technique. The momentum extracted by turbines is simulated by computing the modification of mass and energy fluxes across the line sink of momentum. The region occupied by a turbine array acts as an interface between upstream and downstream conditions of the flow. This interface represents the elevation and velocity discontinuities produced by power extraction. Computation of the discontinuity which initially separates upstream and downstream conditions is equivalent to solving a Riemann problem [13]; consequently, the Riemann solution indicates the flux through the discontinuity and therefore the shock propagation. To complete the line sink of momentum numerical representation, a condition on the head drop across the array is required [14]. This condition is given by the relative change of the head drop provided by the LMAD-OCH. This method was used to assess the tidal-stream potential of turbines configured as a fence in an idealised channel [11]. Such a configuration was used in the Pentland Firth to estimate maximum power extracted, 4.2 GW, and the extractable power at the sub-channels [15]. It was found that power availability varies according to the device operating conditions within the fences. A further refinement of the Pentland Firth estimation is given by [16]; they reported a more conservative upper limit of 1.9 GW based on a large, but a viable blockage ratio (*B* =0.4). Contrary to the line sink of momentum approach, in this research the momentum extracted by the turbines was simulated considering the turbine operating conditions and implementing a method that does not constrain the water elevation change. For this task the momentum sink TOC was used [17], the method is based in LMADT-OCH and computes the axial component of the thrust force exerted by turbines to the flow by solving a sink term in the momentum equations. Consequently, momentum sink TOC incorporates the natural constraints of the coastal tidal-stream, the operating conditions of the turbine, and does not condition the head drop produced by power extraction. In this research, the momentum sink TOC was incorporated in a (1) conventional and (2) an innovative numerical solver to assess the resource. The conventional scheme solves gradually varying flows, which do not experience strong spatial gradients and are characterised by small Froude numbers. On the other hand, the innovative scheme solves rapidly varying flows that enable the simulation of flows, which experience discontinuities such as hydraulic jumps, flood waves, and shock waves. Therefore, a GVF solver simulates a flow that experiences smoother and slower changes than RVF. The rapidly varying flows solver implemented in this research is a shock-capturing model that consists of the algebraic combination of first-order and second-order upwind schemes. In this way, the model uses a low-order scheme to simulate the regions where strong gradients exist and spurious numerical solutions are likely to appear. This solver constitutes an efficient approach for simulating energy capture in rapidly varying flows because the scheme is not required to solve the Riemann problem at each grid, where an array of turbines are defined, to compute the discontinuities produced in the flow due to power removal. Solving the Riemann problem represents an expensive procedure in computational terms [18]

Two turbine configurations were selected to assess the resource at a regional scale: A fence and a partial fence. Partial-fence configuration covers a selected part of a channel cross-section producing flow diversion. Part of the stream passes through the array, experiencing velocity reduction due to the momentum loss, and the rest bypass the array, presenting velocity intensification when circumventing the partial fence. This scenario is referred to as unbounded flow and two analytical theories describe power extraction with this type of flow: (i) The LMAD-OCH and (ii) two scales models, their main characteristics are described below. LMAD-OCH studies a partial fence as a long row [11], and in this way, upstream conditions can be assumed to be uniform, and marine turbines within the array can be represented by actuator discs. This approach estimates the momentum extracted by a partial fence and provides information on energy lost only within the turbine's near field region (*Lv*) [11], i.e., just downstream of the turbine. LMAD-OCH is also referred to as a single-scale model [19] because the quasi-inviscid flow assumption used in the theory allows important turbulent mixing to occur at a far downstream region (*Lh*) where the pressure equilibrates across the channel cross-section. The *Lv* region is smaller than the far downstream region (*Lv* < *Lh*). Within the *Lv* region the quasi-inviscid assumption of LMAD-OCH enables the estimation of turbine-scale wake mixing, which occurs just downstream of the turbine and it is not significant in comparison to the array-scale wake mixing. This quasi-inviscid model is consistent with three-dimensional actuator disc computations of mixing just downstream of the turbines [20]. Therefore, the LMAD-OCH theory enables the parametrisation of turbine-wake mixing within the near-field region, in this region it is assumed that elevation and velocity perturbations due to power extraction occur. This assumption captures vertical flow variations produced by horizontal mixing effects at the turbine scale [21]. In contrast to the LMAD-OCH theory, the two scale model studies energy extraction via a partial fence considering two scales of flow: (a) A turbine scale flow, which describes the flow around a single turbine and the wake of the turbine, studied by [8] and (b) a array scale, which is the large-scale flow around the turbine array and the array wake [22]. The two scales' separation implies that the flow bypassing the turbine occurs faster than the horizontal expansion of the flow bypassing the entire array, and analyses power extraction as two quasi-inviscid open channel flow problems. The two scales are linked via the upstream boundary condition of the turbine-scale flow because power extraction produces a flow diversion around the array and reduces the mass flux through the array and this mass flux reduction provides the upstream boundary to the turbine scale [23]. Initially, a single long turbine row and a rigid lid surface approximation which did not account for the effect of changes in water depth were considered [22]. Later, a better accounting for the interaction of the device- and array-scale flow events enabled the analytical and numerical study of short rows [20]. Recently, [24] included a deformable free-surface in the two-scale analysis by incorporating finite Froude (*Fr*) numbers in the analysis and consequently accounting for a deformable surface.

Analytical models provide an understanding of the physics involved in tidal-stream power extraction; however, numerical simulations are required to provide a more complete analysis of the effects of tidal energy extraction. In this context, the two-dimensional approach, the line sink of momentum was developed to assess tidal-stream resource with a long partial fence, through the solution of shallow water equations. The method was used to numerically estimate the force applied by a porous disc in scale experiments [25], and to identify the power available at an idealised coastal headland [26] and at Anglesey Skerries, off the Welsh coast [12]. On the other hand, attempts to numerically couple the turbine scale with the array scale are limited to three-dimensional simulations of small turbine arrays [27], which require high-fidelity turbine scale simulations; this methodology has not been applied to regional scale [21].

In this research, tidal-stream energy resource at regional scale was evaluated in two scenarios: bounded and unbounded flows, which were represented with a fence and a partial-fence turbine array configuration, respectively. The energy capture by the arrays was simulated with the momentum sink TOC method, which numerically implements the LMAD-OCH theory. To investigate the role of numerical solution schemes for simulating tidal energy capture, the momentum sink TOC was implemented in two hydrodynamic models that use a conventional (GVF) and up-to-date (RVF) solution scheme. The novelty of this research lays in the numerical representation of energy capture by arrays of turbines considering (i) a flow more representative of coastal regions, (ii) not conditioning of head drop across a turbine array, (iii) use of hydrodynamic models which are able to represent realistic coastal scenarios, and (iv) implementation of an efficient RVF solver that uses a shock-capturing scheme. The methodology proposed will enable the identification of a less computationally expensive numerical tool that provides a reliable evaluation of the resource in realistic scenarios.

#### **2. Methodology**

The assessment of tidal-stream resource at regional scale requires the specification of three aspects: A numerical model approach, the tidal energy extraction representation, and domain size [21]. These specifications follow.

#### *2.1. Modelling Approach*

Momentum sink TOC was implemented in two numerical hydrodynamic models based on the depth integrated velocity and solute transport (DIVAST) model [28]. DIVAST was developed to simulate hydrodynamic solute and sediment transport processes in rivers, estuaries, and coastal waters, and it incorporates a flooding and drying capability [29–31].

#### 2.1.1. Gradually Varying Flows

The first model uses an alternating direction implicit (ADI) methodology, solves two-dimensional shallow water equations (2D-SWEs), and simulates gradually varying flows, which are characterised by small Froude numbers. Hereafter ADI, represents the conventional strategy used to assess tidal-stream energy. The 2D-SWEs were used to describe the evolution of a tidal-stream and power extraction, the depth integrated continuity (Equation (1)) and *x*-component momentum (Equation (2), a similar expression was used for the *y*-component) equations to solve an inviscid flow, neglect Coriolis force, and to omit wind forcing. The inviscid character of the governing equations allows the neglection of vortical structures associated with the turbine's blade root and tip. The viscous terms omission is justified if bottom friction prevails over viscosity effects [32]. Moreover, a scaling analysis indicates that viscous terms can be disregarded if the horizontal scale of the domain is large and tidal currents change smoothly over the length of the domain [14]. Furthermore, viscous terms neglection has been considered in the numerical simulation of regional scale tidal-stream energy extraction [11,16].

$$\frac{\partial \mathcal{L}}{\partial t} + \frac{\partial q\_x}{\partial x} + \frac{\partial q\_y}{\partial y} = 0 \tag{1}$$

$$\frac{\partial q\_x}{\partial t} + \frac{\partial (\frac{\beta q\_x^2}{H})}{\partial x} + \frac{\partial (\frac{\beta q\_x q\_y}{H})}{\partial y} = -gH\frac{\partial \zeta}{\partial x} - \frac{gq\_x \sqrt{q\_x^2 + q\_y^2}}{H^2 \mathcal{C}\_\varepsilon^2} - F\_{T\ge} \tag{2}$$

where *qx* = *UH* and *qy* = *UH* indicate depth-integrated velocity flux component in the *x*- and *y*-direction; *t* stands for time and *β* is the momentum correction factor for non-uniform vertical velocity profile. The surface elevation change with respect to mean water depth *h* is represented by *ζ*; therefore, total water depth is defined as *H* = *h* + *ζ*. The depth-integrated turbulence model considers only the effects of the bed shear stress, which is a function of Chezy roughness coefficient (*Ce*) and gravity (*g*).

#### 2.1.2. Rapidly Varying Flows

The second model combines a MacCormack and a symmetric five point total variation diminishing (TVD) schemes to solve strong spatial gradients in the flow using an efficient shock-capturing method [30,33]. Henceforth, TVD allows the computation of discontinuities likely to appear in rapidly varying flows simulation without treating the discontinuity as a Riemann problem. To simulate flows which experience discontinuities TVD solves the conservative form of 2D-SWEs, [17] describe this form of the governing equations in detail. Some examples of this kind of flows are hydraulic jumps, storm surges, flood waves, and shock waves.

#### *2.2. Numerical Methods*

ADI and TVD use a finite difference spatial discretisation to approximate a solution of the governing equations. The models discretise the equations onto a square structured grid, where the edges of the grid cells are oriented parallel to the Cartesian coordinates [34,35].

The GVF solver uses an ADI's semi-implicit finite difference scheme which splits a single time-step solution into two time steps [36]. The final finite difference equations for each half time step (HFDT) are solved using the method of Gaussian elimination and back substitution [34]. Meanwhile, the RVF solver uses an explicit scheme to approximate the solution of the governing equations at the current time and convergence to the solution is conditioned to a maximum time step. RVF solver uses an algebraic combination of the first-order and second-order upwind schemes, where a kind of artificial viscosity is included to smooth the solution close to the shock [30,37]. The proportion of the contribution of each scheme is adjusted depending on the nature of the flow; if the solution is sub-critical (smooth), the second-order scheme is implemented otherwise, if the solution is trans- or super-critical, the lower-order scheme is implemented [30].

The stability of the numerical models solution was assured by satisfying Courant-Fredrichs-Lewy condition. To identify the most efficient spatial and temporal resolution a range of time steps and grid sizes were tested. The elected spatial discretisation was Δ*X* = 150 m in both models. The semi-implicit character of ADI allowed a relatively large time step (Δ*t* = 12 s), but the explicit scheme of TVD did not tolerate time steps larger than 1.50 s.

#### *2.3. Turbine Representation*

The method used to simulate energy capture by marine turbines is momentum sink TOC and it implements numerically the LMADT-OCH theory using a sink approach. The mathematical representation of energy capture was introduced as a sink term *F*-*<sup>T</sup>* in the momentum equations (see *x*-component, Equation (2)). This term represents the turbine thrust force (*T*) responsible for tidal stream momentum loss (Equation (3)). The thrust force is a function of the turbine cross-sectional area *<sup>A</sup>*, thrust coefficient *CT*, and depth-average velocity −→*<sup>U</sup>* .

$$T = \frac{1}{2} \rho \,\overrightarrow{\mathcal{U}}^2 A \mathcal{C}\_T \tag{3}$$

Figure 1a illustrates the cross section of the turbine area *A*, indicated by a dotted-red line, and the thrust force exerted on the flow. The components of the thrust force are given by:

$$F\_{Tx} = T|\sin\theta|\tag{4}$$

$$F\mathbf{r}\_{\mathbf{\bar{y}}} = T|\cos\theta|\tag{5}$$

where it is assumed that the thrust force makes an angle *θ* with the *y*-axis. Numerically, the sink term is the thrust force per unit-grid and per unit-mass. The components of the thrust force introduced as sink terms in the momentum equations and solved by the models were:

$$\begin{aligned} F\_{Tx} &= \frac{\mathbf{T}}{\Delta x \Delta y} |\sin \theta| = \frac{1}{\Delta x \Delta y} \frac{1}{2} \rho A\_x C\_T U^2 \\ F\_{Ty} &= \frac{\mathbf{T}}{\Delta x \Delta y} |\cos \theta| = \frac{1}{\Delta x \Delta y} \frac{1}{2} \rho A\_y C\_T V^2 \end{aligned} \tag{6}$$

Furthermore, the thrust force computed with momentum sink TOC is a function of the thrust force coefficient (Equation (7)), which depends on the specification of wake-induction (*α*4) factor and the bypass flow coefficient (*β*4).

$$C\_T = \beta\_4^{\cdot 2} - \alpha\_4^{\cdot 2} \tag{7}$$

LMADT-OCH theory allows to relate turbine operating conditions to momentum captured by the turbine and to link changes produced by power extraction to thrust forces within the turbine's near field length *Lv* [14]. LMADT-OCH analytical model assumes that the upstream flow passes through the fence, mixes, and returns to a vertical profile similar to that upstream over a length *Lv*. A sub-critical tidal-stream is assumed and energy capture produces a tidal flow division into core-flow and bypass-flow components, a head drop across the array (Δ*h*), and energy loss by turbine wake mixing. The core-flow refers to the stream that passes through the turbine, which experiences a velocity reduction due to energy extraction, represented by the turbine velocity coefficient *α*2. Core-flow presents a further velocity reduction downstream of the turbine due to turbine wake mixing dissipation, denoted by the wake induction factor *α*4. On the other hand, the bypass-flow circumvents the turbine and presents a velocity magnitude intensification denoted by the bypass induction factor *β*4. Based on LMADT-OCH energy capture within *Lv* region is parameterised by: bypass induction factor (*β*4), wake induction factor, turbine velocity coefficient (*α*2), and water depth drop across the array (Δ*h*).

Calculation of thrust force required: (i) election of turbine operating conditions, (ii) identification of upstream conditions of the flow, (iii) and parametisation of changes produced by power extraction within the turbine near-field region. These requirements are sketched in Figure 1a. Turbine operating conditions were defined by the blockage ratio *B* and wake induction factor *α*4. Factor *α*<sup>4</sup> is an indicator of velocity rate reduction due to wake mixing at the turbine scale [38]. In ADI and TVD models, the blockage ratio implemented is the ratio of the cumulative turbine area per cell-grid (*A*) over the grid-cell cross-section area (*H* Δ*X*) [17]. Parameterisation of flows' depth and velocity changes within *Lv* required the specification of *α*4, *B*, and upstream Froude number (*Fr* = *U*/ *gH*). These parameters enabled the calculation of *β*4, which is the physical admissible root of a quartic polynomial [10]. An eigenvalue method was used to find the roots of the polynomial [17]. Finally, is possible to calculate *CT*(*α*4, *β*4) and thrust force (Figure 1b).

**Figure 1.** Two-dimensional representation of the thrust force (*F*- *<sup>T</sup>*) exerted by the turbine to the incident tidal-stream (**a**) and the momentum sink turbine operating conditions (TOC) computation procedure (**b**).

#### *2.4. Tidal Channel and Turbine Array*

A tidal fence and a partial-fence configuration were implemented in a large domain, which corresponded to a large channel connecting two large basins (Figure 2). The channel is relatively narrow and long with 12 km length (*L*) and 3 km wide (*W*) and it is characterised by an aspect ratio *W*/*L* = 0.25. The basin extension is 4.5 times the length of the channel. The domain configuration was selected based on previous studies of tidal-stream energy assessment [11]. The channel was forced with a standing wave with semi-diurnal frequency *M*<sup>2</sup> on the western boundary (Equation (8)). The amplitude of the forcing was introduced as a ramped-up wave over two *M*<sup>2</sup> tidal periods (Equation (9)); the gradual

introduction of the forcing removed the initial numerical noise in the simulation and helped to reach a time varying steady state faster. This state was reached when tidal circulation became regular.

$$\mathcal{L} = A(t)\cos\left(\omega t\right)\cos\left(k\mathbf{x}\right) \tag{8}$$

$$A(t) = A\left(-2\left(\frac{t}{2T\_{M\_2}}\right)^3 + 3\left(\frac{t}{2T\_{M\_2}}\right)^2\right), t < 2T\_{M\_2} \tag{9}$$

At the eastern boundary, water elevation was set to constant and velocity was set to zero. At the walls, the boundaries satisfied a no-slip condition. The water depth is relatively shallow, it was set to 40 m = 2.5 *d*, where *d* = 16 m is the turbine diameter. A small and constant bottom friction, defined as a function of a non-dimensional drag coefficient, was implemented (*Cd* = 0.0025). This coefficient has been shown to best reproduce field measurements of velocity and elevation of *M*<sup>2</sup> tidal currents in the vicinity and far-field of Rathlin Sound [39].

Simulations started from quiescent initial conditions and a steady periodic flow was reached after the second tidal period. The simulation time was 50 h, equivalent to four *M*<sup>2</sup> tidal periods (4 *TM*<sup>2</sup> ). It is worth mentioning that evaluations of 8 *TM*<sup>2</sup> were consistent with conditions reported at the fourth tidal cycle. The consistent boundary conditions and forcing implementation between the models were assured by obtaining consistent hydrodynamic conditions at the natural state. This state refers to marine turbines omission and nil energy extraction.

Turbine configurations were deployed in the middle of the channel (Figure 2). The fence corresponds to an array of turbines distributed as a single row that fully extends across the channel cross-section. On the other hand, the partial fence covers 40% of channel cross-section; the length of the array is 1200 m defined over eight grid-cells, where each grid-cell contains a cluster of turbines.

**Figure 2.** Tidal channel model domain. A channel zoom out (upper-right location) shows the array of turbines deployed. Forcing implemented is represented by the ramping up elevation in magenta.

#### **3. Resource Assessment**

Though a significant amount of power could be extracted if high blockage ratios were used, this procedure presents adverse effects such as (i) significant flow rate reduction [40], (ii) tidal hydrodynamics affectation, and (iii) mixing and transport processes impact at turbine scale and regional scale [6,21,41,42]. As a result of this evidence, evaluation of tidal-stream resource is performed within realistic blockage ratios (0 < *B* ≤ 0.4) as a post-processing stage.

Indicators used to assess the energy resource were the following metrics: Total power extracted, power dissipated by turbine wake mixing, and power available for electrical generation. The procedure used to calculate them follows. The initial step was the election of turbine operating conditions and calculation of turbine velocity coefficient (Equation (10)). Subsequently, thrust (Equation (7)) and power (Equation (12)) coefficients associated to the turbine were estimated. The next step was identification of head drops across the array and turbine efficiency (Equation (13)). These

two parameters enable estimation of total power extracted (Equation (14)), power dissipated by turbine wake mixing (Equation (15)), and power available for electrical generation (Equations (11) and (16)).

$$\alpha\_2 = \frac{2(\beta\_4 + \alpha\_4) - (\beta\_4 - 1)^3 (B\beta\_4^{\
\alpha\_4} - B\beta\_4 \alpha\_4)^{-1}}{4 + (\beta\_4^{\
\alpha\_4} - 1)(\alpha\_4 \beta\_4)^{-1}}.\tag{10}$$

$$P = \frac{1}{2}\rho l l^3 A \mathcal{C}\_P;$$

$$\mathbb{C}\_{P} = a\_{2} (\beta\_{4}^{2} - a\_{4}^{2}) \\ \tag{12}$$
 
$$\left( \begin{array}{c} 1 \ \Delta h \end{array} \right) \tag{13}$$

$$
\eta \approx a\_2 \left( 1 - \frac{1}{2} \frac{\Delta h}{h} \right) \tag{13}
$$

$$Pr = \rho g l I \frac{A}{B} \Delta h \left(1 - F\_r^2 \frac{1 - \Delta h / 2h}{(1 - \Delta h / h)^2} \right) \tag{14}$$

$$P\_W = P\_T(1 - \eta) \tag{15}$$

$$P\_\* = \eta P\_T \tag{16}$$

The operating conditions of the turbine and the turbine configurations evaluated with ADI and TVD models are specified in the four scenarios reported in Table 1.


**Table 1.** Scenarios simulated and initial parameters specification.

#### *3.1. Thrust and Power Coefficients Time Series*

Calculation of *CP* and *CT* required specification of the wake induction factor, which depends on the location of the turbine within the array and the local tidal dynamics [26] as both affect turbine wake mixing. However, to maximise power available for electrical generation the downstream core flow (*α*4*U*) was adjusted by setting the wake induction factor to *α*<sup>4</sup> = 1/3 [43,44]. A constant *α*<sup>4</sup> was used during the simulation; however, a slightly better power available is obtained if a variable wake induction factor is used [12].

Thrust and power coefficients are functions of the upstream conditions through the Froude number used in *β*<sup>4</sup> calculation; therefore, *CT* and *CP* are modulated by the tidal cycle. Thrust and power coefficients time series were obtained from both models using a fence configuration constituted by turbines with a relatively large blockage ratio (*B* = 0.3). The time series from an *M*<sup>2</sup> tidal cycle and taken from the middle of the channel, where the array was deployed, are shown in Figure 3a,b. *CT* and *CP* time series present two maximum values that correspond to maximum flood tide and maximum ebb tide velocities. Comparison between the solutions obtained from the models indicate that ADI produces slightly higher magnitudes of *CP* (1.21 < *CP* < 1.238) and *CT* (2.38 < *CT* < 2.43) than TVD's 1.21 < *CP* < 1.235 and 2.36 < *CT* < 2.40. Time-average of *CT* and *CP* obtained from ADI and TVD are consistent with LMADT-OCH analytical model and they indicate upper limits for these coefficients.

Additionally, time series obtained with ADI present a phase delay evident in the maximum values. ADI's phase delay was evident at stream-wise velocity component of *M*<sup>2</sup> tidal current at natural state. This behaviour was consistent along the domain, i.e., it was not restricted to the channel, suggesting that the delay may be related to the different solution scheme employed by the models. However, further understanding of *M*<sup>2</sup> phase simulated by both models would require examination of field observations of free surface elevation and depth-averaged flow velocity. It is worth noticing that

the parameters to be presented in the following sections refer to values averaged over a tidal period; therefore, the phase delay reported by ADI time series does not influence the results presented.

The larger magnitude of *CT* reported by ADI is associated with the smaller affectation of the momentum extraction on the bypass velocity reduction (represented by *β*4). ADI underestimation of velocity reduction due to power extraction was reported by [45]. In the case of *CP*, the turbine velocity coefficient reduces the difference between the magnitudes obtained from the models.

**Figure 3.** Time-series of *CP* (**a**) and *CT* (**b**) coefficients obtained from alternating direction implicit (ADI) (continuous-line) and total variation diminishing (TVD) (dash-line) schemes for *B* = 0.3.

#### *3.2. Upstream Velocities*

Contrary to a tidal fence configuration, where energy extraction produce uniform conditions along the array and a middle-fence cell capture this dynamic, the implementation of a partial fence produces a scenario where energy extraction occurs along a limited section of the cross-channel, indicating the existence of array-bypass flows. To identify the effects of a partial fence on the upstream velocities, the time-averaged stream-wise component of velocity was calculated at each grid cell that constitutes the partial fence. Upstream velocities normalised by natural state velocities and obtained for increasing values of blockage ratio computed with ADI and TVD are presented in Figure 4. Upstream velocities along the fence tend to be uniform for small blockage ratios (*B* ≤ 0.15), but further increase of *B* influence differently the grid-cells located at the middle and at the edges of the partial fence. Large values of *B* produce higher velocity decreases in the middle of the array, meanwhile the edges exhibit smaller velocity reduction. A comparison of results obtained from both models shows that upstream velocities from ADI present smaller reductions with increasing *B* than TVD. In addition, for small *B*, velocities at the edges of the partial fence simulated with ADI show a slight magnitude enhance. This suggests that the faster array-bypass flow computed with ADI model is influencing the velocities at the partial-fence edges. The smaller velocity reduction presented in ADI's velocity solutions is consistent with [45]. They indicated that a gradual varying flow scheme under-estimate the velocity decrease produced by momentum extraction; meanwhile, a rapidly varying flow scheme approximates more accurately this reduction due to its capacity to solve the velocity spatial gradients generated during energy extraction.

**Figure 4.** Blockage ratio effect on normalised upstream velocities along a partial fence, ADI (**a**) and TVD (**b**). Natural state flow condition is indicated by NT.

#### *3.3. Head Drop and Turbine Efficiency*

Estimation of turbine efficiency required calculation of head drop across the array. Δ*h* was calculated using an analytical and numerical approach. Analytical solutions were obtained by solving a cubic polynomial, which represents a condition on the head drop across the array obtained from LMADT-OCH [8]. Coefficients of the polynomial are a function of blockage ratio, thrust coefficient, and Froude number. *CT* and *Fr* were obtained from TVD simulations, as this scheme approximates the velocity reduction more accurately for a fence scenario [45]. Meanwhile, numerical solutions correspond to head drops obtained from depth differences across the array at upstream and downstream locations. These locations correspond to the neighbouring upstream and downstream cell. In the case of a 16 m diameter device, the mean axial upstream velocities are recovered to about 80% [46,47] within 160 m; therefore, the grid size used in this research (Δ*X* = 150 m) is of the order of the turbine wake length suggesting that the neighbouring cell approximates the upstream conditions of the flow. Maximum head drops within a tidal cycle were calculated for increasing *B*, Δ*hmax* obtained analytically and numerically using a fence (Figure 5a) and a partial fence (Figure 5b) were compared. For bounded-flows, analytical and numerical solutions are similar; however, a GVF scheme (Figure 5a.1) produces a solution more consistent with the analytical solution (Figure 5a). In the case of unbounded flows, the analytical solution of LMAD-OCH indicates a homogeneous water drop along the partial fence for given blockage ratio (Figure 5b). This solution is consistent with the 1D assumption of the theoretical model [9,43]. Conversely, numerical solutions indicate non-uniform head drops along the array for increasing *B* values. Within the partial fence, maximum values of head drop are found in the middle cells, while smaller depth change are exhibited towards the edge of the array. Comparison of Δ*hmax* obtained from TVD and ADI models for a partial fence indicate that TVD simulates larger head drops (Figure 5b.2). Additionally, for small *B*, Δ*hmax* at the middle-cells of the array obtained from TVD are similar to analytical solutions; however, for *B* ≥ 0.30, the head drop magnitudes become smaller than the analytical solution. These results indicate that rapidly varying flow solver provides a reasonable insight into the effect of small blockage ratios on the head drop across the partial fence, where Δ*h* is influenced by the bypass flow producing a head drop reduction at the edges of the partial fence.

Turbine efficiencies were estimated for increasing blockage ratios. Analytical and numerical solutions of *η* for flows characterised by small Froude numbers were calculated, as the solutions were found to be consistent only numerical solutions are reported. These solutions are functions of the head drops obtained with ADI and TVD. Time-averaged and array-averaged turbine efficiencies calculated for both configurations: fence and partial fence are presented in Figure 6a,b, respectively. Turbine efficiencies were plotted against the time-averaged, array-averaged thrust coefficients. Increasing *B* values indicate a gradual efficiency reduction but thrust coefficient augmentation. Maximum turbine efficiency found (*η* = 0.63) correspond to the smaller blockage ratio tested (*B* = 0.05) and this value is associated to *CT* = 1.0; meanwhile, blockage ratio increases to *B* = 0.4 produce a turbine-efficiency reduction (*η* = 0.47) and thrust coefficient increase (*CT* = 3.5). Efficiency decrease with *B* augmentation is explained by the larger energy dissipated due to turbine-wake mixing [48]. On the other hand, the analytical study of further blockage ratio increase (*B* > 0.4) indicates a slow turbine-efficiency decrease rate [14]. This trend indicates that for a very large *B*, subsequent dimension augmentation produces a small efficiency reduction; such a scenario is attractive; however, large blockage ratios are not a practical option.

**Figure 5.** Maximum head drops across a fence (**a**) and partial fence (**b**); analytical (filled markers) and numerical solutions (un-filled markers) obtained from ADI (\*.1) and TVD (\*.2) models.

Comparison of ADI and TVD solutions from a partial fence shows similar results; however, slightly higher efficiency values are reported with TVD. As the efficiency is a function of Δ*h*, higher *η*'s are explained by the larger head drops obtained with RVF scheme. Furthermore, comparison of the solutions obtained from both configurations indicate a high degree of similarity. This consistency suggest that the partial fence is long enough to resemble conditions similar to a fence at the turbine's near-field region. It is worth mentioning that turbine efficiencies calculated for a partial fence do not capture the power dissipated by the array-wake mixing, because LMAD-OCH theory only describes the energy extraction dynamic within the turbine's near-field region [20,23]. In this region, wake mixing produced by individual turbines occurs and flow's uniform conditions are recovered; therefore, wake mixing generated by the array is not included in the analysis.

**Figure 6.** Turbine efficiency and thrust coefficient for partial fence (**a**) and fence (**b**) configurations obtained from ADI and TVD.

#### *3.4. Power Metrics*

To compare the performance of a partial fence with a fence, time-averaged and array-averaged power metrics were obtained from both configurations. Total power extracted estimated with ADI and TVD are presented in Figure 7(a.1,a.2), respectively. *PT* is a function of water drop and the increase of blockage ratio produces a head drop augmentation which is responsible of larger values of total power extracted. For bounded flows (F), [45] indicated that *PT* is accurately approximated by an ADI model because a GVF solver approximates better the head drops across a fence. For unbounded flows scenario (PF), *PT* magnitudes obtained from ADI and TVD are smaller than the fence. This result is explained by the head drop reported by the models for the partial fence; smaller Δ*h* magnitudes were obtained for this configuration due to the reduced head drops at the edges of the array. The comparison of *PT* obtained by the models for the partial fence indicate that TVD scheme produces larger magnitudes of total power extracted. These solutions are more reliable as TVD approximates better the head drop across a partial fence. The power dissipated by turbine-wake mixing estimated with both ADI and TVD models are shown in Figure 7(b.1,b.2), respectively. *PW* depends on the turbine efficiency and the total power extracted, which in turn is strongly influenced by the head drop across the turbine array. This dependency indicates that *PW* solutions are also explained by head drop values. For bounded flows, *PW* values obtained with ADI are reliable solutions [45]. In the case of unbounded flows scenario, *PW* magnitudes obtained from the models are slightly smaller than the fence. This result is expected as smaller Δ*hmax* were obtained for the partial fence. The comparison of *PW* obtained from the models for the partial fence, indicate that TVD produces larger magnitudes of power lost. TVD solutions are more accurate due to the better approximation of the head drop for unbounded flows provided by the model.

Regarding the power removed by the turbine, two metrics were calculated. First, in function of turbine efficiency *P*<sup>∗</sup> and second, in terms of cubic velocity *P*. These metrics are expected to be consistent as they describe the amount of power available for electricity generation.

Power available in terms of turbine efficiency estimated with ADI and TVD are presented in Figure 8(a.1,a.2), respectively. *P*<sup>∗</sup> is a function of total power extracted and consequently depend on the head drop. For bounded flows, *P*<sup>∗</sup> is better approximated by an ADI scheme; however, for unbounded flows scenario, this metric is better simulated by an TVD scheme.

On the other hand, power available in terms of cubic velocity calculated with ADI and TVD are presented in Figure 8(b.1,b.2), respectively. *P* magnitude depends on the upstream velocity and for bounded flows, [45] reported the incapacity of ADI model to satisfactorily approximate the velocity reduction due to power extraction. This underestimation of velocity decrease for increasing values of blockage ratio explains the large values of *P* reported by ADI for both configurations. Meanwhile, *P* obtained with TVD for a fence and a partial-fence configuration present magnitudes similar to *P*<sup>∗</sup> obtained from both ADI and TVD models. This result indicates that TVD scheme better approximates the consistency between *P* and *P*<sup>∗</sup> for both configurations. ADI fails to represent this consistency as ADI underestimate the velocity reduction due to power extraction and consequently report higher values of Power than TVD. On the other hand, TVD simulates better this velocity reduction, so powers in function of velocity are better solved by TVD and these values are more consistent with the power in function of turbine efficiency.

**Figure 7.** Effect of blockage ratio on *PT* (**a**) and *PW* (**b**) for a fence (dash-line) and partial-fence (continuous-line) configuration. Solutions obtained from ADI (\*.1) and TVD (\*.2).

**Figure 8.** Effect of blockage ratio on *P*<sup>∗</sup> (**a**) and *P* (**b**) for a fence (dash-line) and partial-fence (continuous-line) configuration. Solutions obtained from ADI (\*.1) and TVD (\*.2).

#### **4. Discussion and Conclusions**

This research implemented momentum sink TOC method in two hydrodynamic models that can solve realistic coastal scenarios to simulate marine turbine configurations and to perform tidal-stream resource assessment. Contrary to the conventional consideration of constant thrust coefficient [4,6,49], the use of turbine operating conditions to calculate the thrust force, imparted by the turbine on the stream, enabled the calculation of non-constant thrust and power coefficients. These coefficients undertake the tidal regime specific to a potential site and allow the use of tailored values of *CT* and *CP* to assess the resource, therefore providing a more realistic numerical representation of the turbine.

In contrast to [11], the head drop across the array used to assess the resource was obtained from the depth difference between the array upstream and downstream location. The use of a no-fixed head drop in a partial fence enabled the identification of non-uniform conditions along the array. This configuration represented the simulation of an unbounded flow, which favours the existence of

two regions: array-bypass flow and array-wake flow. These regions describe the velocity intensification when bypassing the array, and the significant velocity diminution downstream of the array due to array-wake mixing. Flow conditions along the partial fence are not uniform as array-bypass flows influence the upstream velocity and water depth at the edges of the array. These features illustrates the effect of a partial fence on the free flow as reported by [23,50] and indicate the advantage of using numerical tools to evaluate the tidal-stream resource and complement analytic models.

In terms of turbine efficiency, both turbine configurations showed consistent efficiencies. The similarity of the results could be explained by the use of the same wake-induction factor (*α*<sup>4</sup> = 1/3), which is used for the turbine velocity coefficient (*α*2) calculation [45]. The turbine efficiency used in this research considers the turbine-wake mixing within the near-field region. Over the near-field length scale, it is assumed that flow passing through the fence and partial fence mixes and leads to a smooth vertical profile similar to upstream profiles [43]. This assumption captures vertical flow variations produced by horizontal mixing effects at the turbine scale [21]. The models compute the turbulence generated by turbine-wake mixing correctly by implementing (i) turbulence induced by the bottom friction and (ii) inviscid flow assumption [51].

This research assesses upper limits for tidal-stream resource assessment in a semi-narrow tidal channel using the *M*<sup>2</sup> tidal constituent. Additionally, evaluation of two turbine configurations allowed the identification of the numerical scheme that performs more accurate tidal-stream energy resource evaluation. Opposite to conventional methodologies to assess the resource, where tidal energy extracted by marine turbines was approximated with a bottom roughness [40,52,53] and a momentum sink [6,47,49,54], the numerical implementation of LMADT-OCH enabled the calculation of the turbine efficiency. This parameter provided a distinction between total power extracted, power available for electrical energy generation, and power dissipated by turbine-wake mixing. Evaluation of the conventional scheme to assess tidal-stream (ADI model which solves GVF) and a scheme with shock-capturing capability (TVD model which solves RVF), indicate that unbounded flow scenarios are better assessed with a TVD scheme. The better approximation of the head drops and velocities reduction by a TVD model for a partial-fence scenario implies that a rapidly varying flow solver is required to calculate an upper bound for tidal-energy extraction when a realistic configurations of turbines is used. The conclusions of this work are as follows:


**Author Contributions:** Conceptualization, L.F.M. and M.H.; methodology, L.F.M.; validation, L.F.M.; formal analysis, L.F.M.; resources, M.H.; writing–original draft preparation, L.F.M.; writing–review and editing, M.H.; visualization, L.F.M.; funding acquisition, M.H.

**Funding:** This material is based upon works supported by the Science Foundation Ireland under Grant no. 12/RC/2302 through MaREI, the National Centre for Marine Renewable Energy Ireland.

**Acknowledgments:** The authors thank the Irish Centre for High-End Computing for the provision of computational facilities and support.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


c 2019 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 (http://creativecommons.org/licenses/by/4.0/).
