*Article* **Development of a Hydrokinetic Turbines Backwater Prediction Model for Inland Flow through Validated CFD Models**

**Chantel Monica Niebuhr 1,\*, Craig Hill <sup>2</sup> , Marco Van Dijk <sup>1</sup> and Lelanie Smith <sup>3</sup>**


**Abstract:** Hydrokinetic turbine deployment in inland water reticulation systems such as irrigation canals has potential for future renewable energy development. Although research and development analysing the hydrodynamic effects of these turbines in tidal applications has been carried out, inland canal system applications with spatial constraints leading to possible blockage and backwater effects resulting from turbine deployment have not been considered. Some attempts have been made to develop backwater models, but these were site-specific and performed under constant operational conditions. Therefore, the aim of this work was to develop a generic and simplified method for calculating the backwater effect of HK turbines in inland systems. An analytical backwater approximation based on assumptions of performance metrics and inflow conditions was tested using validated computational fluid dynamics (CFD) models. For detailed prediction of the turbine effect on the flow field, CFD models based on Reynolds-averaged Navier–Stokes equations with Reynolds stress closure models were employed. Additionally, a multiphase model was validated through experimental results to capture the water surface profile and backwater effect with reasonable accuracy. The developed analytical backwater model showed good correlation with the experimental results. The model's energy-based approach provides a simplified tool that is easily incorporated into simple backwater approximations, while also allowing the inclusion of retaining structures as additional blockages. The model utilizes only the flow velocity and the thrust coefficient, providing a useful tool for first-order analysis of the backwater from the deployment of inland turbine systems.

**Keywords:** hydrokinetic; computational fluid dynamics; backwater; inland hydrokinetic; axial flow turbines

#### **1. Introduction**

Research and development of hydrokinetic (HK) devices in canal systems is increasing in popularity due to increasing electricity costs and the drive towards finding renewable energy sources with unconventional applications [1–3]. Although most development has focussed on tidal applications, multiple opportunities exist for the deployment of HK systems within inland water infrastructure (e.g., canal systems) [1]. However, the placement of such a device can have significant water level and hydrodynamic energy loss effects [4].

Prediction of the hydrodynamic effects of hydrokinetic turbines in canal systems remains an important pre-development objective. Due to the nature of canal design, these systems usually have flat slopes and subcritical flow regimes. Therefore, the analysis of backwater effects from blockages is critical for the prevention of flooding and water loss. This is especially important in array schemes where the cumulative effect of multiple devices can exceed the top of the channel and cause it to overtop.

**Citation:** Niebuhr, C.M.; Hill, C.; Van Dijk, M.; Smith, L. Development of a Hydrokinetic Turbines Backwater Prediction Model for Inland Flow through Validated CFD Models. *Processes* **2022**, *10*, 1310. https:// doi.org/10.3390/pr10071310

Academic Editors: Santiago Lain and Omar Dario Lopez Mejia

Received: 23 May 2022 Accepted: 20 June 2022 Published: 4 July 2022

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

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

Generally, backwater calculations utilize a blockage size (e.g., typical backwards-facing step, weir or pier shape) or energy loss function (quantified energy losses) to predict the backwater effect. Due to the novelty of HK energy, a streamlined procedure for determining the backwater effect has not yet been determined. This may be attributed to the variability of turbine types, operational conditions, and efficiencies, all of which result in a different effective blockage.

Previous studies have investigated the hydrodynamic effects of horizontal-axis hydrokinetic turbines (HAHTs) both experimentally (e.g., in canals [5], investigating boundary layers [6] and varying Reynolds numbers [7]) and computationally (e.g., both CFD applications in [4,8]). However, most of the studies are performed under constant operational conditions, and are site-specific (e.g., three-bladed [9] and two-bladed [9,10] turbines under optimal conditions). Some attempts have been made in the past to develop backwater models using roughness values [4] or analytical relationships [11], but the lack of experimental results over a range of turbine designs and operational conditions has resulted in site-specific models. Additionally, use of the models without in-depth knowledge of input variables limits their utility.

This study aimed to develop a simplified method for calculating the backwater effect of HK turbines in canal systems. An analytical approximation based on assumptions of performance metrics and inflow conditions was tested using validated computational fluid dynamics (CFD) models. These models allow a larger dataset and, thus, validation of the analytical approximation recommended.

#### **2. Background**

Placement of an HK device extracting energy in confined flow may affect water surfaces and water surface profiles. This is especially true in array schemes, and must be considered to ensure that the clearance between the rotor blades and the surface is sufficient. Myers and Bahaj [12] investigated a 1:30 scale HAHT model (rotor diameter, *dT* = 0.4 m), and observed a clear difference in the water surface once energy was extracted. Water depths increased immediately upstream of the rotor and decreased downstream for about 2 *dT*. Details of the water surface profiles can be seen in Figure 1. The results observed a standing wave 7–8 *dT* downstream (it should be noted that this was for the high-freestream-velocity case).

**Figure 1.** Water surface profile through a scaled turbine operating at 2 different velocities, compared to the no-energy-extraction stage [12].

Prediction of such occurrences requires understanding of the specific energy and flow regime (Froude number) to predict flow behaviour after HK deployment; here, a Froude number based on turbine diameter (*FrD* = √*<sup>U</sup> gdt* , where *U* is the mean velocity and *dt* is the turbine diameter) may be more useful, which has also been found to govern the free-surface effects [13].

