Computational Design Analysis of a Hydrokinetic Horizontal Parallel Stream Direct Drive Counter-Rotating Darrieus Turbine System: A Phase One Design Analysis Study

Abstract

This paper introduces a phase one computational design analysis study of a hydrokinetic horizontal parallel stream direct-drive (no gear box) counter-rotating Darrieus turbine system. This system consists of two Darrieus rotors that are arranged in parallel and horizontal to the water stream and operate in counter-rotation due to the incoming flow. One of the rotors directly drives an armature coil rotor and the other one a permanent magnet generator. A two-dimensional (2-D) and three-dimensional (3-D) computational fluid dynamic (CFD) simulation study was conducted to assess the hydrokinetic performance of the design. From a high computational cost and time perspective, the simulation setup was reduced from a 3-D to a 2-D analysis. Although useful information was obtained from the 3-D simulations, the output performance could be assessed with the 2-D simulations without compromising the integrity of the turbine output results. A scaled experimental design prototype was developed for static (non-movement of the rotors with dynamic fluid flow) particle image velocimetry (PIV) studies. The PIV studies were used as a benchmark for validating and verifying the CFD simulations. This paper outlines the prototype development, PIV experimental setup and results, computational simulation setup and results, as well as recommendations for future work that could potentially improve overall performance of the proposed design.

1. Introduction

Presently, there is an ever-growing interest in renewable energy extraction from biomass, ocean, wind, geothermal, and the sun (solar) due to the increase in the cost of energy produced from fossil fuels and the resulting detrimental effects on the environment [1,2,3,4,5,6,7]. These sources of sustainable energy have immense potential for providing the energy required to power electric grids and water desalination plants [8]. Of these source options, lakes, rivers, and oceans provide a renewable source of energy production with the advantage of being predictable many days in advance, reliable in terms of being stable during day and night, and substantially superior in energy density compared to wind and solar energies [9]. This source of clean and renewable energy can be used either directly for electricity generation [10] and/or for seawater desalination [11]. Renewable ocean energy can be divided into two main forms: wave and tidal energy [12]. Tidal energy takes advantage of the natural tide and current of seawater caused by gravitational fields [13,14,15]. Tidal energy convertors can be utilized for converting the kinetic energy of water to mechanical energy of shafts, similar to wind turbines [16,17,18,19,20].

Several turbine systems have been developed that can use hydrodynamic lift effectively in the efforts of extracting energy from a stream of water or wind [21,22,23,24,25]. Of these types, Darrieus turbines have demonstrated to be effective for both wind and waterpower [26,27,28,29]. The key benefits of Darrieus turbine system designs are: (1) the design of the system components and overall system is relatively easy as compared to other turbine types, (2) the system components can be adjusted and tuned at ground level, and (3) contrary to wind turbines, Darrieus turbines do not require a yaw system [30]. Several computational and experimental design evaluation attempts have been made, specifically in the United Kingdom (U.K.), the United States (U.S.), Japan, and Germany, to develop efficient hydrokinetic Darrieus turbines for small and large scale energy extraction [31].

From these works, it is widely accepted that parametric analysis of Darrieus turbines is necessary to optimize the design in the efforts of reducing the cost of manufacturing, maintenance, and maximizing performance. In comparing other turbine types, an apparent weakness of tidal turbines is the price of manufacturing. These turbines are required to be submerged in seawater, which is quite costly and there have been attempts by engineers to reduce the price of the turbine [32]. When designing hydrokinetic turbines, there should be a trade-off between reducing the price of manufacturing and maintaining or improving the performance of the turbine.

From a performance perspective, Darrieus turbines are typically installed in rivers flowing at velocities ranging from 0.5 to 2.5 m/s. These turbines spin at its rated speed (typically 100–150 rpm) to achieve the designed power coefficient while the generator turns at a different speed (usually higher) via a gearbox. The gearbox is used to manage the rotor speed and torque before the generator and is responsible for a large amount of energy loss due to friction. An additional weakness of Darrieus turbines is the issue of self-starting. There have been several studies conducted to evaluate the performance of hydrokinetic Darrieus turbine systems designs and configurations to resolve these issues.

In the work of Clarke et al. [33], a horizontal three-blade hydrokinetic rotor was coupled to a four-blade rotor. The hydrokinetic performances were evaluated, and it was reported that the overall coefficient of hydrodynamic performance was 0.43. Usui et al. [34] developed an experimental model where the same coefficient of hydrodynamic performance was observed. Both studies suggest that series-stream counter-rotating turbines (SS-CRTs) do not improve the hydrokinetic performance coefficient over conventional designs because the rotors work as a single system extracting energy from a single flow stream. Didant et al. [35] designed a direct-drive system for wind energy applications to improve the output performance and rotational speed over conventional designs. The reported results showed an overall improvement in output performance. Unlike SS-CRTs, parallel stream systems utilize the rotor capacitor more efficiently and the system swept area is the total of both rotor areas. This theory has been investigated in only a few studies and more investigations are needed.

