Model of a Ducted Axial-Flow Hydrokinetic Turbine – Results of Experimental and Numerical Examination

Abstract The article presents the numerical algorithm of the developed computer code which calculates performance characteristics of ducted axial-flow hydrokinetic turbines. The code makes use of the vortex lattice method (VLM), which has been developed and used in IMP PAN for years, to analyse the operation of various fluid-flow machines. To verify the developed software, a series of model tests have been performed in the cavitation tunnel being part of IMP PAN research equipment.


INTRODUCTION
The process of electric energy generation in hydrokinetic turbines assumes the use of kinetic energy collected in water currents of rivers or other inland waterways, as well as in sea currents. Unlike conventional water turbines, hydrokinetic turbines do not require artificial water damming systems, as a result of which the operation of these devices disturbs the natural flow of the water current to a minimal extent.
Studying possibilities of the use of hydrokinetic energy began in the 1970s, after the energy crisis caused by the embargo imposed by OPEC countries on the USA and Western European countries [24]. Dramatic price increase of crude oil in those times, followed by similar price increase of coal and natural gas, led to the stagnation and recession of the world economy. Countries importing energy-producing raw materials were forced to intensify the research oriented on modernisation of technology of energy production from conventional sources and search for new energy sources. At present, a basic motivation for taking actions which aims at reducing the percentage of mining resources in the entire energy production in favour of ecologically clean technologies based on renewable energy sources is environmental protection. That is why great interest towards hydrokinetic turbines is observed, and wide-ranging studies of devices of this type are performed worldwide. The overview of existing commercial designs of systems which utilise hydrokinetic energy to produce electric current, and their classification can be found in [18,24,29].
Research studies upon hydrokinetic turbines are performed in two ways: (1) as theoretical analyses, based on computational fluid dynamics and making use of relevant computer codes, and (2) as experiments, performed both on laboratory scale on research rigs, and in real, natural conditions. In the former case a number of calculation methods are used. One of the most popular methods used in the fluid-flow machine design and optimisation process is the Blade Element Momentum (BEM) method. It was initially developed by H. Glauert for airplane propellers [8], as a combination of the liquid flow theory (operating disc model) developed by W. Froude and W. J. Rankine [27] for marine propellers, and the model of blade segment developed, separately, by S. Drzewiecki and W. Froude [7]. A description of the BEM method applied to wind and hydrokinetic turbines can be found in [3,13,20,23,31].
Another group of methods used to analyse the operation of fluid-flow machines comprises methods which are based on surface distribution of singularities. Those methods apply the model of perfect fluid and assume the potential nature of the flow. As a consequence, the equations of motion are reduced to the Laplace equation, the solution of which is a superposition of elementary solutions representing the uniform flow and hydrodynamic singularities of vortex, source, and/or dipole type. A broad description of those methods, with their differentiation depending on the type of applied singularities and boundary conditions, can be found in [28]. This group of methods includes the Vortex Lattice Method (VLM), used in this publication, which is widely applied to analyse the operation of such fluid-flow machines as: marine propellers [21,22], centrifugal pumps [17], wind turbines [6,14], and hydrokinetic turbines [12,25,32].
In recent years, the increasing computing potential of newly built computers has provided good opportunities for the use of advanced and demanding, in terms of the required computer parameters, methods to study the flow through hydrokinetic turbines. These methods consist in solving the Reynolds Averaged Navier-Stokes (RANS) equations [16,32]. Very popular solvers, such as ANSYS/FLUENT, ANSYS/CFX, or NUMECA/Fine, which are used for this purpose make use of the finite volume method. These widely applicable solvers are used to perform comprehensive numerical analyses of issues related to hydrodynamics and compressible fluid mechanics, taking into account multiphase and transonic flows.
Along with theoretical-numerical analyses of the flow through hydrokinetic turbines, experimental investigations are also performed, mainly in laboratory conditions. Based on the available literature, two types of laboratory research rigs can be named. The first type is comprised of closedcircuit water tunnels [2,11]. In those tunnels, the water flow is generated by circulation pumps and the water flows in a closed space in which the examined model turbines are placed. In the other research rig type, the model turbines are placed in open water tunnels, or towing tanks (channels) with free water surface [4,15,30]. Experiments are also performed in real, natural conditions [5].
The object of numerical and experimental examination, the results of which are presented in this article, is a model axial-flow hydrokinetic turbine equipped with a duct as part of the flow system. The duct improves the energy efficiency of the turbine by increasing the velocity of the water flowing through the rotor. These solutions are sometimes used in hydrokinetic turbines [19,30] and wind turbines [1,26].
As compared to the earlier work by the authors [12], this article presents new elements introduced to computing procedures of the VLM based software developed in IMP PAN in recent years. The introduced elements mainly consist of: • extending the VLM by a procedure which calculates the shape of the vortex wake downstream of the rotor [9], • distributing the vortex lattice on both the suction and pressure sides of the rotor blade [9], instead of the average surface of the blade, a method which was used in earlier versions, • introducing a numerical procedure which calculates the effect of walls bounding the measuring space of the cavitation tunnel on the obtained results of calculations [10]. As a result of the above actions, better compliance between the theoretically predicted performance characteristics of the examined model hydrokinetic turbines and those experimentally recorded in the IMP PAN cavitation tunnel was obtained.