The flow effects observed can be explained using the specific energy of the flow section, which is a function of water depth and velocity. When ignoring friction losses (e.g., from the channel sides and bed), the specific energy may be defined as follows:

$$E = y + \frac{aQ^2}{2gA^2} \tag{1}$$

where *y* is the water depth, *α* is the energy coefficient, *Q* is the volumetric flow rate, *A* is the cross-sectional area, and *g* is the gravitational acceleration. The energy coefficient can then be defined as follows:

$$\infty = \frac{\sum \mu^3 \Delta A}{l l^3 A} \tag{2}$$

where *A* is the total flow area and *u* is the velocity measured within an elemental area, Δ*A*. The flow regime (sub- or supercritical) and, thus, the Froude number of the flow govern the behaviour of the flow [14].

A parameter of specific interest is the critical depth of the channel in which the turbine is placed. When the water surface decreases to critical or subcritical depth, flow phenomena such as hydraulic jumps may form downstream to allow recovery to normal flow depth.

Simplification of crucial free-surface parameters in inland HK installations can be summarized by two fundamental effects:


#### *2.1. Free-Surface Effects of HK Turbines*

Free-surface effects are a critical aspect in riverine and tidal turbine array design, as the standing wave formed downstream of the turbine affects additional downstream turbines. Previous studies have concluded that the depth of the downstream water surface is strongly dependant on the *FrD*. Additionally, the blockage ratio also affects this free-surface change, albeit not as strongly as *FrD* [13].

Myers and Bahaj [12] found that when imposing the typical wake expansion on the water surface, this coincided with the increased elevation observed 4–5 *dt* downstream of the turbine rotor (as shown in Figure 2). In addition, due to the wake expansion coincident with the free surface, cumulative turbine placement at intervals smaller than the recovery length may cause the flow to approach critical depth, causing severe undulations in the water surface profiles (WSPs). Turbine operation and efficiency may also vary due to decreasing fluid velocity over the blades during operation. Accurate quantification of the WSPs around an array may be a challenge due to the multiple effects of turbulence, wake mixing, and superposition of WSP effects [12].

**Figure 2.** Wake expansion effect with free surface [12].

The presence of a support structure (tower/stanchion) also strongly affects the freesurface effect [12] (Figure 1). The same experiment also indicated the strong possibility of the formation of a hydraulic jump downstream of the turbine when flow is forced to a supercritical level due to the presence of the turbine and support structure.

Free-surface effects may be more pronounced for shallow turbines compared to turbines installed well below the free surface. It is also important to consider possible cross-sectional changes in the infrastructure where the turbine is placed, as this may alter/dampen/exaggerate these effects. Specified clearance coefficients have been investigated to limit the severity of decreased depths downstream of the turbine, or possible exposure of the turbine. Birjandi [15] proposed a clearance coefficient, *Ch*, defined as follows:

$$\mathbb{C}\_{h} = \frac{H}{L} \tag{3}$$

where *H* is the turbine submergence depth (i.e., the height between the top of the turbine and the water surface) and *L* is the rotor diameter. The recommended clearance coefficients for commercially available turbines can be seen in Table 1.


**Table 1.** Clearance coefficients for commercial HAHTs.

#### *2.2. Backwater Effect*

A turbine acts as a blockage in the channel and results in energy loss in inland flow infrastructure (where flow is constrained). When Fr << 1 and subcritical flow is prevalent, the backwater effect (damming upstream) may extend a large distance upstream as well as causing significant damming. This depends greatly on the blockage ratio, which is a function of the turbine swept area (*AT*), channel flow area (*Ao*) (*BR*(%) = *AT Ao* ), and additional constrictions [16], as well as the theoretical to actual efficiency [17].

In a study on a pilot HK installation in an irrigation canal, the backwater effect from the presence of the turbine extended up to 2.7 km upstream, due to the flat slope and subcritical flow present in the channel [1] (Figure 3). The clearance coefficient for this installation had not yet been defined. The specific turbine studied in that project contained grids upstream of the turbine, which trapped debris and caused a further increase in the backwater effect, to the point of channel overtopping. The blockage ratio for the installation was around 12.5%.

**Figure 3.** Backwater effect due to turbine blockages [1].

Additional to the blockage ratio, the Froude number (*Fr*) or Froude number based on turbine diameter (*FrD*) of the flow can influence the backwater effect. A previous study analysing this effect drew the following conclusions [13]:


Previous studies have attempted to quantity the effective blockage of an HK device. Some have addressed this through the relationship of power extracted to total power dissipated by the devices [18], analytical relationships [11], and even enhanced Manning n-values quantifying the energy loss as a friction loss [4]. However, a simple formula quantifying the effective blockage for different turbine types and operational conditions is yet to be determined, and provides the motivation for this paper.

#### *2.3. Backwater Calculations*

The extraction of energy resulting from the HK device may also be analytically incorporated through the use of the momentum equation [19], where the power extraction term is added as a shear stress component (added to the effective shear stress caused by bed friction). Assuming gradually varied steady-state flow, the conservation of mass and momentum can be used to adjust the standard open-channel flow equation [20] with the addition of a term for artificial energy extraction [19]:

$$\left(1 - \frac{Q^2}{h^3 b^2 g}\right) \frac{\partial h}{\partial x} = \frac{\partial h}{\partial x} \frac{Q^2}{g h^2 b^3} - \frac{1}{\rho g b h} P \,\tau\_{eff} \tag{4}$$

The effective shear (*τeff*) is defined as a combination of the bed shear (*τo*) and power extraction added as a shear term (*τadd*):

$$
\pi\_{eff} = \pi\_o + \pi\_{add} = \rho \frac{\mathcal{g}}{\mathcal{C}^2} \mathcal{U}^2 + \frac{P\_\chi R}{\mathcal{U}} \tag{5}
$$

where *P* is the wetted perimeter, *C* is the Chezy friction coefficient, and *Px* is the term added for power extraction, which may be more useful to express in terms of *PA* being the power extracted per unit area, as the flow passes through a plane where *Px* = *PA* <sup>Δ</sup>*<sup>x</sup>* (Δ*x* being the change in distance). Such an effect can be seen graphically in Figure 4, where a 10% energy extraction term has been added.

**Figure 4.** Influence of artificial energy extraction on speed and depth of flow [19].

A previous study [21] attempted to determine the backwater curve for instances when either the cross-section varied, the channel slope changed, or there was an obstacle in the channel and gradually varied flow was present. The model was based on the Bernoulli equation between two cross-sections. The energy loss between the cross-sections (related to distance) was termed the hydraulic loss, *I*, which can be calculated as follows:

$$I = \frac{n^2 \mathcal{U}^2}{R\_H^{\frac{4}{3}}} \tag{6}$$

where *n* is the Manning roughness (s/m1/3), *Rh* is the hydraulic radius of the channel (m), and *U* is the velocity of the water (m/s). The change in water levels (Δ*z*) between two sections can then be determined between two significant cross-sections (e.g., 0 and 1) and calculated as shown in Equation (7), where *α* is the Coriolis coefficient and *U0* and *U1* are the average velocities over distance Δ*L*:

$$
\Delta z = \Delta L \left(\frac{I\_0 + I\_1}{2}\right) + \frac{\alpha}{2g} \left(\mathcal{U}\_0^2 - \mathcal{U}\_1^2\right) \tag{7}
$$

Very few models have been developed to attempt to predict the backwater effect mostly 1D analytical models [4,11]. Most studies have focussed on tidal arrays and using free-surface effects to determine the optimal number and placement of turbines [22,23]. Within tidal applications, free-surface effects are only of concern for tip clearance; therefore, these models have limitations within the application in steady inland channels where spatial constraints are of primary concern.

In a study by Kartezhnikova and Ravens [4], an increased Manning roughness coefficient was used on the channel section representing the hydrokinetic device. The *n*-value used was a function of the actual channel *n*-value, slope, water depth, device efficiency, blockage ratio, and device deployment density. This method can then be used to determine the hydraulic impact, as well as the impact of various device configurations.

The head loss associated with the channel friction (*hLt*) (used in the energy conservation equation) can be written as shown in Equation (8), as a function of cross-sectional area (*A*0), channel hydraulic radius (*Rh*), discharge (*Qn*), and the length over which the loss is applied (*L*).

$$h\_{Lt} = \left(\frac{Q\_n}{A\_0 R\_h^{2/3}}\right)^2 L \tag{8}$$

Based on the assumption that the upstream and downstream velocity and pressure heads are equal, and assuming that the drag loss is negligible, the following equation for an enhanced bottom roughness (*nt*) can be derived as a function of the total power dissipated (*hp*), change in elevation (Δ*z*), and channel Manning roughness coefficient (*n*).

$$m\_l = n\left(1 - \frac{h\_p}{\Delta z}\right)^{-1/2} \tag{9}$$

Lalander and Leijon [11] investigated the use of numerical and analytical models to determine the effects on upstream water levels in a river. The analytical models are dependent on the channel blockage (of the HK device) determining energy loss from the energy capture, as well as the energy losses in the wake. The numerical models are based on the same theory—that energy is removed, causing a power loss; thus, energy capture is a component of the total friction in the channel. The total head loss can be determined as the sum of the friction loss (Δ*hf*) and head loss caused by the turbine (Δ*ht*), and *Pt* is the total power of the turbine (W). The formulation of the stress term *τ<sup>f</sup>* is shown in Equation (11), where *f* is equal to the Darcy–Weisbach coefficient (unitless) and *U* is the velocity of water (*m/s*).

$$
\Delta h\_{\text{tot}} = \Delta h\_f + \Delta h\_t = \frac{L}{\rho g R\_h} \ast \tau\_f + \frac{P\_t}{\rho g \mathcal{Q}} \tag{10}
$$

$$
\tau\_f = \frac{f\rho}{8} \ast \mathcal{U}^2 \tag{11}
$$

It is also important to consider the blockage effect, which can increase the turbine power output [22,24].

#### *2.4. Summary of Literature*