In the work of Jannon and Boonsuk [36], the authors investigated the starting time profile of a direct-drive vertical axis Darrieus hydrokinetic turbine with an axial-flux permanent magnet generator. The turbine system was designed to reduce the overall frictional loss by eliminating the gearbox. The turbine and generator radiuses were varied, and the authors reported that the variation in radii influenced the total moment of inertia, the hydrodynamic torque, and the relative speed at which the winding coils cut through magnetic fields. A MATLAB code was developed for assessing the starting period and performance of the turbine design. It was reported that the starting period could be split into three phases: swing rotation, transient rotation and steady state rotation. It was concluded that the larger turbine radius produces more torque and shortens the starting time.

In the work of Jannon, a 1-D TMATLAB code was developed to assess the torque coefficient and ripple profiles for a hydrokinetic counter-rotating Darrieus turbine system of the same design. The author reported that friction loss and low angular velocity were the two major issues with both the wind and hydrokinetic Darrieus turbines. Hydrokinetic Darrieus turbine designs have been developed and reported incorporating two or more horizontal turbines placed together in series, where a geared mechanism enabled the turbines to turn in counter rotation with controlled angular velocity which only allows the two turbines work within one swept area.

There have been several studies conducted in relation to decreasing frictional losses in the rotor system, increasing the rotational speed of the rotors, and investigating starting profiles [37,38,39,40]. There have also been studies that have delved into the torque characteristics and power output profiles of single stream systems, however there are no studies to the authors knowledge that have gathered experimental or CFD data on the hydrodynamic performance of direct drive parallel-stream counter-rotating turbine systems [41,42,43]. There is a need for investigative work to improve the design of Darrieus turbine systems to extract energy from a stream of water efficiently. Although there have been many computational and experimental investigative attempts to develop single axis rotor systems [44,45,46], dual axis rotor turbine systems, and hydrokinetic Darrieus turbines placed in series to extract energy [47], there is little to no work on improving existing systems by developing a direct drive parallel axis Darrieus hydrokinetic turbine system that can extract energy efficiently from a stream of water.

In this work, we hypothesize that a direct-drive counter rotating system would aid in eliminating a large magnitude of friction losses while operating at low angular. We further hypothesize that a direct-drive system on a pair of vertical-axis Darrieus turbines that turn in counter-rotation to double the angular velocity could potentially extract more energy from a separate stream of water. Therefore, the system becomes a parallel stream direct-drive vertical axis Darrieus hydrokinetic turbine (DD-VADHT). A two-dimensional (2-D) and three-dimensional (3-D) computational fluid dynamic (CFD) simulation study was conducted to assess the hydrokinetic performance of the design. A 2-D and 3-D study will provide more accurate hydrokinetic performance data than 1-D studies. The next section provides an overview of the methodology of the present work.

2. Materials and Methods

A detailed overview of the DD-VADHT prototype development, PIV setup, and CFD numerical methodology is provided in the following sections.

2.1. DD-VADHT Prototype

A scaled prototype of the DD-VADHT design system was designed and developed for PIV experimental analyses in this work. A computer aided design (CAD) model of the scaled DD-VADHT prototype was developed using Solidworks (2020) CAD software and is shown in Figure 1. The PIV experiments were conducted on the prototype as a benchmark for validating the CFD simulations conducted in this work. The PIV studies were conducted in the University of North Carolina at Charlotte’s (UNC-C) Motorsports Center water channel. The dimensions of the total test section of the water channel are 1 m wide, 1 m tall, and 9 m long. The setup specification, settings, and procedures used in the PIV experimental studies along with results obtained are discussed. The DD-VADHT design prototype system is composed of two parallel axis rotors with four overall turbines (top and bottom) consisting of three NACA 0018 air foil blades placed 120° relative to each other as shown in Figure 1. The DD-VADHT prototype design parameters that were used in both experimental and computational studies are shown below in

Table 1. DD-VADHT prototype design specifications.

Description	Notation	Value
Chord length	c	1.3in
Turbine radius	R	8.5in
Turbine height	H	10in
Blade profile	-	NACA 0018
Number of blades	N	3
Tip speed ratio	λ	2