MODEL TESTS
The examined physical model of hydrokinetic turbine consists of a 5-blade rotor with a hub, closed in a duct -see Fig. 1. According to the adopted assumption, both the rotor blades and the duct are 3D shaped. The model turbine was examined in the IMP PAN water tunnel. The research rig is schematically shown in Fig. 2. The cross section of the measuring chamber has the square shape of size 0.425 m.
The measuring instrumentation of the rig enables to measure: the average velocity ∞ V of the water flow in the tunnel chamber, the turbine rotor revolutions n, the torque M on the turbine shaft, and the axial thrust T generated by the turbine rotor.
The average water velocity ∞ V in the chamber is measured indirectly, using a calibrated tunnel confusor. A quantity which is measured directly for this purpose is the water head difference (∆H) between the confusor inlet and exit cross sections. The measurement is done using a precise differential manometer or pressure transducer, and the average water velocity is calculated using the formula obtained from calibration. The turbine rotor rotational speed n is measured using the pulse frequency meter which counts pulses generated by a disc with properly machined grooves.
The shaft torque and the hydraulic thrust of the examined model turbine were measured using dynamometers, installed in the rig, which were equipped with special balancing scales. In order to eliminate the effect of shaft resistance on the torque measured during the tests, an additional measurement was performed, which aimed at experimental determination of the dependence of the shaft line resistance torque on the rotational speed, in the case when the turbine rotor is dismantled. The model turbine torque which was then measured in the basic tests was corrected by the value of the shaft line resistance torque corresponding to the current rotational speed.
In case I the rotor was placed exactly in the inlet section of the duct, while in case II it was shifted by 5 cm towards the duct exit - Fig. 3.
The tests were performed for comparable parameters of water flow in the tunnel, with the average water flow velocity approximately equal to V ∞ =3.4 m/s. The required rotational speed of the examined turbine model rotor was set using the electric current generator load control system.

Fig. 2. Scheme of the research rig (cavitation tunnel) [11], 1-turbine model, 2 -torque measurement, 3 -thrust measurement, 4 -measurement of pressure difference at measuring confusor inlet and exit, 5 -rotational speed measurement, 6 -circulation pump
The object selected for tests described in the article was a 5-blade rotor of D = 148 mm in diameter -detailed data in Tab  The performed tests made the basis for calculating the below defined dimensionless coefficients: Water energy utilisation factor: Thrust coefficient: TSR -tip speed ratio: The results of the performed tests, along with their uncertainty, are shown in the latter part of the paper (Figs. 9 and 10), and compared to the results of calculations obtained using the developed code.