When HK devices are placed in array schemes in inland channels/rivers, the cumulative effect and inter-effect of these devices should be well understood to avoid unfavourable free-surface effects. As shown in the organogram (Figure 5), the blockage resulting from the HK device (in typical subcritical, flat-sloped channels) may have multiple subsequent effects influencing downstream installations (within the array) as well as upstream flow conditions. Neglecting the free-surface effects in high-blockage cases may result in exposed downstream turbines (i.e., freeboard reduced), hydraulic jump formation (enforced critical flow), and upstream damming effects, and may also lead to potential blade-tip cavitation problems. Accurately quantifying the water surface effects is a challenge, and guidelines to avoid unfavourable conditions may be extremely useful.

**Figure 5.** Water surface deformation from HK devices in inland flow infrastructure.

For inland systems, the backwater effect and the calculation thereof are of primary concern. Existing models include multiple unknowns and assumptions of values required as inputs that may not be available to users. Development of a clear, simple, effective blockage approximation of an HK device that can be applied to typical backwater calculations is necessary.

#### **3. Validation of CFD Models**

Comprehensive validation of a developed analytical model requires a dataset from a range of turbine styles under various operational conditions. Computational models offer an alternative to physical scale models to simulate the complex flow physics around HK devices. However, the importance of model validation through a physical model dataset should not be neglected. A number of studies have tested and validated modelling approaches for HK turbines [25,26].

Developments in tidal energy have led to multiple large-scale analyses of tidal turbine arrays, where simplified numerical models are used for array schemes [27–30]. With the ever-increasing availability and reliability of computational power, computational fluid dynamics (CFD) simulations are being used to model complex external effects and more accurately resolve fluid dynamic and wall effects [31]. Additionally, blade element momentum (BEM) theory is often used as a rotor modelling technique that also significantly reduces computational load [32].

CFD models may be used to resolve the effects of turbulence at the sub-grid-scale level, and are being used more often for first-order analysis or design. They have the potential to offer more comprehensive solutions and insights when their limitations are understood. Accurately representing turbulent flow in CFD is imperative due to its strong dependence on initial conditions, as well as the wide range of scales (eddies) present in the flow. Most often, statistical approaches based on the Reynolds-averaged Navier–Stokes (RANS) equations, with eddy viscosity models for turbulence closure, are used [33].

CFD models are valuable tools in flow-field analysis, especially at the sub-grid-scale level and for resolving turbulent length scales. When they are correctly applied and their limitations are understood, they offer useful insights and an alternative to laboratory testing. Correctly representing the turbulent flow, along with prescription of initial conditions, is important—especially due the wide range of scaled eddies present in the flow. Approaches based on the Reynolds-averaged Navier–Stokes (RANS) equations are most often used with various turbulence closure models incorporated, such as eddy viscosity models (k-Epsilon and SST k-ω).

A number of recent studies have used the BEM embedded in CFD method [8,34,35], which is also widely used in wind turbine array models, with good representation in terms of experimental results [36,37]. Various authors [38,39] have used this method to analyse the flow field of turbines arranged in arrays—specifically for tidal turbine optimization. A study analysing the accuracy of RANS approaches revealed good correlation with the experimental data found when using a Reynolds stress turbulence model (RSM) coupled with a BEM blade modelling technique (RSM-BEM) [40]. The study also highlights the importance of using an RSM rather than standard eddy viscosity models, due to the strong anisotropic flow in the wake.

#### *3.1. CFD Models*

Three different turbines (shown in Table 2) were modelled using CFD and validated with experimental results (a variation of free-surface, wake, and performance measurements). The primary validation case used was the U.S. Department of Energy's Reference Model 1 (RM1) dual-rotor axial flow turbine, which was modelled in the St. Anthony Falls Laboratory (SAFL) at the University of Minnesota [10]. Free-surface measurements allowed validation of the multiphase CFD model and free-surface deformation.

Additionally, two three-bladed turbine models with experimental results found in the literature [41] and tested on site [1] were modelled as specified in [40], for further measurements and validation of the backwater calculation methods. Where possible, BEM-VD (virtual disk) models were used to reduce the computational expense required. However, for the turbines which customized NACA profiles (Smart Hydropower turbine used in [1]) the full rotor geometry was modelled using a sliding mesh.

**Table 2.** Turbines modelled in CFD.


#### *3.2. RM1 Model Validation*

A 0.5 m diameter dual-rotor axial flow tidal turbine was investigated in a laboratory setup. The relevant details of the experiment can be seen in Table 3, with additional details available in [42]. The experimental results [42] provide high-resolution wake measurements of the near- and far-wake flow field surrounding the turbine from −5 to 10 *dt*. These were used previously to validate the CFD procedure for a single-phase analysis [40].


**Table 3.** RM1 laboratory setup details [42].

Although the velocity profiles and performance metrics were adequately modelled with both the single-phase and multiphase models, correct approximation of the backwater effect through the multiphase modelling required validation. For this, free-surface measurements were collected for the optimal operational case (*TSR* = 5.1, *U* = 1.05 m/s) of the RM1 experimental setup, at a resolution of 1 *dt* (diameters) in the streamwise direction, and 0.4 *dt* in the cross-stream direction. The measurement zone was −5 *dt* to 10 *dt* downstream. Elevation data were sampled at 50 Hz for 120 seconds at each location using a Massa ultrasonic range sensor, allowing for both time-averaged and fluctuating water surface elevation analysis and CFD validation.

