Get your ducts in a row Şireli is current and turbulent.

E. Mete Şireli

Who should read this paper? The oceans are powerful – the potential for generating power with the ocean is of great interest to activists, scientists, technologists, policy analysts, and a barrage of other professionals, not to mention energy consumers. Stakeholders and technology and project developers in hydrokinetic turbine projects will be particularly affected by this research, but anyone involved in marine current energy may want to navigate this work. Why is it important? Marine tidal currents are more predictable than wind and solar power. In general, marine energy appears a good contender for the future of energy production: there is power to be harnessed in the ocean’s waves, tides, temperature changes, and even salt content. Unlike traditional water turbines located on embankment dams, marine currents, both tidal and not, are a generally underdeveloped resource. Both offer appealing predictability and power. Introducing a duct to a traditional horizontal axis turbine is not trivial. The aim of this paper is to present a method for determining the optimal duct/ blade configuration, and maximizing the power derived from a turbine located in a free stream of water. Șireli’s work clarifies the theory behind ducted turbines, which he notes is currently lacking. What flows from his work is an optimal design for ducted blades on horizontal axis turbines. His method is based on the Blade Element Momentum theory, which describes the flows of fluids in relation to the blades of a turbine. Clean Current Power Systems is poised to commercialize the technology within the next year, as it is just completing the first year of testing in Canada, with locations in Manitoba, expanding to Nova Scotia and the Northwest Territories.

About the authors E. Mete Şireli completed a Master of Science degree in Naval Architecture at Istanbul Technical University, and a Master of Applied Science degree in Coastal and Ocean Engineering at University of British Columbia. He currently works for Clean Current Power Systems Incorporated as a Senior Hydrodynamic Design Engineer, and has authored several of Clean Current’s patents.

Copyright Journal of Ocean Technology 2014 68 The Journal of Ocean Technology, Vol. 9, No. 2, 2014

DESIGN IMPLICATIONS OF DUCTS ON HORIZONTAL AXIS TURBINES E. Mete Șireli Clean Current Power Systems Inc., Vancouver, BC, Canada ABSTRACT The Horizontal Axis Turbine (HAT) has been successfully used in wind kinetic energy harvesting for decades. The new technologies emerging in the tidal and river kinetic energy harvesting fields are also using HATs but they are also investigating the implementation of ducted HATs. This is partly due to the smaller sizes of turbines, which make duct implementation relatively easier. This paper investigates the implications of ducted designs for underwater current turbines using the Blade Element Momentum (BEM) Theory and proposes analytical relationships for the design of ducted turbine blades. To do so, axial momentum theory for the ducted turbines is coupled with the blade element theory. Formulas for optimum blade geometries are derived for ducted and unducted turbines. With the help of these formulas, two blades are designed to operate at a given tip speed ratio in ducted and unducted configurations. The ducted and unducted geometries are tested using a commercial Computational Fluid Dynamics (CFD) software package. The implications of ducts suggested by the axial momentum theory are compared to the CFD test results. Structural implications on the blades are also investigated using a commercial Finite Element Analysis (FEA) software package. In addition to all, a method to estimate the back pressure factor of a diffuser at the initial design stage is also proposed and tested using CFD. It is found that the effect of ducts on the blade design is not trivial. The performance increase suggested by the theory is achievable with a blade designed specifically to the effects of the duct of interest. As the theory suggests, at the design tip speed ratio, the thrust load on the blade remains the same for both ducted and unducted blade designs, but the ducted blade generates more power due to the duct’s area ratio and back pressure factor. This causes the blade to have more torque and therefore correspondingly more stress and cavitation. Compared to the unducted turbine, the overall thrust on the ducted turbine is also higher by the duct’s area ratio and back pressure factor just like the increase in power it generates. All of which indicate that, in a farm setting, fewer ducted turbines are needed to harvest the same amount of energy as unducted turbines. Fewer turbines result in fewer moving parts to do the same work. KEYWORDS Blade Element Momentum; Computational Fluid Dynamics; Wind; Tide; River; Turbine Copyright Journal of Ocean Technology 2014 The Journal of Ocean Technology, Vol. 9, No. 2, 2014 69

NOMENCLATURE a

Tb

= Axial induction factor

= Blade thrust (Tt = Tb for



unducted) [ N ]

Ab

= Area of the blade disk [

c

= Local chord length [ m ]

Tt

copt

= Optimum local chord length

Tb for unducted) [ N ]

Cd

m2

]

[m]

= Drag coefficient

v∞