CALCULATION METHOD
The calculation method developed to analyse the flow through the examined hydrokinetic turbine model ( Fig.  1), bases on the VLM, which, as mentioned above, is one of surface singularity distribution methods. It assumes that in the entire analysed domain the flowing fluid is incompressible and inviscid, except some singular fragments. After assuming that the flow is steady and at the same time without vorticity, , the equations of motion are reduced to the Laplace equation: where f is the scalar function which bears the name of velocity potential and describes the velocity field as: (5) In the case of liquid flow through the examined turbine, the areas where singular points are located are: (1) boundary layers of solid elements (duct, rotor blades, hub) coming into contact with the fluid, and (2) vortex wakes forming downstream of these elements. In those areas, dominated by viscous forces, strong vorticity concentration is observed [17]. The adopted method consists in discrete distribution of the so-called horseshoe vortices in the abovementioned singular areas [22]. The horseshoe vortex consists of the bound vortex, free vortices, and trailing vortices modelling the vortex wake forming the downstream of the washed turbine elements due to the motion of the boundary layer - Fig. 4. This way, the continuous vorticity distribution in the flow is substituted by a discrete distribution of vortex filaments.
As mentioned in the introduction, certain modifications have been introduced to the initial version of the developed in-home code the results of which had been presented in [12]. Firstly, the vortex wake, modelled by trailing vortices, had initially been approximated by helices with the assumed pitch angle. Proper selection of this angle was essential for the correctness of the obtained results, and it required huge experience from the code user. Therefore a decision was made to determine the shapes of the trailing vortices in a numerical, iterative way. Justification for this approach is given further in the article, in the form of a comparison analysis.
The next innovation introduced to the earlier code version referred to the way of distribution of the vortex lattice over the rotor blades. Initially, the vortices had been distributed on the average blade surfaces. However, during attempts to modify the code it was found, based on comparison of numerical and experimental results of ship propeller tests [9], that the distribution of vortices over the suction and pressure side produce more realistic pressure distributions. This opinion became the motivation for introducing the above change.
The third essential modification was extending the algorithm with a correction taking into account the effect of the measuring chamber walls on the characteristics of the turbine model placed in this chamber during laboratory tests. This correction was introduced in the way described further in the article. Its introduction was necessary, as neglecting this effect had led to remarkable overstatement of model turbine performance, compared to the machine operating in unbounded space.
The developed calculation method neglects flow turbulence and some other phenomena, in particular those resulting from fluid viscosity, such as flow separation for instance. In the presented calculation task, first the circulations Γ of the vortex filaments distributed over the surfaces of fluidflow elements of the examined machine are to be determined. For this purpose the Neumann condition was used, which says that normal velocities to the surface of the washed body are zeroed. This condition can be written as follows: (6) In the examined case this conditions is expressed as the sum of projections of the undisturbed flow velocity and the velocities induced by vortices of circulation Γ onto the direction normal to the wall at control points situated between the bound vortices: . The impact factor Kn ij is interpreted as a projection of the velocity induced by the j-th horseshoe vortex with unit circulation onto the normal direction at the i-th control point. Its formula is: The code assumes constant and uniform velocity of the fluid flowing into the turbine. As a consequence, the computational domain can be divided into characteristic segments, in which the circulations of the vortices situated at the same (corresponding) points are the same. This assumption leads to a remarkable reduction of the equation system and shortens the computational time. In this case, one main segment of the computational mesh is selected, and control points are placed within this segment. A sample single segment of the turbine flow system is shown in Fig. 5.