Siemens STAR-CCM+ Commercial modelling software was used to simulate the turbines. The computational domain representing the RM1 laboratory model is shown in Figure 6. A wall-bounded model was used, extending from −14 *dt* upstream to 16 *dt* downstream of the axis of rotation. The specified inlet length allowed full flow development prior to reaching the turbine axis of rotation. The outlet length ensured that no effects from the downstream boundary condition affected the near-wake behaviour. Previous studies have found that around 15 *dt* is usually adequate for the outlet boundary length [43,44].

A velocity inlet and downstream pressure outlet were specified as boundary conditions. The laboratory test turbulence and velocity values that were measured experimentally were specified at the inlet. Full development of the boundary layer on all surfaces (i.e., boundary walls, blades, and stanchion) was ensured through the specific turbulence model wall treatment and mesh resolution for each test case (details seen in [40]).

A virtual disk (VD) rotor modelling technique was used, and a BEM model was employed over the VD. This VD-BEM modelling approach has demonstrated good accuracy in the past, also significantly reducing computational costs [8,34,35,45,46]. A BEM tip-loss correction was incorporated using the Prandtl tip-loss correction method [47].

A Reynolds stress model (RSM) was used to allow more accuracy in prediction of possible flow anisotropy. The RS linear pressure strain two-layer (RS-LPS2) model [48,49] was found to model the near wake accurately, with a low y<sup>+</sup> wall treatment on the turbine and turbine structure (y<sup>+</sup> < 1).

Only transient simulations were performed, where unsteady terms were discretized using a 2nd-order implicit scheme. A time step ensuring Courant numbers of less than 1 over the domain was ensured as far as possible (some cells exceeded 1 at the blade tips, where smaller time steps did not change results and, therefore, larger Courant numbers were allowed to reduce computational time). The time steps were around 0.003 s for the BEM-CFD models. However, the time steps varied over each approach, depending on the results of the grid convergence index (GCI).

**Figure 6.** Computational domain with grid refinements: (**A**) near wake, (**B**) blades, and (**C**) free surface.

The computational domain (Figure 6) consisted of a polyhedral mesh with grid refinements in the free surface, near wake, far wake, and surrounding the turbine structure. In the past, results indicated the importance of fine grids to track the tip vortices when using simplified RANS models [31]. For grid refinement, an adapted GCI method [50] was used, due to the variance in grid sensitivity in the different regions. Separate regions' mesh sizes were incrementally decreased, until no changes in turbine performance, free surface, or wake behaviour were observed. As depicted in Figure 6, the final mesh sizes were around 14 million cells.

The smallest cell sizes were specified at the near wake and free surface, where a 10 mm minimum size proved to be adequate. Refinement in the near-wake region is imperative for accurate development of the complex near-wake behaviour, where both separated and attached flow exist [51]. Gibson and Launder [52] investigated the pressure fluctuation effects of capturing the boundary layer surrounding the turbine, and noted the importance of an accurate capture of boundary layer formation. For the VD rotor modelling technique, a minimum of four cells over the blade thickness was ensured for all meshes, which is the recommended minimum when using the VD method. A mesh base size of 14 mm (2.8% *dt*) proved adequate in the far-wake region.

The channel walls around the turbine were modelled as non-slip walls using a high-y+ wall treatment to ensure that the effects of the wall boundary layer were included in the simulation. A two-layer formulation may be applied to the linear pressure strain model (LPS2 model). All analyses used the two-layer formulation.

Simulation of the air–water interface (multiphase flow) may be approached in various ways. A multiphase analysis ensures a robust approach, but demands a higher computational load; therefore, a symmetry boundary condition is often used in a single-phase model [40].

Due to the necessity of free-surface measurements in this analysis, the free surface was modelled using a volume-of-fluid approach. This approach uses a 2nd-order discretization to compute a clear interface between the air and water. A volume-fraction variable is used to specify the spatial distribution of each phase. Cells with multiple phases are treated as mixtures, and the method is highly dependent on adequate mesh resolution. A high-resolution mesh (Figure 6C) was ensured on the free surface, with the final cell size determined through a GCI test focussed on a free-surface profile analysis.

Validation of the modelled wake and performance can be seen in [40]. Due to the importance of the backwater approximation, and correct modelling of the backwater effect caused by the turbine blockage, the free-surface measurements from a multiphase model were compared in this paper. Comparison of the experimental and CFD results can be seen in Figure 7. The results correlated well, with the computed free-surface behaviour fitting well with the experimental tests. The CFD results were recorded over one rotation and plotted as an envelope. The backwater effect was predicted well, with a maximum of 12 mm damming occurring upstream.

**Figure 7.** Comparison of experimental and computational water surface profiles for the RM1 tests: (**a**) experiment and CFD water surface graphics; (**b**) lateral WSE comparison; (**c**) longitudinal centreline WSE comparison.

#### **4. Methods**

The primary objective of the development of a backwater model is to allow a usable model with only basic inputs required. Thus, a mathematical formulation based on existing models, basic hydrodynamic principles, and a set of recommended assumptions (where/if information is not yet available) was carried out. The workflow for the backwater model's development can be seen in Figure 8, and is explained in the following subsections.

**Figure 8.** Workflow of backwater model development.

#### *4.1. Assumptions and Exclusions*