.

The DD-VADHT prototype is a (1:2) scaled model that includes two white cylinders in the center axis of the rotors as a surrogate for axial-flux permanent magnet generators. The cylinders are used to serve as a surrogate for geometry flow disruption and mounting point connections for the top and bottom turbine rotors. The turbines and generator surrogates were developed using an additive manufacturing process (3-D printing). The 3-D printer used to print the generator surrogates and the NACA 0018 turbine blades was a Project 3500 HDmax Printer. An image of the 3-D printed NACA 0018 turbine blades and generator surrogate are shown below in Figure 2.

A major design constraint of the DD-VADHT prototype from a PIV experiment perspective is that the size of the water tunnel test section restricts the scaling of the prototype. An additional design constraint is the blockage effect in the water channel and resulting wake at the turbine region of interest (PIV measurement area) area due to the channel walls. The prototype design consisted of two parallel counter-rotating turbine rotors with four overall turbines (top and bottom on each rotor). The rotors and turbines were position and enclosed in a silver anodized aluminum T-slotted frame rail cage system. Cross rails were added to the top and bottom of the cage system with corrosion resistant bearings to secure the rotors in place. The turbine blade supports are made of carbon fiber and were cut to shape and size using a Wazer waterjet desktop machine. The carbon fiber supports were used to hold the NACA 0018 blades in place with hex screws. The carbon fiber supports were secured in place with a carbon fiber cylindrical base and central connecting hub that connected the supports to the turbine rotor as shown in Figure 3.

All four turbines were fabricated using the same manufacturing process. A photograph of the assembled DD-VADHT prototype is shown in Figure 4. The black cylinder (pictured on the left side of Figure 4) and the white cylinder (pictured on the right side of Figure 4) are the 3-D printed generator surrogates.

As mentioned previously, the PIV setup for this work was conducted in the UNC-C Motorsports Center water channel. A schematic of the PIV setup is shown in Figure 5. The DD-VADHT scaled prototype was mounted inside the flow area of the water channel. A Dantec Dynamics Raypower 2000 laser is positioned under the test section generating a vertical laser light sheet that spans the streamwise direction of the flow after being routed 90° via a custom developed mirror assembly as shown in Figure 6. A Dantec Dynamic Flow Sense USB 2M-165 camera is also mounted on the underside of the water channel to capture images of the flow patterns around the turbines.

After installation and a preliminary PIV test run, design changes to the original scaled DD-VADHT prototype were required due to the size, unsteady nature of the flow inside the test section as a result of the rotating blades and channel walls, and reflection of laser light due to surrounding objects. The original prototype of the turbine system was designed for transient experiments where the turbine blades would rotate, and the flow pattern would be analyzed. Another major drawback observed were shadows created by surrounding objects on the prototype when illuminated by the laser. High resolution PIV data is only obtainable in areas where the fluid tracer particles are illuminated, and reflection and refraction are at a minimum. When the laser light hits a metal object on the prototype, a shadow is created behind the object and data capture in this area is difficult. This problem was rectified by coating the blades of the model with counter florescent orange paint. The paint was used in conjunction with a low pass filter on the camera lens to generate the clearest image of the particles while eliminating reflection. The upper turbines were also removed to avoid multiple objects to create scattering light and shadows. The revised physical model of the DD-VADHT installed in the water channel with the PIV laser is shown in Figure 6.

The PIV settings for the PIV configuration and experiment setting parameters for this work are shown in

Table 2. PIV configuration and setting parameters.

Step	Description	Parameter
Camera pixel	CMOS camera pixel specifications	1920 × 1200 pixels
Filter	Band pass filter fitted with the CMOS lens	532 nm
Inlet velocity	Average inlet velocity on the water tunnel	1 m/s
Laser energy	Energy emitted by the laser during acquisition	700–832 mJ
Particle type	Seeding particle diameter and type	1 µm Glass hollow
Frequency	Acquisition frequency of photographs	15 Hz
Time	Time interval between frames	2 ms

. The PIV settings are similar to the setting used in the work of Arpino et al. [48] and yielded the best results during trial runs. Opensource PIVlab (MATLAB) was used to process all PIV images using a Fast Fourier Transform (FTT) and cross-correlation technique. As mentioned previously, the rotors are fixed in place (not rotating) while water moves around the turbine blades. The PIV technique in this work is based on laser light illumination (pulsating) at twice the rate of the flow seeded by micron-sized glass hollow particles. The illuminated particles allow one to obtain two distinct images. The FFT and cross-correlation software allows the local velocity vectors to be reconstructed [49]. This technique also allows 2-D instantaneous velocity field measurements in a time interval that is short compared to the flow time scales [50]. Details about such non-intrusive measurement technique are available in the scientific literature [48,51,52,53,54] and are not reported here, but in the authors thesis [55].