Fig. 5. Single segment of the computational domain
When applying this simplification, the impact factor formula in the equation system is to be modified to the following form: (10) where n is the number of the segment with a single blade, Z is the total number of blades, i' is the sequential number of the control point in the selected main segment, and j' is the sequential number of the vortex in the n-th segment.
Once the circulations of the horseshoe vortices are known, the velocities tangential to the walls of the analysed flow element are calculated from the formula: (11) where represents the tangent vector to the surface of the washed element at the i-th control point, is the width of the surface element on which the bound vortex is situated. The sign "-" refers to the suction sides of the rotor blade surface, hub, and duct while the sign "+" refers to the pressure side of the rotor blade surface. The matrix of Ks coefficients is calculated from the formula: Based on the Bernoulli equation, for the obtained tangential velocities the dimensionless pressure coefficient values are calculated as: (12) where ∞ p is the pressure in the undisturbed flow, i p is the pressure at i-th control point, and, ρ is the density of the liquid. The obtained distribution of the dimensionless pressure coefficient is used to calculate the torque M of the turbine rotor and its axial thrust T. These quantities are calculated from the following integral equations: In these equations S represents the area of the rotor blade surface, r z , r y , are components of the vector connecting the rotation axis with the control point, t x , t y , t z and n x , n y , n z are components of the tangent and normal vector to the rotor blade surface, respectively, Cp s and Cp p are the dimensionless pressure coefficients on the suction and pressure side of the rotor blade, respectively, and C f is the viscous drag coefficient, defined, according to [21], as: (15) The obtained results of calculations enabled to determine the dimensionless values of the water energy utilisation factor C E , thrust coefficient C T , and tip speed ratio TSR, defined by the formulas (1), (2) and (3), respectively.
The experience gained when developing the software here described revealed high sensitivity of the obtained values of the above coefficients to the assumed shape of the trailing wake, which, as already mentioned, in the initial versions of the code had been approximated by helices with the assumed pitch angle. In the recent version, a decision was made to determine the shape of the trailing wake in an iterative way. For this purpose a procedure was developed, the basic concept of which is given in Fig. 6. Fig. 6. Scheme illustrating how the shape of vortex wake is determined [31], -node points of the wake line, -auxiliary point, local flow velocities, calculated at the node point and at the auxiliary point , respectively The calculation process begins with introducing the vortex wake in the form of helices with an assumed pitch angle to the computational domain. The shapes of these lines are approximated by a sequence of rectilinear line segments connected at nodes. In further iterations, the shape of the vortex wake is modified by new vortices shed sequentially from trailing points situated on particular turbine elements (rotor blades, hub, and duct). The coordinates of the vortex generation points, being the first nodes of particular vortex filaments, are kept constant in subsequent iterations, while the locations of the following nodes are found based on local flow velocity vectors, calculated as the sum of the undisturbed flow velocity and the velocities induced by the system of vortices having the shape determined in the previous iteration: The calculation process for a single mesh node is shown in Fig. 8. It consists of two main steps. In the first step, according to formula (16), the local velocity is calculated at the earlier node , making the basis for determining the point from the following formula: (17) At this step, the local flow velocity is also calculated from (16). In the second step, the coordinates of the trailing wake node , are calculated using the average flow velocity vector obtained from the earlier calculated vectors and as: (18) In order to improve the stability of the iteration process, the calculation procedure has been complemented by the relaxation formula having the following form: (19) where i * is the node number, s is the iteration number, and ε is the weight coefficient, with the value of 0.5 assumed in the reported calculations.
Each time, for the newly determined shape of the vortex wake, the dimensionless values of water energy utilisation factor and thrust coefficient are calculated. When the values of those coefficients differ by less than 0.2% or oscillate around a constant value in successive iterations, the calculation process is terminated.
The calculation algorithm described above, and the resulting developed code, can be used to calculate performance parameters of hydrokinetic axial-flow turbines working in unbounded space. Since the experimental examination of the model turbine was performed in the bounded tunnel space, the calculation method had to be properly extended to improve the compliance of the numerical results with the experiment. The method extension by a correction which takes into account the earlier neglected effect of walls is described below.
The algorithm used to calculate this correction bases on the image method [7]. In Fluid Mechanics this method is used for modelling the flow past a body (here: through a hydrokinetic turbine) at the presence of surfaces bounding the flow area, when the flow is modelled using a system of hydromechanical singularities, vortices for instance. The method consists of substituting the Neumann boundary condition Vn = 0 on immobile walls by the images of systems of singularities having proper strengths and situated at proper points outside the analysed flow area, to get the resultant normal velocity component induced by the system of singularities and their images equal to zero. The method is illustrated in Fig. 7, for a single vortex in the flow bounded by a single wall. The correction taking into account the effect of the presence of walls bounding the measuring space of the cavitation tunnel is included in the code by adding the velocities induced by the images of the vortex systems modelling the flow through the turbine. This is done when the coefficient of the equation system (8) and the tangential velocities are calculated, and the new shape of the trailing wake is determined.
In order to take into account the effect of walls bounding the measuring space of the tunnel, the wall impermeability condition (7) was extended by adding the velocities induced by the images, to get the following form: (20) where stands for velocities induced by images at the i-th control point.
After applying the above equation to all control points, one obtains the equation system in the matrix form, which can be written as: (21) where is the impact factor matrix from images.
The earlier code version assumed constant velocity of the fluid flowing into the turbine, which enabled to reduce significantly the number of unknown variables due to periodicity of turbine geometry. Although justified for flows in unbounded areas, this assumption is too simplistic for the turbine flow analysed here. Indeed, the examined turbine model was equipped with a 5-blade rotor and a duct being a body of revolution, but the test chamber had a square cross section, which made the entire system aperiodic. In order to preserve the earlier reduced equation system and, simultaneously, keep the accuracy of calculations at an acceptable level, a decision was made to average in radial direction the velocity field generated by the created images. The quantity which was averaged for this purpose was the projection of the velocity vector onto the normal direction to the surface of the washed element at the control points marked by the same sequential number in each system geometry segment. The impact factors from images were calculated using the following formula: (22) where k is the sequential number of the segment in which the control point is situated, p is the sequential number of the image, n is the sequential number of the segment in which the vortex is situated, L_O -is the number of images, and: • is the vector connecting the control point i ' situated in segment k with the vortex j ' situated in segment n of image p, • is the vector normal at the control point i ' of segment k.
The tangential velocity calculation also includes velocities induced by the images, calculated from the formula: (23) in which L_W_S stands for the number of vortices in a single segment, Γ i' is the circulation of the vortex situated at control point i ' , is the distance between bound vortices, is the tangent vector to the surface of the washed turbine element at control point i ' of segment k. The values of the coefficient matrix K t are calculated from the formula: Note: the sign in front of the last term in formula (23) is assumed in the following way: "-" for the suction side of the rotor blade surface, and the hub and duct surfaces, "+" for the pressure side of the rotor blade surface.
The velocities induced by the created images were also taken into account when calculating the vortex wake. In that case, the local velocity at mesh nodes in the area occupied by the wake was calculated using the following formula: (24) As mentioned above, flow aperiodicity in the analysed case is the reason why the shape of the trailing wake shed from each turbine segment is not the same. In order to preserve the reduced equation system, the velocities induced by the images were calculated at wake node points in each geometry segment. Then, they were rotated by a given angle in such a way that the calculated vectors were situated at the mesh node of the main segment. Finally, the components of the calculated velocity vectors were averaged.