It is important to state the limitations of the model. Operational condition boundaries were set to within the typical canal operating conditions, where deployment of HK devices would be considered. Scenarios outside of these boundaries were not considered, as it was assumed that the use of such a model would not be necessary outside of these conditions (e.g., at supercritical flow conditions) [1,14,53].

The following limitations on flow conditions were set:


#### *4.2. Mathematical Formulation*

Energy losses in a channel are categorized and included with various approximations. For channel roughness, the friction losses can be accounted for by the Manning equation. Additionally, sudden losses due to channel features such as piers, bends, and drop structures have also been defined/estimated empirically, and can be included as form losses [14]. These are typically defined as eddy losses (*he),* included as an energy loss:

$$h\_{\mathfrak{c}} = \mathbb{C}\_L \frac{V^2}{2g} \tag{12}$$

where *CL* is the loss coefficient predefined for typical losses in a channel. The drop in water level due to a particular loss can be quantified/included by applying either the momentum or energy equation over a channel section, and the upstream and downstream sections (in which the energy loss exists). Additionally, an empirical approach may be used, where experimental results are used to determine an empirical relationship, such as that done by Yarnell in 1934 for bridge piers [54].

#### 4.2.1. Approach 1: Momentum Approach

A possible approach often used to determine the effect of an object/structure on the free surface (backwater effect) is the momentum approach. Energy losses occur due to flow separation, vortex generation, friction, and turbulence—all associated with the changes in velocity due to the presence of the turbine. In some cases, the presence of the turbine may also result in the formation of a hydraulic jump on the water surface, resulting in additional energy losses. Applying the momentum approach avoids inclusion of these individual energy losses by considering the change in momentum between an upstream and downstream section. Additionally, if drag can be quantified, the momentum approach may be favourable.

The simple momentum formulation between two flow sections can be used with the resistance of the turbine represented as a drag coefficient, as shown in Figure 9 and Equation (13), where the change in momentum is quantified by the hydrostatic forces upstream (*F1*) and downstream (*F2*) of the device (water level change), as well as the friction from the channel bed and walls (*Ff*) and the force due to the turbine (*FD*). This can then be rewritten to Equation (14), in terms of the drag force (*FD*) due to the presence of the turbine.

$$F\_1 - F\_2 - F\_D - F\_f = \rho Q (\beta\_2 V\_2 - \beta\_1 V\_1) \tag{13}$$

$$F\_D = \frac{1}{2}\rho\text{g}\text{Bl}\_1^2 - \frac{1}{2}\rho\text{g}\text{Bl}\_2^2 - \rho\text{Q}(\beta\_2 V\_2 - \beta\_1 V\_1) \tag{14}$$

**Figure 9.** Momentum approach schematic (adapted form [55]).

The drag coefficient for the HK device can be rewritten, which allows the determination of *h*<sup>1</sup> (the upstream water level) through knowledge of the drag coefficient and downstream flow conditions. However, this requires knowledge of the downstream conditions, which is not always possible in feasibility studies in the design phase. An alternative conservative analysis would be to use the normal flow depth as the downstream value.

#### 4.2.2. Approach 2: Energy Approach

The energy approach is more often used to determine backwater effects in typical open-channel flow scenarios. This is also a common approach in bridge pier modelling [14] when modelling the backwater effect of arch bridges [56], bridge piers, and even irregular structures such as wood jams [57].

The energy approach can be seen in Equation (15) and Figure 10. All energy losses between a point upstream and downstream of an HK device or blockage are quantified as terms that contribute to either the friction losses (*hf*) or local losses (*hl*).

$$\frac{V\_1^2}{2g} + h\_1 + Z\_1 = \frac{V\_4^2}{2g} + h\_4 + Z\_4 + \sum h\_{f1-4} + \sum h\_{l1-4} \tag{15}$$

All terms are defined in terms of velocity *U*), water depth (*h*), and distance from a datum (Z). The total energy (*TEL*), energy grade line (*EGL*), and water level/hydraulic grade line (*HGL*) are shown in Figure 10.

**Figure 10.** Energy approach.

Approximations of the loss due to the presence of the turbine have been attempted, such as those mentioned in Sections 2 and 3. However, preliminary tests of these methods indicated inaccuracies over a range of varying turbines and operational scenarios.

A method of including the turbine loss in the energy equation includes quantifying it in terms of a loss coefficient that has been calibrated to the turbine type and operating conditions. The energy loss due to the presence of the turbine (*ht*) can be written as a function of a loss coefficient (*α*), the freestream velocity (*U*), and the blockage ratio of the turbine, as shown in Equation (16).

$$h\_t = \kappa \frac{\mathcal{U}^2}{2g} \times \frac{A\_t}{A} \tag{16}$$

Although this approach is most often used in the literature, the loss *(ht*) may also be quantified as a pressure drop, which is then directly converted to an energy loss as follows:

$$h\_t = \frac{\Delta P\_t}{\rho \text{g}} \times \frac{A\_t}{A} \tag{17}$$

The pressure change (Δ*Pt*) is measured in the computational models as the pressure drop across the turbine, and applied (with consideration to the blockage ratio) to the channel area. This approach was followed to allow the use of CFD models to approximate energy loss due to the presence of the turbine. Other than the losses due to the direct pressure drop over the turbine, there are additional losses in the near wake, due to the turbulent flow. Using the Δ*Pt* approximation (Equation (17)) allows measurement of the pressure drop over the turbine and near-wake area, thus including additional losses.