PIV images of the turbine blades are taken in the water channel and PIV setup and are reported. The glass hollow particles and instantaneous velocity vectors are constructed using PIVlab software and are shown in Figure 7. Prior to analyzing images, a mask was created around the air foil blades and surrounding exterior features. For all PIVlab image processing, the FFT cross-correlation with multiple passes (two passes) and large deforming windows for enhancing the signal-to-noise ratio were used along with velocity vector smoothing and large magnitude velocity vector value outlier rejection. Figure 7a shows the reconstructed velocity vector and contour plots for a steady flow run in the water channel with a 1m/s flow setting. Figure 7b shows a 2-D CFD simulation steady state run with a 1m/s velocity boundary condition setting and a wall setting on the top and bottom and sides of the flow domain which has the same dimension of the water channel.

Figure 8 shows a plot of the PIV and CFD velocity profile extractions between the turbine blades from the region of interest (ROI) site of the CFD and PIV data. The plot also shows the standard deviation (error bars) for the PIV data as compared to the CFD data. The largest error occurs around the curves of the blades which is believed to be due to shadow regions and laser glaring that makes resolving velocity vectors a challenge. The results are within 84% agreement when comparing the overall average of the data points.

2.2. CFD Modelling and Meshing Methodology

The replication of the real-time flow physics observed in the PIV experiments for the DD-VADHT prototype CFD simulations along with obtaining quality numerical results fundamentally depend on the development of the correct geometric model, mesh generation, and appropriate boundary conditions. In this work, the commercial CFD package Star-CCM+ was used for flow model development, meshing, and numerical simulations. The geometry model of the DD-VADHT system was deployed from Solidworks into the meshing platform of Star-CCM+. Prior to deployment for meshing, faces of the CAD model were healed in sections where there are sharp edges and intricate features such as the aluminum T-slotted cage system, couplers and fasteners, the top half of the turbine rotors were removed. The couples and screws were also removed from the turbine rotors and blades to produce a simplified model for flow analyses as shown in Figure 9.

A top view of the three-dimensional fluid flow domain is defined and illustrated in Figure 10. The domain is split into a global stationary domain and subdivided into two rotating domains [56]. The rotating domains are defined by the simplified turbines as shown in Figure 9. The inlet, outlet, and surrounding faces of the stationary region are dimensioned at a length to accommodate the full development of the upstream and downstream flow from affecting the results for the analysis. The inlet and outlet boundaries are located both upstream and downstream from the centroid of the turbine system. A symmetry boundary condition is applied at the side, top, and bottom walls. For all flow simulations, water was modelled as the working fluid and the boundary conditions of the flow domain were set to: (1) 1 m/s for the inlet velocity, the pressure outlet was set to 0, and the rotating domains were set to counter-rotate at a speed of 120rpm for a TSR of 2.75 based on previous works [37,57].

As mentioned previously, the quality of the computational mesh directly influences the results in terms of the rate of convergence and grid independence. For this work, the cell sizes of the mesh were generated using an edge and face sizing in conjunction with a tetrahedral meshing scheme unstructured in both stationary and dynamic (rotating) regimes. The selection is based on the rationalizations that unstructured tetrahedral grids have the capabilities to discretize complex geometries with fast and minimum user intervention. The first coarse mesh produced for this work is shown in Figure 11.

For grid independence of the 3-D model, the first mesh was generated with a minimum face size of 0.004m and decreased in length scale increments of 20% until the delta change in terms of comparing current to previous mesh parametric results were minimized to less than five percent error. The same was also conducted for the 2-D simulations. For grid independence of the 2-D model, the power coefficient of the left and right turbines was used as the parameter of interest in evaluating and determining grid independence. The power coefficient is the ratio of power generated by the turbine to the entire amount of available energy flowing through the swept area. The grid independent study was initiated by selecting and running a coarse baseline mesh. The baseline element size was then reduced in steps until the difference in power coefficients was less than 1 percent. A 1 percent difference was chosen over a 5 percent difference due to the lower computational cost of a 2-D simulation.

Table 3. 2-D Model grid independent study results.

Mesh	Base Size	Refinement Ratio	Cell Count	Left CP	Left CP
Base	10.0in	-	85,995	0.09298	0.09291
M1	5.0in	1.53	131,750	0.09742	0.09772
M2	2.5in	1.61	212,454	0.10189	0.10276
M3	1.0in	2.01	427,539	0.10183	0.10243