COMPARISON OF EXPERIMENTAL AND NUMERICAL RESULTS
The results of experimental tests and numerical calculations are compared in Fig. 9 for case I, when the rotor is situated in the inlet plane of the duct, and in Fig. 10 for case II, when it is shifted downstream of the inlet by 5 cm. The results are presented in dimensionless form, as water energy utilisation factors C E and thrust coefficients C T , calculated for selected TSR values. Based on the experimentally recorded points, for which total measuring uncertainties were also marked, the approximating curves were constructed using the method of least squares. The distance between the numerically calculated points and this curve is considered the measure of the difference between calculations and experiment. Based on the obtained results, a general conclusion can be formulated that for both examined turbine model cases good compliance between the measured and calculated values of the dimensionless energy utilisation factor, C E , has been obtained. The largest discrepancies can be observed for extreme TSR values. In particular: for case I relatively large differences are observed at characteristics points representing maximal TSR values. For case II, on the other hand, the calculated values of both coefficients are slightly overestimated for high TSR values, as compared to the experiment based approximating curve, while the C E values are slightly underestimated for small TSR values.
At the present stage of research, only probable causes of the observed discrepancies can be indicated. A very likely cause, especially in the areas of lower C E values or those going beyond the optimal operating points of the examined models, is significant simplification of the developed calculation method, mainly resulting from neglecting fluid viscosity and flow separation phenomena in the model. What is also noteworthy is the vibration of the turbine model rotor observed during the experiment. This vibration was much higher than that the one observed when the rotor was removed to measure the dependence of the rotor shaft resistance on rotational speed of the turbine without rotor. Therefore, the real resistance torque of the vibrating shaft with rotor could be slightly higher than the values assumed in the calculations and, consequently, the resultant difference could be larger for higher TSR.
For the dimensionless thrust coefficient C T in case I, large discrepancies between the calculated and experimental values were observed for high TSR. In case II, relatively good compliance of C T values obtained from calculations and experiment was observed. A probable cause for the observed differences is disregarding separation and cavitation phenomena in calculations.

BASIC CONCLUSIONS
Comparing the results obtained from the experiment with those calculated numerically, one can observe: • good compliance of values of the dimensionless water energy utilisation factor C E in the vicinity of the optimal point of turbine operation, • relatively large differences of the dimensionless thrust coefficient for high TSR values in case II. Good compliance between the calculated and experimentally recorded values of the dimensionless energy utilisation factor C E and thrust coefficient C T within the range of optimal operating conditions of the hydrokinetic turbine model suggests that the developed software can be effectively used in designing turbines of this type.
The presented results of computer simulations of the flow through the hydrokinetic turbine neglect such phenomena as flow separation and cavitation. Further steps of software development will include preparing relevant procedures which will enable to take into account these phenomena, which were likely to reduce the presently observed differences between experimentally recorded and calculated values of coefficients C E and C T .

NOMENCLATURE
c -profile chord length C E -dimensionless water energy utilisation factor C f -viscous drag coefficient C p -dimensionless pressure coefficient