#### 4.2.3. Validation of Pressure Drop Measurement in CFD Results

CFD has previously been used to measure the backwater effects from blockages such as bridge piers [58], with computed and measured levels showing almost identical results. Multiple methods have also been analysed and validated for determining the backwater effects of common structures found in river channels [59–61]. To validate whether the approximation for *h*<sup>t</sup> shown in Equation (17) holds true, the Δ*Pt* was measured in the CFD model for the RM1 validation case. The subsequently calculated loss (*h*t) was then compared to the measured backwater effects in the laboratory tests (as well as multiphase CFD analysis). Inclusion of the support structure blockage was incorporated using the Yarnell approximation. The Yarnell approximation for a single circular bridge pier (similar to the support stanchion) was implemented:

$$\left[\frac{\Delta y}{y}\right]\_{empirical} = K \left(K + 5F\_r^{\ 2} - 0.6\right) \left(15a^4\right) Fr^2\tag{18}$$

where Δ*y* is the backwater generated by the pier, *y* is the undisturbed flow depth, *F*<sup>r</sup> is the downstream Froude number, and α is the ratio of the flow area obstructed by the pier to the total flow area downstream of the pier (also referred to as the blockage ratio). *K* is used as a coefficient reflecting the pier's shape. To ensure that the Yarnell approximation and pressure loss (Δ*Pt*) calculation work independently, the RM1 model free-surface deformation was measured with and without the stanchion structure (Figure 11), and the results were compared to the backwater calculation using only the pressure drop, as well as including the stanchion through the Yarnell approximation.

**Figure 11.** (**a**,**b**) Velocity and (**c**,**d**) surface water measurements graphics for the RM1 full model vs. the RM1 rotor and nacelle only.

The pressure drop due to the presence of the HK device was found to be maximal when measured over the size of the turbine-swept area from 1 D upstream to 1.5 D downstream (Figure 12). The calculated backwater (Equation (17)) was then compared to the measured backwater (Table 4), which was calculated from the average disk pressure drop, as shown in Figure 13.

**Figure 12.** Pressure measurements over the horizontal and vertical planes (at the turbine hub height centerline).

**Figure 13.** Pressure measurements over the disk and planes upstream and downstream of the RM1 turbine and retaining structure.



As shown in Table 4, the calculated analytical backwater values from Equation (17) and the computational reading for disk Δ*P* result in a very similar *h*t, as would be the result of the RM1 device in channel flow.

#### 4.2.4. Lambda Approximation

To allow a simple empirical model for the determination of the energy loss due to the turbine*, λ<sup>T</sup>* was selected as the energy loss coefficient used in the energy equation:

$$
\lambda\_{\rm fl} = \lambda\_T \times \frac{A\_{\rm fl}}{A} \frac{\mathcal{U}^2}{2g} \tag{19}
$$

where *h*<sup>t</sup> is included as a loss in the energy equation (Equation (15)), and *λ*<sup>T</sup> is calculated as a function of the thrust coefficient (*Ct*):

$$
\lambda\_T = \frac{\mathbb{C}\_t \times \mathcal{U}^2}{2\mathfrak{g}} \tag{20}
$$

where *Ct* is a value that can be obtained from the manufacturer, calculated, or assumed in the pre-feasibility stage. For HAHTs, these thrust coefficients (*Ct*) usually range from 0.52 to 0.89 [9,62–64]. According to the actuator disk theory, *Ct* may be written in terms of the induction factor *a* [65]. It is also known that ideally, according to the Betz limit, *a =* <sup>1</sup> 3 ; therefore, the ideal and highest attainable *Ct* would be 0.88. Theoretically, according to the BEM theory, this should result in the highest velocity deficit in the near wake and, therefore, the "worst case" scenario for the operational conditions. Realistically, the values lie at an

upper limit of *Ct* = 0.8. The thrust coefficient can be calculated directly if the thrust force (*T*), inlet velocity (*U*), and swept area (*A*) are known:

$$C\_t = \frac{T}{\left(\frac{1}{2}\right) \rho U^2 A} \tag{21}$$

To justify the use of the *λ<sup>T</sup>* approximation, the validated CFD models were analysed, the pressure drop/total thrust was measured, and the subsequent backwater effect was determined. The calculated *ht* (through Equation (19)) was then compared to the *ht* determined through the Δ*Pt* (Equation (17)) results, as validated in Section 4.2.3.

Two approximations for *λ<sup>T</sup>* were included (calculated and assumed *Ct*). The model should be usable with only basic knowledge of the turbine installation and operating parameters; therefore, simple available metrics could be used to obtain a conservative result. Acceptable correlation between the experimental and calculated values created confidence to proceed with the model and build a larger dataset to analyse the model's accuracy at a larger operational variance from optimal conditions.

A regression to the mean approach was used to accumulate the necessary dataset for the analysis of the aforementioned calculation procedure and assumptions. This was required to reduce computational costs, and due to the lack of available data on various input parameters. The dataset was created through results from three models of turbines typically used for inland installations (Table 2). These CFD models were validated through benchmark validation using experimental results obtained at optimal performance points. The models were then varied in five primary operational states, namely:


These primary variables have been previously investigated and shown to influence the turbine thrust imposed on the flow area, and may therefore influence the backwater effect. Constraints were set to the variation of these variables to ensure that the computational models remained within realistic scenarios, whilst allowing insight into the effects of changes. The primary objective of the model validation was to ensure that a relative level of accuracy was obtained and, more importantly, a conservative approach to predicting the possible backwater effects caused by such a device.

A measure of accuracy of the methods is indicated by the mean square absolute error. The absolute error (MAE), rather than relative error (RMSE), was used to place greater emphasis of the larger backwater values (at higher blockage ratios) rather than uniform predictions over the range of backwater predictions (*ht*). Additionally, as the sample size changed in the analysis, the strength of the sample size effect was minimized when comparing MAE. The variance was also included to give an indication of the test conditions with greater variability, and under which test conditions the model (and assumptions) performed best.

$$
\sigma = \sqrt{\frac{1}{N} \sum \left( \frac{h\_t}{\mathcal{Y}\_{\text{exp}}} - \frac{h\_t}{\mathcal{Y}\_{\text{calc}}} \right)^2} \tag{22}
$$

The MAE (*σ*) and variance values calculated for various scenarios are shown in Table 5. The variations between experimental and approximated *ht* for variations in blockage ratio (*BR*) velocity (*U*) and FrD are shown in Figure 14.



**Figure 14.** Effects of *U*, *BR*, and *FrD* on the determined *ht* values.

From the results in Table 5 and Figure 14, the following observations can be drawn:


was also higher than the 0.8 approximation, indicating that the turbine operates closer to the Betz limit and ideal induction factor (*a*), which could be further tested and calibrated. The *Ct* calculation performed better in this case, predicting *Ct* = 0.89. Therefore, utilizing this assumption may be favourable for avoiding errors—especially when turbines with higher operational tip speed ratios are used.

Based on the small dataset obtained from the three turbine models, a *Ct* value for each turbine can be determined empirically, which could be improved with a larger dataset.

The recommended model also performed significantly better than models found in the literature, as well as needing less input data and background knowledge on specific turbine performance. This indicates the usefulness of this approximation.

#### **5. Conclusions**

Quantifying the backwater resulting from the effective blockage caused by the operating turbine remains a challenge for the deployment of HK turbines in inland infrastructure, such as canal systems. This paper shows a simple analytical model to estimate the effective blockage and backwater effect from HK devices. The development of the model followed a similar approach to what has been previously used for bridge pier modelling and quantifying blockages from such structures.

Three variations of typical examples of commercially available HK devices were modelled using a recommended CFD approach, and conditions were varied to allow testing over a range of operational conditions. The water surface profile measurements from the Reference Model 1 (RM1) scaled model experiments allowed validation of the water surface deformation obtained when using a VOF approach coupled with an RS-LPS2 Reynolds stress closure model and BEM-VD blade modelling approach. This highlights the usefulness of CFD models in HK energy development.

The developed backwater model allows a conservative approach with various levels of certainty attainable, depending on the input parameters installed. Although the recommended procedure is an extremely simplified approximation, the results obtained were significantly closer to available dataset of backwater effects than the methods found in the literature. The ease of use also makes this method useful for engineers and developers when detailed numerical models are not feasible.

This model also allows cumulative estimation of backwater, with simple inclusion of blockage structures or multiple turbines using available approximations such as the Yarnell equation. Additionally, due to the nature of flow in canals (flat slopes), subcritical conditions govern, and backwater effects extend a large distance upstream. This allows simple cumulative inclusion of blockages without complex computational modelling, making this approach simple and relatively accurate (or at least favourably conservative).

This approach allows room for further development and determination of calibrated thrust coefficients for typical turbines or typical operational conditions (similar to what was previously done for pier shapes, etc.). Additionally, a similar approach may be investigated for cross-flow turbines where experimental results are available.

Although this paper focussed on backwater determination, other free-surface effects such as the water level drop over the turbine, or possible hydraulic jumps downstream—are important considerations for array designs. Recommendations for clearance coefficients and turbine spacing (due to wake recovery) should be carefully considered prior to deployment.

**Author Contributions:** Conceptualization, C.M.N. and M.V.D.; methodology, C.M.N., M.V.D. and L.S.; software, C.M.N. and L.S.; validation, C.M.N. and C.H.; formal analysis, C.M.N., C.H., M.V.D. and L.S.; investigation, C.M.N. and C.H.; resources, C.H.; data curation, C.M.N. and C.H.; writing original draft preparation, C.M.N.; writing—review and editing, C.H., M.V.D. and L.S.; visualization, C.M.N. and C.H.; supervision, M.V.D. and L.S.; project administration, C.M.N. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding. Experimental results used are open-access and can be found in [65].

**Data Availability Statement:** Not applicable.

**Acknowledgments:** This work was supported by the University of Pretoria. The computational capabilities were made possible due to academic hours allocated by the Center for High-Performance Computing (CHPC) South Africa. Siemens STAR-CCM+ Simcenter software and support were provided by Aerotherm Computational Dynamics (Pty) Ltd. Experimental results were obtained from the Sandia Laboratories repository. All contributors are thanked for their kind assistance and support.

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

#### **References**