shows the results for the 2-D model grid independent results.

Figure 12 shows the change in mesh size in terms of cell counts and the corresponding power coefficient for each steady state run for the model. As shown in Figure 11, the model becomes grid independent around 211,454 cells.

Figure 13 show images of the changes of the mesh around the hydrofoil blades. The number of prism layers was kept constant along with the first layer height and the total thickness of the prism layer. 40 prism layers were used with a total prism layer total thickness of 0.1 inches. This ensured that the boundary layer was contained in the prism layer for all mesh simulations to reduce variability. The first cell height off the wall was 5.0 ×10⁻⁶ m to ensure that the wall y-plus value was under 1. The target cell size on the surface of the blades was selected as 1 percent of the base. The minimum cell size was selected as 0.1 percent of the base size. The trailing edge of a hydrofoil requires a good quality mesh due to the turbulent flow in this area.

The mesh chosen for the CFD simulations is the M3 mesh. To ensure that the selected mesh was suitable, and that the boundary layer is resolved, the wall y⁺ value was calculated over the entire surface of the rotor for the 2-D and 3-D simulations. The prism layer of cells next to the walls was selected so that the wall y+ value was less than 5. The wall y⁺ value is a non-dimensional wall distance and is defined in Equation (1) below where u_τ is the friction velocity, y is the distance to the surface, and ν is the kinematic viscosity. The friction velocity can be defined as given in Equation 2 below with τ_w is the shear stress on the wall and ρ is the fluid density. The y⁺ value is calculated as:

$$y^{+}=\frac{u_{T}y}{\nu }$$

(1)

where,

$$u_T=\sqrt{\frac{{\tau }_w}{\rho }}$$

(2)

Figure 14 and Figure 15 show the y⁺ contours for both the 2-D and 3-D turbine blades. The maximum wall y⁺ values are on the edges of the hydrofoil blades which were areas that required mesh refinement. For turbulence modelling, the first cell should not fall outside of the log-law layer and the value should also not be small relative to the first node placed in the viscous sublayer within the boundary layer [58].

The next sections discuss the CFD flow modelling methodology.

2.3. CFD Flow Modelling Methodology

All CFD simulations in this work are evaluated using the Reynolds-Averaged Navier–Stokes (RANS) and k-ε solver in Star CCM+. The relations that describe conservation of mass, momentum and energy are provided and discussed in this section. The Reynolds decomposition process is discussed first. In this process each of the instantaneous variables in the governing equations can be separated into two components: a mean part and a fluctuating one with zero mean. For the velocity components:

$$u_i={\stackrel{\_}{u}}_l+u_i{}^{\prime }$$

(3)

for the pressure components:

$$P=\stackrel{\_}{P}+P^{\prime }$$

(4)

where u_i and P are instantaneous quantities, $\stackrel{\_}{u_i}$ and $\stackrel{\_}{P}$ are the mean parts of the variable and 𝑢_𝑖′ and 𝑃′ are the fluctuating components of the quantity. The average form of the conservation of mass and Navier–Stokes equations are shown below in cartesian form [59,60]:

$$\frac{\partial \rho }{\partial t}+\frac{\partial \left(\rho u_i\right)}{\partial x_i}=0$$

(5)

$$\frac{\partial \rho }{\partial t}+\frac{\partial \left(\rho u_i\right)}{\partial x_i}=-\frac{\partial P}{\partial x_i}+\frac{\partial }{\partial x_j}\left[\mu \left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}-\frac{2}{3}{\delta }_i{}_j\frac{\partial u_l}{\partial x_l}\right)\right]+\frac{\partial \left(-\rho \stackrel{\_}{u_l^{\prime }}\stackrel{\_}{u_j^{\prime }}\right)}{\partial x_j}$$

(6)

Equations (3) and (4) are the RANS equations where all the pressure and velocity terms involved are replaced by their mean value equivalents. Additionally, Equations (3) and (4) assume a viscous and compressible fluid. The last term on the right-hand side of Equation (4) however is not composed of mean velocity components. This expression, $-\rho \stackrel{\_}{u_l^{\prime }}\stackrel{\_}{u_j^{\prime }}$, is known as the Reynolds stress tensor. It represents the effects of turbulent motions (eddies) on the mean flow stresses. This tensor effectively adds six independent unknowns to the equations due to its three-dimensional nature, which prevents the computation of the RANS equations. This issue is known as the turbulence closure problem. The Reynolds stress tensor is defined in the following way:

$$\tau =-\rho \left[\begin{array}{ccc}\stackrel{\_}{u^{\prime }}\stackrel{\_}{u^{\prime }}& \stackrel{\_}{v^{\prime }}\stackrel{\_}{v^{\prime }}& \stackrel{\_}{u^{\prime }}\stackrel{\_}{w^{\prime }}\\ \stackrel{\_}{v^{\prime }}\stackrel{\_}{u^{\prime }}& \stackrel{\_}{v^{\prime }}\stackrel{\_}{v^{\prime }}& \stackrel{\_}{v^{\prime }}\stackrel{\_}{w^{\prime }}\\ \stackrel{\_}{w^{\prime }}\stackrel{\_}{u^{\prime }}& \stackrel{\_}{w^{\prime }}\stackrel{\_}{v^{\prime }}& \stackrel{\_}{w^{\prime }}\stackrel{\_}{w^{\prime }}\end{array}\right]$$

(7)

Many turbulence models were created to model the Reynolds stress tensor term and solve the RANS system of equations. The Boussinesq approximation assumes that the momentum transfer caused by the turbulent eddy interactions can be modeled with the concept of an eddy viscosity. This is analogous to how the momentum transfer resulting from the molecular interactions inside a gas can be related to a molecular viscosity. The Boussinesq approximation provides a new equation for turbulence closure and is formulated in the following way:

$${\tau }_{ij}=2v_t{S}_{ij}-\frac{2}{3}k{\delta }_{ij}$$

(8)

where τ_ij is the Reynolds stress tensor, v_t is the eddy viscosity, k is the turbulent kinetic energy, and is the Kronecker delta and S_ij is the strain-rate tensor of the mean field expressed as:

$$S_{ij}=\frac{1}{2}\left(\frac{\partial {\stackrel{\_}{u}}_i}{\partial x_j}+\frac{\partial {\stackrel{\_}{u}}_j}{\partial x_i}\right)$$

(9)

As mentioned previously, the standard k-ε turbulence model is also used in this work for CFD simulations. This model is the simplest of turbulence modelling and is a two-equation model in which the solution of two separate transport equations allows the turbulent velocity and length scales to be independently determined. The turbulence kinetic energy, k, and its rate of dissipation, ε, are obtained from the following transport equations [56,61,62]:

$$\frac{\partial \left(\rho k\right)}{\partial t}+\frac{\partial \left(\rho k{u}_i\right)}{\partial x_i}=\frac{\partial }{\partial x_j}\left[\left(\mu +\frac{{\mu }_t}{{\sigma }_k}\right)\frac{\partial k}{\partial x_j}\right]+G_k+G_b-\rho \epsilon -Y_M+S$$

(10)

and

$$\frac{\partial \left(\rho \epsilon \right)}{\partial t}+\frac{\partial \left(\rho \epsilon u_i\right)}{\partial x_i}=\frac{\partial }{\partial x_j}\left[\left(\mu +\frac{{\mu }_t}{{\sigma }_{\epsilon }}\right)\frac{\partial \epsilon }{\partial x_j}\right]+C_{1\epsilon }\frac{\epsilon }{k}+\left(G_k+C_{3\epsilon }{G}_b\right)-C_{2\epsilon \rho }\frac{{\epsilon }^2}{k}+S_{\epsilon }$$

(11)

In these equations, G_k represents the generation of turbulence kinetic energy due to the mean velocity gradients, G_b is the generation of turbulence kinetic energy due to buoyancy, and Y_M represents the contribution of the fluctuating dilatation in compressible turbulence to the overall dissipation rate. C_1ε, C_2ε, and C_3ε are constants, σ_k and σ_ε are the turbulent Prandtl numbers for k and ε, respectively. S_k and S_ε are user-defined source terms. The turbulent eddy viscosity, µ_t, is computed by combining k and ε as shown below:

$${\mu }_t=\rho C_{\mu }\frac{k^2}{\epsilon }$$

(12)

where C_µ is a constant. The model constants C_1ε, C_2ε, C_µ, σ_k and σ_ε are 1.44, 1.92, 0.09, 1.0 and 1.3, respectively. More information can be found in the Star CCM+ user manual. The next section presents the results of this work. The next section presents the results of this work.

3. Results

This section of the paper provides a discussion of the CFD results obtained in this work. The authors were mindful that it is difficult to validate 2-D and 3-D CFD results with PIV experimental measurements. Especially in the case of a Darrieus turbine system with high solidity and operating at TSRs above 1. Nevertheless, the results presented in the previous section provided a baseline understanding of the physics behind this turbine geometry and prepared the authors for improved numerical modelling and a high confidence level with the results presented in this section.