= Total thrust on the system (Tt = = Undisturbed free stream velocity



Cdopt = Optimum drag coefficient

va

C1

= Lift coefficient



Clopt

= Optimum lift coefficient

vi

Cn

= Local normal force coefficient



[ m/s ]

= Axial velocity at the blade disk [ m/s ]

= Velocity of the fluid at station i = 0, 1, 2, 3, 4 [ m/s ]

Cnopt = Optimum local normal force

α

= Angle of attack [ degrees ]



αopt

= Optimum angle of attack

C P

coefficient

= Power coefficient



[ degrees ]

CPmax = Maximum theoretical power

β

= Diffuser area ratio



γ

= Back pressure factor

coefficient

CQ

= Blade torque coefficient

η

= Turbine efficiency

CT

= Total thrust coefficient (CT =

λ

= Tip speed ratio

CT-Blade for unducted)

λr

= Local speed ratio

CT-Blade = Blade thrust coefficient (CT =

ν

= Kinematic viscosity of the fluid

CT-Blade for unducted)



p∞

ξ

= Undisturbed free stream static

pi

= Static pressure at station i = 0,

Pb

pressure [ Pa ] 1, 2, 3, 4 [ Pa ]

= Blade power or harvested



power [ N ]

= Maximum diffuser flare angle

ξef

at the trailing edge [ degrees ]

= Estimated effective flow angle

ξl

[ m2/s ]

at the diffuser exit [ degrees ]

= Approximated local flow angle



at the diffuser exit [ degrees ]

Qb

= Blade torque [ Nm ]

ρ

= Density of the fluid [ kg/m3 ]

r

= Local radius [ m ]

φ

= Flow angle [ degrees ]

R

= Maximum radius [ m ]

φopt

= Optimum flow angle [ degrees ]

Rb

= Blade tip radius [ m ]

ω

= Angular speed [ rad/s ]

Copyright Journal of Ocean Technology 2014 70 The Journal of Ocean Technology, Vol. 9, No. 2, 2014

INTRODUCTION The need for energy sources is increasing and the need for clean sources even more so. There have been recent developments in the tidal and river energy extraction devices that are bringing this industry closer to commercialization. These devices usually take the form of a horizontal axis turbine. Some of the technologies also utilize a duct or a shroud, in the form of a diffuser or an augmenter, in search of enhancing the flow and increasing the energy extracted per device. In the past, some experimentalists believed that flow augmentation will increase the power output of a turbine by the cube of the velocity increase through the annulus; yet in practice, the power increase was only proportional to the increased mass flow through the turbine. Such beliefs still persist today [Chaudhari et al., 2013; Ehsan et al., 2012]. A look at the development of the BEM theory may shed some light on why such a misunderstanding is still widespread: The optimization of a free stream turbine is achieved by maximizing its power harvesting capability. Therefore, the first challenge is to identify whether a theoretical limit for maximizing the power output of a horizontal axis turbine exists. Using the axial momentum theory, Lanchester [1915] came very close to finding this theoretical limit for a free stream turbine. It was Betz [1920], however, who was able to find this limit definitively. In all of these attempts, an unducted device was the envisioned baseline. Later, Glauert [1935] introduced the blade element theory to the momentum theory and hence established the BEM theory for the unducted rotor.

At this point, since the BEM theory linked the power harvested by the blades to the free stream velocity by a cubic relationship as a result of mass and momentum equations, a mistaken belief started to emerge that by increasing the velocity through the blades via flow augmentation, the power can be increased by the cube of this velocity increase. This belief is false because it ignores the effects of such flow augmentation devices on the pressure field and focuses only on the velocity changes. Since the original BEM theory lacked the inclusion of the effects of a duct, various attempts have been made to understand the effects of ducts on the performance and design of a free stream turbine accurately. De Vries [1979] summarized the fluid dynamics effects of wind energy conversion devices, where he also touched on the subject of flow augmentation by the so called “wind concentrators.” In his report, a systematic approach towards understanding the effects of diffuser augmented wind turbine was presented. However, this approach was not developed further to a complete theory for ducted turbines. Using the principle of superposition in axial momentum theory, van Bussel [1999] managed to shed a significant light on the ducted turbine theory. He was able to show clearly that the expected power and overall drag increase from a ducted turbine may only be proportional to the increase in mass flow through the blades. He also showed that the turbine rotor had to be loaded to the same magnitude regardless of it being ducted or unducted. Later Lawn [2003] and Jamieson [2008] came to the same

Copyright Journal of Ocean Technology 2014 The Journal of Ocean Technology, Vol. 9, No. 2, 2014 71

conclusions using similar yet different models for ducted turbines. The goal of this paper is to use the simple BEM theory (without wake rotation) to shed further light on the effects of a duct on the performance of a horizontal axis turbine and, most importantly, on the design of the optimum blade geometry for maximum power generation. The simple BEM theory combines van Bussel’s axial momentum theory for ducted turbines and blade element theory together. Because axial momentum theory is used, the relationships derived are for the case without wake rotation. Since the rotational component in reality is very small, this approximation is simple yet accurate enough to understand the phenomenon happening in ducted turbines. The further intent of this paper is to show that the effect of ducts on the blade design is not trivial [Aranake et al., 2013], and that extensive CFD and FEA analyses [Luquet et al., 2013] required in the search for optimum ducted blade designs can be shortened significantly by the proposed method. AXIAL MOMENTUM THEORY FOR THE UNDUCTED HAT The axial momentum theory for the unducted HAT is well known. It is repeated here for completeness and as a foundation for the ducted HAT theory. The schematic of the stream tube of an unducted HAT during optimum operation condition is shown in Figure 1. We write the mass and momentum conservation for this stream tube [Lanchester, 1915; Betz, 1920; Glauert, 1935].

Figure 1: Stream tube of an unducted turbine during optimum operation (Blade disk: 1–2).

Mass conservation: ρ ∙ v0 ∙ A0 = ρ ∙ v1 ∙ A1 = ρ ∙ v2 ∙ A2 = ρ ∙ v4 ∙ A4 (1) Equation 2 and Equation 3 refer to the velocities and the areas for an infinitesimally thin disk plane, respectively. va = v1 = v2 (2) Ab = A1 = A2 (3) Hence, Equation 1 simplifies to Equation 4. v0 ∙ A0 = va ∙ Ab = v4 ∙ A4 (4) Momentum conservation: 2

2

p0 ∙ A0 + ρ ∙ A0 ∙ v = p4 ∙ A4 + ρ ∙ A4 ∙ v 4 0

(5)

p∞ = p0 = p4 (6) For A4>A0, there must be a net axial force = thrust Tb. This must be exerted on the blade disk. Using Equations 5 and 6: 2

2

A4 ∙ v ) (7) Tb = p∞ (A4 – A0) = ρ(A0 ∙ v – 0 4

Copyright Journal of Ocean Technology 2014 72 The Journal of Ocean Technology, Vol. 9, No. 2, 2014

Rewriting the thrust equation using the axial induction factor, we get:

Further, using Equation 4: Tb = ρ ∙ Ab ∙ va (v0 – v4 ) (8) The relationship between the velocities v0, va and v4 can be obtained by using the Bernoulli equation. Because the streamlines break at the blade plane, it should be written before and after the blade plane, separately. p0 + 0.5 ∙ ρ ∙

2 v = 0

p2 + 0.5 ∙ ρ ∙

2 v = 2

0.5 ∙ ρ ∙

2 v + 1

p1 (9)

0.5 ∙ ρ ∙

2 v + 4

p4 (10)

After summing Equations 9 and 10, substituting Equations 2 and 6, and some rearrangement, we find the following relationship: p1 – p2 = 0.5 ∙

2 2 ρ(v – v ) 0 4



(16)

Since the flow is assumed frictionless, all the power due to thrust is converted into rotational shaft power. Thus, Pb = Tb ∙ va (17) Pb = 4a(1 – a)2 ∙ 0.5 ∙ ρ ∙ Ab ∙v 03



(18)

By definition, the thrust and power coefficients, respectively CT and Cp are: 2 CT = Tb ⁄ (0.5 ∙ ρ ∙ Ab ∙ v ) 0



(19)



(20)

And

(11)

The pressure drop in Equation 11 is the reason for the thrust in the system. It happens at the blade plane and therefore: Tb = (p1 – p2) Ab = 0.5 ∙

Tb = 4a(1 – a) ∙ 0.5 ∙ ρ ∙ Ab ∙ v 02

2 2 ρ(v – v ) 0 4 Ab (12)

Using Equations 8 and 12, we find the following important relationship: va = 0.5(v0 + v4) (13) Here, we introduce the axial induction factor a, such that:

3 CP =Pb ⁄ (0.5 ∙ ρ ∙ Ab ∙ v ) 0

Substituting Equations 16 and 18: CT = 4a (1– a)

(21)

And CP = 4a(1 – a)2 (22) Here, we define the efficiency as the ratio of the CP to CT, namely: η = CP ⁄ CT (23)

va = (1– a) ∙ v0 (14)

η = (1– a)

And from Equations 13 and 14:

This efficiency definition is important from the turbine farm perspective because while CP corresponds to the harvested power, CT to the

v4 = (1–2a) ∙v0 (15)







(24)

Copyright Journal of Ocean Technology 2014 The Journal of Ocean Technology, Vol. 9, No. 2, 2014 73

extracted, as the extracted power is Tb multiplied by the free stream velocity v0. One needs to maximize the harvested power while keeping the overall thrust as low as possible to increase the farm’s overall efficiency. AXIAL MOMENTUM THEORY FOR THE DUCTED HAT: THE GENERAL CASE In the past, some researchers looked into the effects of diffusers on HATs systematically [De Vries, 1979]. Others tackled the problem using axial momentum theory, which is simple yet elegant [van Bussel, 1999; 2007]. The axial momentum theory for the ducted HATs provides a general case from which the unducted turbine solution can be obtained as a special case. The theory is further expanded here to include the back pressure estimation for the initial design, and later coupled with the blade element theory to be able to calculate blade geometries.

Figure 2: Stream tube of a duct (shaded) without a turbine (Inlet: 1–2).

va = β ∙ v3 (28) If the diffuser is designed correctly, the expansion of the stream tube will continue even after the flow has left the diffuser’s outlet. Therefore, the velocity at the diffuser outlet is related to the undisturbed free stream velocity by the back pressure factor γ, which is a function of the diffuser’s geometry. Thus: v3 = γ ∙ v0 (29) And, for an empty diffuser:

Duct without Turbine: Empty Diffuser The schematic of the stream tube of a diffuser without a turbine is shown in Figure 2. The conservation of mass will provide the relationships between the velocities inside and outside the duct. ρ ∙ v0 ∙ A0 = ρ ∙ v1 ∙ A1 = ρ ∙ v2 ∙ A2 = ρ ∙ v3∙ A3 (25) Using the notation in Equations 2 and 3: va ∙ Ab = v3 ∙ A3 (26) For β denoting the area ratio of the diffuser: β = A3 ⁄Ab (27)

va = βγ ∙ v0 (30) So, for increasing the velocity at the narrowing of the diffuser, one has to not only keep β>1, but also γ≥1, as demonstrated later. Estimating γ at the Initial Design Stage The presented method is for estimating the back pressure factor γ at the initial design stage. The designer should strive to achieve and surpass the γ value calculated by using this method. The current method should only be used for estimating the initial exit angle of the diffuser, which needs to be revised and improved though CFD analyses and iterations. Copyright Journal of Ocean Technology 2014

74 The Journal of Ocean Technology, Vol. 9, No. 2, 2014

We assume an axially symmetric duct with a sharp trailing edge to sustain the rear stagnation point at the trailing edge (Kutta condition). Due to the Kutta condition, the flow inside and outside the duct will be forced to have the same velocity at the trailing edge of the empty duct. If the flow inside the duct is assumed to be attached at all times, the flow should leave the duct at the same angle that the trailing edge makes with respect to the symmetry axis.

The flow outside the duct, which is parallel to the symmetry axis, is the undisturbed free stream flow with the velocity v0, which makes the angle ξef with the flow leaving the duct at the velocity v3. In an ideal fluid, far downstream of the duct (Figure 2), v4 = v0, which makes v3 a component of v0, such that: cos (ξef ) = v0 ⁄ v3 (34) Substituting Equation 29 in Equation 34:

We further assume that the flow is perfectly parallel to the symmetry axis at the symmetry axis and that it gradually and linearly conforms to the trailing edge angle as the flow fans out at the outlet (Figure 3), such that the local angle of the flow is given by Equation 31.

Figure 3: Schematic showing the change in the flow angle with radius at the duct outlet.

ξ1 = ξ ∙ r ⁄ R







(31)

The effective flow angle must therefore be the weighted average of the local flow angle: R

R

O

O

ξef = ∫ ξl ∙ 2πr ∙ dr ⁄ ∫ 2πr ∙ dr

γ = 1 ⁄ cos (ξef ) (35) As stated before, this relationship is for the initial estimate of γ. It assumes that the flow is inviscid and that there are no losses in the duct. In practice, net back pressure factors up to 1.2 are achievable with a good duct design – despite viscous losses – without increasing the overall thrust. The calculation of ξef can be further improved with the inclusion of a hub and the corresponding radius as the lower bound of the integrals, and by considering viscous losses at the boundaries. Duct with Turbine The schematic of the stream tube of a diffuser with a turbine is shown in Figure 4. Due to the existence of the turbine, the axial flow through the turbine is reduced by (1 – a) as in Equation 14, but due to the influence of the diffuser, the axial flow is increased by βγ as in Equation 30. Thus:

(32) va = βγ ∙ (1 – a) ∙ v0 (36)

With Equation 31, Equation 32 yields:

And, from Equation 28: ξef = ξ ∙ 2⁄3







(33) v3 = γ ∙ (1 – a) ∙ v0 (37)

Copyright Journal of Ocean Technology 2014 The Journal of Ocean Technology, Vol. 9, No. 2, 2014 75

Tb = (p1 – p2) Ab (43) 2

Tb = 4a (1 – a) ∙ 0.5 ∙ ρ ∙ Ab ∙ v 0

(44)

The total thrust Tt on the ducted turbine system can be obtained by the momentum conservation:

Figure 4: Stream tube of a duct (shaded) with a turbine (Blade plane: 1–2).

However, since there is not any radial expansion in the far wake, Equation 15 still remains valid. Using the Bernoulli equation on either side of the blade plane, one can find the pressure drop across the blade plane:

2 2 p0 ∙ A0 + ρ ∙ A0 ∙ v = p 4 ∙ A4 + ρ ∙ A4 ∙ v 0 4 (45)

After substituting Equation 40 and some rearrangement: 2 2 A4 ∙ v ) (46) Tt = p∞ (A4 – A0) = ρ (A0 ∙ v – 0 4

Mass conservation states:

2 2 0.5 ∙ ρ ∙ v + p1 (38) p0 + 0.5 ∙ ρ ∙ v = 0 1

v0 ∙ A0 = va ∙ Ab = v4 ∙A4 (47)

2 2 p2 + 0.5 ∙ ρ ∙ v = 0.5 ∙ ρ ∙ v + p4 (39) 2 4

Just as before:

Using Equation 15, which is still valid for the ducted turbine case, Equation 36 and Equation 47 in Equation 46, one finds:

p∞ = p0 = p4 (40)

Tt = βγ ∙ 4a (1– a) ∙ 0.5 ∙ ρ ∙ Ab ∙ v 02

va = v1 = v2 (41)

This implies that the thrust on the duct is due to the power harvested by the blades, as for the unducted turbine, β = 1, γ = 1 and Tt = Tb, and hence the ducted turbine equations are for the general case. Since the flow is ideal, all the power due to thrust is converted into rotational shaft power. Thus,

After summing Equations 38 and 39, substituting Equations 40 and 41, and some rearrangement, we find the following relationship: p1 – p2 = 0.5 ∙ ρ (v02 – v42 )

(48)

(42)

This pressure drop is identical to the one given in Equation 11. This means the pressure drop is independent of the existence of a duct. Therefore, the thrust on the blade plane is not affected by the presence of a diffuser. Hence:

Pb = Tb ∙ va (49) 3 Pb = βγ ∙ 4a(1– a)2 ∙ 0.5 ∙ρ ∙ Ab ∙ v 0

(50)

One can also obtain the same solution using the total thrust on the system and the axial Copyright Journal of Ocean Technology 2014

76 The Journal of Ocean Technology, Vol. 9, No. 2, 2014

flow velocity leaving the system, in which case: Pb = Tt ∙ v3 ⁄ γ







For which, a = 1/3 is found for maximizing Cp. This condition is true for both ducted and unducted cases. Thus, for a = 1/3:

(51) CPmax = βγ ∙ 16 ⁄ 27





(61)

Based on the blade plane area Ab, the blade thrust coefficient CT-Blade, the overall system thrust coefficient CT and the power coefficient CP are:

CT = βγ ∙ 8 ⁄ 9





(62)

CT-Blade = 8 ⁄ 9





(63)

CT-Blade = 4a (1– a)

(52)

η = 2 ⁄ 3





(64)

CT = βγ ∙ 4a (1– a)

(53)

The equations clearly show that from a farm perspective, the ducted turbine provides no additional gain in efficiency as far as harvested over extracted power is concerned. However, each turbine’s power harvesting capacity is increased by the factor of βγ, while maintaining the blade thrust the same. This means that the ducted turbines’ main shaft bearings carry the same axial load as unducted turbines of the same size with the added thrust from the duct being supported by stationary parts instead.





CP = βγ ∙ 4a (1 – a)2 (54) And: η = CP ⁄ CT (55) η = (1 – a)







(56)



THE MAXIMUM POWER CONDITION OPTIMUM BLADE GEOMETRY We use the general form of the power coefficient: CP = βγ ∙ 4a (1– a)2 (57) The axial induction factor that maximizes CP can be found by the following condition: dCP ⁄ da = 0 & d2 CP ⁄ da2 < 0

(58)

Optimum blade geometry equations derived here are for the general case without wake rotation. For unducted turbine equations, one has to set β = 1 and γ = 1. Optimum Flow Angle and Pitch For the case without wake rotation, the flow angle at the blade plane due to blade rotation is as follows (Figure 5):

Since βγ is independent of a: dCP ⁄ da = 4 (1– a)(1– 3a) = 0



d2 CP ⁄ da2 = 4(3a – 1) + 12(a – 1) < 0

(59) (60)

tan(φ) = va ⁄ (ωr) = βγ ∙ (1– a) ∙ v0 ⁄ (ωr) (65) For a = 1/3 and local speed ratio λr = ωr/v0, the optimum flow angle φopt can be obtained:

Copyright Journal of Ocean Technology 2014 The Journal of Ocean Technology, Vol. 9, No. 2, 2014 77

In the blade element theory, the thrust on the same circumferential differential element is given by Equation 70, where B denotes the number of blades.

Figure 5: Velocity triangle and flow angle φ at the blade plane without wake rotation.

tan (φopt ) = 2 ∙ βγ ⁄ (3λr) (66) Optimum section pitch θopt of the blade can be calculated using the optimum angle of attack αopt (Figure 6): θopt = φopt – αopt (67) Here, αopt corresponds to the maximum value of Cl/Cd. Optimum Chord Optimum chord can be calculated by equating the blade thrust equations of blade element and momentum theories. For a circumferential differential element: dAb = 2πr ∙ dr





(68)

From Equation 44: 2 dTb = 4a (1– a) ∙ ρ ∙ v ∙ 0 πr ∙ dr

(69)

dTb = 0.5 ∙ ρ ∙ Cn ∙ c ∙ B ∙ v2rel ∙ dr

(70)

Cn = Cl ∙ cos (φ) + Cd ∙ sin (φ)

(71)



From velocity triangle at the blade plane, one finds: sin (φ) = va ⁄ vrel = βγ ∙ (1– a) ∙ v0 ⁄ vrel

(72)

Equating Equation 69 and Equation 70, and substituting Equation 72, we get the local chord length: c = 8 ∙ a ∙ πr ∙ sin2 (φ) ⁄ [(βγ)2 ∙ (1– a) ∙ B ∙ Cn] (73) For a = 1/3, the optimum chord length copt can be obtained: copt= 4 ∙ πr ∙ sin2 (φopt) ⁄ [(βγ)2 ∙ B ∙ Cnopt] (74) Cnopt = Clopt ∙ cos (φopt) + Cdopt ∙ sin (φopt ) (75) Here, Clopt and Cdopt correspond to the maximum value of Cl/Cd.

APPLICATION OF THE THEORY

Figure 6: Schematic representation of section pitch θ, angle of attack α and flow angle φ including blade chord c and lift and drag coefficients C1 and Cd.

To test the claims of the axial momentum theory summarized inclusively in Equations 61 to 64, a series of CFD analyses were performed for underwater current turbines. For this purpose, a three-bladed geometry was envisaged. Two blades, one for unducted and another for ducted operation, Copyright Journal of Ocean Technology 2014

78 The Journal of Ocean Technology, Vol. 9, No. 2, 2014

were designed using the optimum blade geometry formulas previously derived using the BEM theory without wake rotation for the tip speed ratio of 4.2, which is the empirically recommended tip speed ratio for a three-bladed design. First, an unducted blade was designed and numerically tested. Later, a duct with an area ratio β of 1.5 was outfitted on the same blade to observe the performance change. And last, a ducted blade was designed for the same duct and tested. In all cases, the tips of the blades were left truncated and the hub/blade diameter ratio of 0.35 was kept for simplicity and consistency. The hub geometry was identical in all cases (Figure 7). The blade diameter was chosen to be 10 m and the CFD analyses were performed at 2.7 m/s free stream flow velocity. Foil Section NACA 4415 foil section was used throughout the blade profiles due to its publicly available data (Figure 8). Standard surface roughness properties of NACA 4415 were used to account for the fully turbulent boundary layer simulated in the CFD environment. The fluid dynamic properties of NACA 4415 for standard surface roughness at the Reynolds Number of 6,000,000 (Figure 9 and Figure 10) are readily available in the literature [Abbot and von Doenhoff, 1959]. Using the information in Figure 9, the optimum angle of attack αopt, at which Cl/Cd reaches its maximum value of 64.541, was calculated to be 5.82°. With this optimum angle of attack, the Cd value of 0.015 was found from Figure 10.

Figure 7: Hub geometry with main dimensions.

Figure 8: NACA 4415 foil profile normalized about chord.

Figure 9: C1/Cd vs angle of attack α for NACA 4415 with standard roughness at Re: 6,000,000 [Abbot and von Doenhoff, 1959].

Figure 10: Cd vs angle of attack α for NACA 4415 with standard roughness at Re: 6,000,000 [Abbot and von Doenhoff, 1959].

Copyright Journal of Ocean Technology 2014 The Journal of Ocean Technology, Vol. 9, No. 2, 2014 79

Table 1: Unducted blade chord and section pitch distribution with Reynolds Numbers for v∞ = 2.7 m/s, a = 1/3, λ = 4.2, ν = 1.004∙10-6 m2/s, ρ = 998.5 kg/m3, βγ = 1.

Table 2: Ducted blade chord and section pitch distribution with Reynolds Numbers for v∞ = 2.7 m/s, a = 1/3, λ = 4.2, ν = 1.004∙10-6 m2/s, ρ = 998.5 kg/m3, βγ = 1.533.

The 10 m blade diameter and 2.7 m/s free stream velocity were purposefully chosen to keep the Reynolds Numbers in all cases about 6,000,000 along the blade span to agree with the published test data (see Table 1 and Table 2 for details). The Empty Duct Determining the Geometry To be able to design the ducted blade geometry, one has to determine the undisturbed axial velocities at the intended blade plane, and hence in the empty duct’s annulus.

Equation 30 gives the relationship between the undisturbed free stream flow velocity and the undisturbed axial velocity through the duct’s annulus. As mentioned earlier, a duct in the form of a diffuser with an area ratio β of 1.5 is envisaged for this exercise. The diffuser’s flare is a perfect cone for simplicity making an 18° angle with the symmetry axis to agree with the Falkner-Skan solution for maximum wedge angle for attached flow to minimize the duct’s length. Since the flare is conical, the trailing Copyright Journal of Ocean Technology 2014

80 The Journal of Ocean Technology, Vol. 9, No. 2, 2014

edge angle ξ is 18°, as well. Therefore, from Equation 33, ξef is found to be 12°. Thus, from Equation 35, the back pressure factor γ is estimated to be 1.022. Therefore, the overall duct factor βγ is estimated to be 1.533. During operation, the axial induction factor needs to be 1/3 for optimum operation as shown before. Using the duct factor βγ of 1.533 and a = 1/3 in Equation 36, it is found that the axial velocity through the annulus during operation would be 1.022 times larger than the undisturbed free stream. In this particular case, the stream tube will resemble to that shown in Figure 11. Therefore, the inlet section of the diffuser can be kept a perfect cylinder to minimize the drag of the diffuser during operation (Figure 12). CFD Tests of the Empty Duct Geometry The duct geometry was tested using CFD to determine the exact back pressure factor γ to be able to design the ducted blade accurately. The duct was subjected to 2.7 m/s free stream velocity in a medium of fresh water (ν = 1.004∙10-6 m2/s, ρ = 998.5 kg/m3) and the velocities at the intended blade plane were measured (Figure 13). The duct factor βγ was calculated based on Equation 30. The CFD calculated area averaged value of the duct factor βγ was found to be 1.533 (Figure 14). This value is identical to the estimated duct factor indicating that the calculated back pressure factor γ was, for this particular case, accurately estimated using the proposed method. One has to be cautious, however, in solely relying on this method to calculate the back pressure factor γ for ducts with complex curvatures.

Figure 11: Stream tube of a duct with turbine in operation (Duct factor βγ = 1.533).

Figure 12: Diffuser geometry with rounded main dimensions together with the hub.

Figure 13: Axial velocity contours at the axial symmetry plane at free stream velocity of 2.7 m/s.

Copyright Journal of Ocean Technology 2014 The Journal of Ocean Technology, Vol. 9, No. 2, 2014 81

Figure 14: Duct factor βγ distribution along the radius with area averaged value.

Figure 15: Unducted blade (left) and ducted blade (right) geometries.

As a case in point, Figure 13 and Figure 14 show strong non-uniformity in the flow close to the duct, which may accentuate vibrations and noise. In an ideal design, this aspect of the diffuser needs to be addressed and improved. As a result, the back pressure factor γ may be significantly improved from the initial estimation as long as the exit angles and diffuser curvatures are chosen such that the flow separation is minimized. Such an optimization of the duct is beyond the scope of this paper, and the conical diffuser configuration will be continued to be used throughout the paper.

Ducted vs Unducted Blade Geometry Equations 66 and 67 were used to calculate the optimum flow angle and section pitch distribution of the blades. Optimum values of NACA 4415 were used. Equations 74 and 75 were used to calculate the optimum chord distribution of the blades. For the unducted blade, the duct factor βγ = 1 was used. For the ducted blade, the duct factor βγ = 1.533 was utilized. In both cases, the number of blades B was three and the tip speed ratio λ was 4.2. The local tip speed ratios λr were calculated accordingly. For both blades, a uniform velocity profile was used to calculate blade geometries. The calculated chord and section pitch distributions are tabulated for the unducted and ducted blade in Table 1 and Table 2, respectively. The final blade geometries are shown in Figure 15 for comparison. As can be seen from the tables and the figures, both blades have similar average chord lengths (Unducted: 0.892 m vs Ducted: 0.836 m). The ducted blade, however, has a less pronounced taper but a more pronounced twist and pitch distribution. CFD ANALYSES Analysis Setup and Convergence A rectangular domain of H: 30 m, W: 50 m and L: 140 m was created to simulate the flow about the turbine (Figure 16). The turbine was situated approximately 30% of the total domain length away from the intended inlet. The domain was subdivided into a cylindrical near filed region of 35 m in diameter to account for the near field flow structures immediately surrounding the turbine. The near Copyright Journal of Ocean Technology 2014

82 The Journal of Ocean Technology, Vol. 9, No. 2, 2014

field domain was further subdivided into the rotational domain encompassing the blades of the turbine. The rotational domain is cylindrical in shape extending from the hub all the way to the diffuser. It is 1.7 m in axial length extending from 0.8 m upstream to 0.9 m downstream of the blade plane. In the far field region, an unstructured hex dominant mesh was preferred to reduce the number of elements (Figure 17). In the near field region, an unstructured tetrahedral mesh was used due to its ease of conformation to complex shapes and surfaces (Figure 18). By this means, the mesh density was also increased in the near field region without the need of specifically changing the element size in the region (Figure 19). In the rotational domain, the minimum necessary mesh density was used to maintain geometric properties of the NACA 4415 blade profile (Figure 20). Prismatic sublayers on the blade surfaces expanding into the rotational domain were further used to calculate the boundary layer formation about the blades accurately. A minimum of 15 layers were intended to be within the boundary layer on the blade surface (Figure 21). Far field, near field and rotational domain regions had approximately 200,000, 300,000 and 3,900,000 elements, respectively making the total count approximately 4,400,000 elements. Opening boundary conditions were used for the walls of the domain to simulate unbounded free stream flow behaviour (Figure 22). Properties of fresh water (ν = 1.004∙10-6 m2/s, ρ = 998.5 kg/m3) were used with Shear Stress Transport (SST) turbulence model in all Copyright Journal of Ocean Technology 2014

Figure 16: Overall dimensions of the computation domain.

Figure 17: Far field (Hex), near field (Tet) and rotational (Tet and prism) domain mesh.

Figure 18: Surface mesh definition on the turbine.

Figure 19: Near field (Tet) and rotational domain (Tet and prism) mesh close-up.

The Journal of Ocean Technology, Vol. 9, No. 2, 2014 83

Figure 20: Surface mesh definition on the blade tip.

Figure 21: Mesh about the blade tip (top) and root (bottom).

Figure 22: Boundary conditions of the overall domain (single arrows: velocity, double arrows: opening). 84 The Journal of Ocean Technology, Vol. 9, No. 2, 2014

analyses. The analyses were first run for 1,000 iterations in the steady state mode. For the rotational analyses (λ = 4.2), the results of the steady state analyses were used as initial conditions for the subsequent transient analyses. The transient analyses were setup to be run for 3,000 time steps. Each time step corresponded to a 6° incremental rotation of the blades. Both rotational and stationary analyses were run until a sufficient convergence was reached for each monitored physical property such as blade torque, thrust, etc. (Figure 23). Results Qualitative Comparison In Figure 24, the wake formations behind the three turbines are compared. The pressure field at the symmetry plane shows the pressure distribution caused by the turbines in operation. The spiral shaped pressure surface leaving the blades indicates the rotational wake for each case as calculated by CFD. In all three cases, the pressure surface tracks -600 Pa, which roughly corresponds to the pressure drop caused by the tip vortices shed after the blades. While ducted and unducted blades show similar wake structures immediately after the turbine, the unducted blade outfitted with the duct shows a very strong wake formation with a large pitch due to suboptimal operation condition. As expected, the ducted turbine shows smoother pressure gradient in the far wake than the unducted turbine due to reduced tip vortices and added diffusion. The smoother flow field of the wake of the ducted turbine is also visible in the streamlines shown in Figure 25 when compared to the unducted turbine and its duct-outfitted version. Copyright Journal of Ocean Technology 2014

Figure 23: Convergence criteria vs accumulated time step.

Figure 24: Wake formation behind turbines. From top to bottom: Unducted, Unducted+Duct, Ducted.

Figure 25: Streamlines through turbines. From left to right: Unducted, Unducted+Duct, Ducted.

Theory vs CFD The theoretical values of the hydrodynamic coefficients CPmax, CT-Blade, CT and η are calculated and tabulated in Table 3 (see Equations 61 to 64 for details). The column labeled Ducted/Unducted shows the ratio of the calculated values. As the theory suggests, the power harvesting capability of the ducted turbine as well as its overall drag are expected to increase by the duct factor βγ = 1.533 for the particular case at hand, while the blade thrust and overall efficiency remain the same. The theory considers an ideal fluid with no losses with a blade plane making use of the entire momentum disk having no hub. Because of these assumptions, the CFD based hydrodynamic coefficients were calculated either by using the blade swept area (radius spanning from the blade root to the blade tip) or the disk area (radius spanning from the centre of rotation to the blade tip), and for some other cases either by including or excluding the drag of the hub. No method of calculation is deemed superior here, yet provided for clarity, as one may arguably find one or the other more compliant with the theory. A definition list for the calculated hydrodynamic coefficients is given in Table 4. CFD calculated values of these hydrodynamic coefficients are provided in Table 5. To show the consequences of not matching the blade design to the duct’s effects, the Unducted+Duct case is included in Table 5. The Unducted+Duct case clearly shows suboptimum performance characteristics when its efficiency values are compared to that of the ducted and unducted turbines regardless of the hub drag component in the calculation Copyright Journal of Ocean Technology 2014 The Journal of Ocean Technology, Vol. 9, No. 2, 2014 85

Table 3: Theoretical values of hydrodynamic coefficients during optimum operation.

Table 4: Definition list for the calculated hydrodynamic coefficients.

Table 5: CFD calculated values of hydrodynamic coefficients during operation.

(η-1 and η-2). For this particular case, one could argue that the blade’s lower pitch will allow it to reach its maximum performance at a higher tip speed ratio due to higher velocities at the blade plane incurred by the duct. However, this would be the same mistake as looking at the velocity field without considering what is happening in the pressure field. To look at it correctly, one has to investigate the thrust coefficients of the blade. As the theory suggests, in a free flow

environment, whether the blade is ducted or not, the turbine will reach its peak performance point when the blade thrust coefficient reaches the value of 8/9. This, in fact, is the universal limiting value of reaching the peak performance point, more so than the Betz’ limit [Jamieson, 2008]. The blade thrust coefficients of the Unducted+Duct case indicate that the blade is already loaded close to its optimum value of 8/9. This means that the blade has already reached its peak Copyright Journal of Ocean Technology 2014

86 The Journal of Ocean Technology, Vol. 9, No. 2, 2014

performance point. In a free flow environment, the blade cannot be more heavily loaded without at its peak performance point without a corresponding drop in performance. Therefore, this configuration cannot perform any better at higher tip speeds as the blade will further drift away from its optimum loading condition. In fact, based on the CT-Blade-1 value of 0.938, one may argue that the peak performance point of this configuration resides at a lower tip speed, at which the blade is loaded closer to its optimum value of 8/9. When the CFD calculated Ducted/Unducted values are compared to the theoretically estimated values, a clear agreement becomes visible. The power coefficients (CP-1 and CP-2) show that with the blade correctly matched to the duct at the design tip speed, the performance is increased 1.516 times – close to the theoretical estimate of 1.533 times. Blade thrust coefficients (CT-Blade-1 and CT-Blade-2), however, indicate that the ducted blade is loaded slightly less with 0.974 times that of the unducted version, but still very close to the target of maintaining the blade load identical to that of the unducted blade as the theory suggested. The ducted turbine’s total thrust coefficients (CT-1 and CT-2) are also higher by the factor of 1.504, a value close to the theoretical estimate of 1.533, when the hub drag is included in the calculation. When the hub drag is not included (CT-3), the increase in the total thrust becomes 1.553. The blade efficiency of the ducted turbine is virtually the same as the unducted turbine: The Ducted/Unducted efficiency ratios are 1.008 (η-1) and 0.976 (η-2) with and without the hub drag, respectively. Overall,

the results are in good agreement with the theory in terms of estimating the effect of the duct on the hydrodynamic properties of an optimized turbine. Cavitation Cavitation percentages of the blades were calculated by dividing the cavitation bubble area (Figure 26) to half of the blade wet surface area. To calculate the cavitation percentage at different water speeds for the same tip speed ratio, the cavitation number (Equation 76) for that particular water speed was first calculated. Using this cavitation number and the water speed at which the analysis was run (2.7 m/s), the vapour pressure corresponding to the water speed of interest was then back calculated. σ = (Patm + ρgh – Pv) ⁄ (0.5 ∙ ρ ∙ v2) (76) The ducted blade clearly exhibits an increased amount of cavitation compared to the unducted blade (Figure 27). The unducted blade outfitted with the duct appears to be slightly worse than the ducted blade. Throughout the examined cavitation percentage range, however, the ducted blade shows an earlier inception at about 90% of the water speed at which the unducted blade reaches the same cavitation level. Therefore, it can be said, that in this particular case, the ducted blade has a 10% earlier cavitation inception of any cavitation percentage for a power increase of about 51.6%. Blade Stress (FEA) Blade stresses for all three blades were calculated using FEA. The pressures from the CFD calculations were mapped on the blade

Copyright Journal of Ocean Technology 2014 The Journal of Ocean Technology, Vol. 9, No. 2, 2014 87

Figure 26: Cavitation bubble on the blades at 4.0 m/s. From left to right: Unducted, Unducted+Duct, Ducted.

surfaces. A fixed support boundary condition was applied to the root of the blade where it interfaced with the hub. CFD calculations and theory suggest that the axial loading on the blades remains about the same when ducted and unducted blades are compared. However, due to a 51.6% increase in the harvested power at the common tip speed ratio λ of 4.2, the torque on the ducted blade is also higher than that of the unducted blade by the same amount. Due to this increase in torque and geometric differences in the blades a maximum stress of 55.6 MPa was calculated for the ducted blade. This is a 25.6% increase over the 44.4 MPa of maximum stress on the unducted blade for a power increase of 51.6%. Suboptimal unducted blade outfitted with the duct exhibits 53.6 MPa at the same tip speed ratio (Figure 28). Starting Torque The starting torque of the three blades was also investigated. For this comparison, the blade torque coefficient was used (Equation 77). CQ = Qb ⁄ (0.5 ∙ ρ ∙Ab ∙ v2 ∙ Rb)

(77)

Figure 27: Change in cavitation percentage with water speed.

Figure 28: Blade stresses during operation. From left to right: Unducted, Unducted+Duct, Ducted.

Table 6 shows the hydrodynamic coefficients during starting (λ = 0). As indicated by the column Ducted/Unducted, the torque is increased 76.6%, while the blade thrust only 19%. The higher increase in starting torque compared to the operation performance increase is due to the contribution of the duct’s flow augmentation as well as the higher pitch and twist of the ducted blade. The overall drag is increased 38.1%. However, the higher starting torque will allow the ducted blade to engage power extraction Copyright Journal of Ocean Technology 2014

88 The Journal of Ocean Technology, Vol. 9, No. 2, 2014

Table 6: CFD calculated hydrodynamic coefficients during starting.

Table 7: Definition list for the torque coefficients.

at an earlier water speed. Table 7 has the definition list for the torque coefficients used. The definitions of the remaining hydrodynamic coefficients are shown in Table 4. CONCLUSION AND FUTURE WORK The theories of ducted and unducted turbines were summarized. Optimum blade geometry formulas were derived for the case without wake rotation. The applied theory was tested using CFD techniques. The results show that the predictions of the theory are viable even for an unoptimized duct (no shape optimization). The results confirm that the power increase in a ducted turbine is proportional to the mass flow increase through the turbine as the theory suggests. The results also indicate that ducted turbine design results in power increase without any increase in blade thrust, while the harvested/ extracted power indicated by the efficiency η can be maintained at unducted turbine levels. All of which support the conclusion that, in a farm setting, fewer ducted turbines are needed to harvest the same amount of energy

as unducted turbines, which translates into fewer moving parts to do the same work. To further substantiate this theory, additional analyses of the ducted and unducted blades are required at tip speed ratios outside the design point. The actual locations of the peak power points need thus to be determined for each turbine. Since the present application of the theory is based on the case without wake rotation, it is likely that the peak power points will occur at slightly different tip speed ratios than the design tip speed ratio as the CFD includes the rotating wake effects. Also, the NACA 4415’s hydrodynamic properties need to be evaluated in the CFD environment, which may also contribute to any possible shift in the optimum tip speed ratio. Finally, the theory needs to be expanded to include the wake rotation. ACKNOWLEDGMENT I thank Mr. Russell Stothers and Dr. Mihai Platon for their editing and proofreading of the paper.

Copyright Journal of Ocean Technology 2014 The Journal of Ocean Technology, Vol. 9, No. 2, 2014 89

REFERENCES Abbot, I.A. and Von Doenhoff, A.E. [1959]. Theory of wing sections: including a summary of airfoil data. Dover Books on Aeronautical Engineering, ISBN-10: 0486605868, ISBN-13: 978-0486605869. Aranake, A.C.; Lakshminarayan, V.K.; and Duraisamy, K. [2013]. Computational analysis of shrouded wind turbine configurations. 51st AIAA Aerospace Sciences Meeting including the New Horizons Forum and Aerospace Exposition, January 7-10, Grapevine, Texas. Betz, A. [1920]. Das maximum der theoretisch möglichen Ausnutzung des windes durch windmotoren. Zeitschrift für das gesamte Turbinenwesen, Vol. 26, pp. 307-309. Chaudhari, C.D.; Waghmare, S.A.; and Kotwal, A.P. [2013]. Numerical analysis of venturi ducted horizontal axis wind turbine for efficient power generation. International Journal of Mechanical Engineering and Computer Applications, Vol. 1, Iss. 5, ISSN 2320-6349. De Vries, O. [1979]. Fluid dynamic aspects of wind energy conversion. AGARDograph No. 243, AGARD-AG-243. Ehsan, M.M.; Ovy, E.G.; Chowdhury, H.A.; and Ferdous, S.M. [2012]. A proposal of implementation of ducted wind turbine integrated with solar system for reliable power generation in Bangladesh. International Journal of Renewable Energy Research, Vol. 2, No. 3, pp. 397-403. Glauert, H. [1935]. Airplane propellers, In: Durand, W.F. (ed.) Aerodynamic theory, pp. 169-360, Springer, Berlin. Jamieson, P. [2008]. Beating Betz – energy extraction limits in a uniform flow field.

European Wind Energy Conference, Brussels, Belgium. Lanchester, F.W. [1915]. A contribution to the theory of propulsion and the screw propeller. Transactions of the Institution of Naval Architects, Vol. 57, pp. 98-116. Lawn, C.J. [2003]. Optimization of the power output from ducted turbines. Proceedings of the Institution of Mechanical Engineers, Part A: Journal of Power and Energy, Vol. 217, No. 1, pp. 107-117. Luquet, R.; Bellevre, D.; Fréchou, D.; Perdon, P.; and Guinard, P. [2013]. Design and model testing of an optimized ducted marine current turbine. International Journal of Marine Energy, Vol. 2, pp. 61-80. Van Bussel, G.J.W. [1999]. An assessment of the performance of diffuser augmented wind turbines (DAWTs). 3rd ASME/JSME Joint Fluids Engineering Conference, San Francisco, California. Van Bussel, G.J.W. [2007]. The science of making more torque from wind: diffuser experiments and theory revisited. Journal of Physics: Conference Series, Vol. 75, doi:10.1088/1742–6596/75/1/012010.

Copyright Journal of Ocean Technology 2014 90 The Journal of Ocean Technology, Vol. 9, No. 2, 2014