The performance comparison characteristics of the DD-VADHT prototype design was predicted in CFD based on power coefficient, pressure coefficient, self-starting capability with torque coefficient analysis, and flow behavior around the turbines. We begin this analysis with a 3-D CFD simulation of the DD-VADHT. Figure 16 shows the velocity iso-contours of the flow around the turbine blades on the right rotor. From Figure 16, a high-speed zone is suggested on the inner region of the rotating domains where the two hydrofoils meet. There are also low-speed regions on the outer regions of both rotating domains. Although the 3-D simulation provides a three-dimensional perspective of the flow patterns and better insight of spatial dependent parameters such vorticity, for the sake of computational expense, this work evaluates a 2-D model for performance evaluation of the DD-VADHT.

Figure 17 shows a plot of the variation of C_P calculated with tip speed ratio (TSR) for the NACA 0018 hydrofoil blades of the DD-VADHT obtained from the 2-D simulation. It can be seen in Figure 17 that the maximum C_P is found to be approximately 0.48, on the similar lines of the results obtained by Kiho et al. [63]. The DD-VADHT is designed for optimal efficiency at a constant revolution per minute (RPM) of 90 rpm, relating to the TSR of 2.75 at an incoming flow of 1 m/s.

A plot of C_P variation over a range of TSR values was analyzed for the selection of the best TSR, leading to the optimal performance of the NACA 0018 hydrofoil profile on the DD-VADHT turbines, which was used in the blade geometrical design. The C_p versus azimuthal position of different TSRs during rotation (azimuth angle variation) is shown in Figure 18. All the TSR curves are periodic and appear to have an increasing C_p until 50 degrees is reached and then a decreasing C_P approximately until 90 degrees. A nearly sinusoidal curve was obtained with three positive maxima and minima (for a TSR of 2.75) or negative minima. This indicates that during a revolution, there are three periods of time where the turbine produces torque on the fluid. From Figure 19, it can be assumed that the maximum torque is achieved at an azimuth angle around 50-degrees. After the peak, the drag begins to increase as the blade enters a dynamic stall, and the drag starts to be dominant up to an azimuth angle of about 135-degrees. For all plots of C_p vs. TSR, the plots show positive values of TSR, meaning that fluid is providing torque to the turbine.

An interesting phenomenon that occurs in Figure 19 is a secondary peak that occurs after the primary peak between the peak-to-peak regions of azimuth angle positions 45-degrees and 165-degrees. This is believed to be attributed to low pressure recirculation zones that occur between the blades as they meet at azimuth angle 120-degrees. This can be better visualized by pressure coefficient contours which is discussed later. Figure 20 shows the C_P vs. varying azimuth angle for three individual blades on the right turbine.

As shown in Figure 20, the maximum C_P produced is achieved around 68°for the first blade. After the peak, the drag begins to increase as the blade enters dynamic stall and the drag starts to be dominant up to an azimuth angle of 195°. The second blade has the opposite trend in C_P as compared to the first blade and the power production is completed with the same motion of the third blade relative to the first blade (at a phase shift) for one rotation of the turbine. The plot of C_P for a TSR of 2.75 shows positive values indicating that the that the fluid is providing torque to the turbine as mentioned previously. This can be seen in Figure 21 when analyzing the torque coefficient C_T. The torque output of the DD-VADHT design is crucial data for system development and performance analysis. There is little to no data available in literature on torque output for DD-VADHT systems and similar designs. As shown in Figure 21, the torque coefficient shows constant periodic values over the entire revolution of the turbine similar to the C_P plots discussed. The same trend in primary and secondary peaks are present as compared to the C_P plots which may be due to the turbines operating in two swept areas.

There have been several studies that have investigated the torque characteristics of single Darrieus rotors. In these cases, the torque coefficient vs. azimuth angle plots were reported to be sinusoidal with peaks at 120° apart (3-blade rotors) with no primary and secondary peaks like the C_P plots reported in this work. Additionally, the magnitude of fluctuation was represented by ripples, which were as high as 600% (helical rotor) as reported in the work of Niblick [64] and 920% Hall [65]. The ripple due to torque can be calculated using Equation (13) below:

$$\left(\frac{C_{T\max }-C_{T\min }}{C_{Tavg}}\right)\times 100\%$$

(13)

In the case of this work, the ripple was found to be 74% for the three-blade design in the DD-VADHT. This suggests excellent torque smoothing and possibly good self-starting capabilities. To characterize the pressure characteristics, the pressure coefficient contours at 105, 120 and 135-degrees (positions where the two hydrofoils begin to meet) for one single blade each of the right and left turbine rotor are provided in Figure 22.

As shown in Figure 22, the pressure coefficient range trends with the change azimuth angle are identical. The pressure fluctuations become more intense in pressure coefficient magnitude with the increase of the spanwise azimuth angle at 135 degrees. At the cross sections of the 120-degree angle, the pressure fluctuations weakened for both the barotropic region (red areas) and the negative pressure region (blue areas). From the pressure contours of the three angles analyzed, it is noticeable that the maximum absolute values of the pressure coefficient concentrate on the region of the leading edge of the airfoil and the suction surface of the airfoil. Unlike the other two cross sections of the hydrofoil, the negative pressure region of the 120-degree azimuth angle spread over the suction surface from the leading edge to the trailing edge. The negative pressure regions indicate wake regions.

Additionally, the pressure coefficient on the outer area of the profile of the rotating domains is observed to be larger (intrados) than that of the inner zone of the rotating domains (extrados). The difference in pressure and velocity causes the overall lift for the hydrofoil blades. The pressures are negative near the edges of the blade and positive at the exit area of the computational domain (to the right of the contours). A snapshot of the pressure coefficient contours for the right turbine rotating domain region is shown in Figure 23 for various azimuth angles to illustrate the pressure coefficient trends at varying angles. This pressure coefficient contours in Figure 23 highlight the intrados and extrados as well as the suggested wake regions. A better observation of wake and recirculation regions can be obtained from viewing velocity and vorticity contours. The velocity contours are provided in Figure 24.

The high-speed zones are found at the tip of the blade leading along the chord length facing the interface of the rotating and stationary domain. The velocity of the fluid is reduced at the boundary of the blade profile near the regions where the two rotating domains meet. Further, it has been found that the pattern of the pressure and velocity contours for each blade is similar for each profile. There also appears to be a region of minimized flow velocity magnitude (compared to flow around the foils) downstream (to the right in the contours) of where the two blades separate at 135°. From a design perspective, this is a major factor that ensures that the turbine has good flow velocity reduction performance. The fact that the turbines operate in two swept areas allows the extraction of more energy from the flow field as opposed to one swept area.

4. Conclusions

The goal of this work was to conduct a computational design analysis study of a hydrokinetic horizontal parallel stream direct-drive (no gear box) counter-rotating Darrieus turbine system. The design analysis required to reduce the computational domain via removing intricate features of the computational model that introduced meshing issues such as the aluminum T-slotted cage system, couplers and fasteners, the top half of the turbine rotors, and introducing filets in sharp corners and taking advantage of symmetry where possible. Essentially, the simulation setup was reduced from a 3-D to a 2-D analysis. Results such as power coefficient, pressure coefficient, flow pattern behavior, torque coefficient, and ripple effect which provides information on torque smoothing were obtained and discussed. The ripple effect reported in this work was 74% which is much less than the average reported in previous works that used one swept area. The flow analysis also showed that the hydrofoils generate sufficient lift to rotate and extract the energy from the flow and produce the optimum power needed.

Additionally, the study has demonstrated the proposed design is an improvement over previous hydrokinetic Darrieus systems that have incorporated one rotor. However, further work is needed to model different design system configurations such as the number of blades and lag angle in order to determine optimal system performance. Additionally, additional experimental work is needed to validate the elimination of friction by removing the gear box from the DD-VADHT prototype.

The presented DD-VADHT design produces the following advantages:

• It appears to be capable of extracting more energy from a double swept area than a single swept area by comparing to previous designs when analyzing ripple effect, torque coefficient, and power coefficient.
• It promises good torque smoothing and possibly good self-starting capabilities, although more work should be done (i.e., CFD, one-dimensional data code analyses and experimental) to verify this.
• Additionally, it shows good flow velocity reduction performance for smooth flow transitions between blades and generating sufficient lift.

Lastly, the following will be considered for future work (Phase 2 Study):

• A 3-D CFD study should be performed to analyze the vorticity flow pattern behavior between the rotating domains. It is difficult to provide conclusive remarks of vorticity in 2-D simulations as vorticity varies in three dimensions. Better conclusions can be drawn from 3-D simulations regarding vorticity.
• A Large-Eddy Simulation (LES) study will be achieved for analyzing recirculation regions and eddies around curved surfaces and between the rotating domains.
• Analyze multi-blade profiles and lag angles to access performance output using a 2-D simulation over a range of RPMs and inlet flow velocities. A comparison analysis will be helpful in assessing ripple effect for different design configuration and design optimization.