An Unsteady Boundary Element Model for Hydrodynamic Performance of a Multi-Blade Vertical-Axis Tidal Turbine

: An unsteady boundary element model is developed to simulate the unsteady ﬂow induced by the motion of a multi-blade vertical axis turbine. The distribution of the sources, bound vortices and wake vortices of the blades are given in detail. In addition, to make the numerical solution more robust, the Kutta condition is also introduced. The developed model is used to predict the hydrodynamic performance of a vertical axis tidal turbine and is validated by comparison with experimental data and other numerical solutions available in the literature. Good agreement is achieved and the calculation of the proposed model is simpler and more efﬁcient than prior numerical solutions. The proposed model shows its capability for future proﬁle design and optimization of vertical axis tidal turbines.


Introduction
Over-exploitation of natural resources has caused serious ecological and environmental problems globally, threatening human civilization and the world economy [1]. A popular trend in resolving this issue is to increase the use of renewable energy. Ocean tidal energy, as one kind of renewable energy, has many advantages, such as lower pollution and low exploitation costs. It is also more renewable and predictable compared with fossil energy. Therefore, many countries, including China, are increasing their investment in developing tidal energy [2]. The development and utilization of tidal energy begins with transforming the kinetic energy of the tide into electric energy using an energy transducer. Energy transducer devices can be divided into two types: Impeller and non-impeller. The turbine is widely used as an impeller conversion device for ocean tidal energy. Turbines can be divided into the horizontal axis and vertical axis according to the relationship between the flow direction and the spindle. For horizontal axis turbines, the spindle is parallel to the flow direction. In contrast, the spindle of a vertical axis turbine is perpendicular to the flow direction. The vertical axis turbine is not affected by flow direction and has the advantages of a simple blade structure, low working speed ratio and low noise. Moreover, the power systems can be arranged above the water. As a result, vertical axis tidal turbines have recently attracted particular attention [3].
At present, prediction methods for the hydrodynamic performance of vertical axis tidal turbines can broadly be divided into the following four types: Stream-tube model methods [4][5][6], vortex model methods [7][8][9], computational fluid dynamic (CFD) methods [2,10,11] and model test methods [12][13][14][15]. Stream-tube model methods are based on the momentum theorem and have difficulties predicting the flow field characteristics of turbines when the speed ratio increases to a point past which the momentum equation may diverge, leading to no solution. The classic vortex model uses a bound vortex filament to replace the rotor blade and changes its strength as a function of azimuthal position. However, the model cannot fulfill the Kutta condition, resulting in an inaccurate calculation of the forces on blades. The free vortex method combined with finite element analysis (FEVDTM) of the flow surrounding the blades was developed by Ponta [8]. This improves the accuracy of vortex models, but the computational cost is expensive. In recent years, the application of CFD methods to analyze turbine performance has become more popular, but too many factors can affect the results of CFD, thus the results need to be verified by experimental tests. Further, the computational cost of CFD is very high and becomes increasingly prohibitive when used for a full-scale turbine or multiple turbines. Compared with actual conditions in engineering applications, the computational results from a small-scale model may have scaling effects, which are difficult to evaluate with accuracy. Model tests are comparatively more reliable but take a longer time and greater effort and thus have a higher cost.
In this paper, a numerical model is developed based on a combination of the boundary element method (BEM) and vortex theory to predict the hydrodynamic performance of vertical axis turbines. The dimension-reducing technique of the boundary element method is used to reduce the computation cost and vortex theory is applied to wake vortices, which are simulated by a system of discrete vortex cores. The Kutta condition is implemented by forcing zero pressure jumping at the trailing edge between the upper and lower surfaces. The time-stepping method is used to solve the unknown strength of sources and vortices and obtain the velocity distribution along the blades at any time. A model test is also carried out to verify our numerical methods.

Methods
The blade of a vertical axis turbine is a slender body; it can be treated simply as the two-dimensional cross section of the blade. Here, we consider unsteady, two-dimensional, inviscid, incompressible flows around a multi-blade vertical turbine, as shown in Figure 1. A moving coordinate system (OXY) fixed on one of the blades is established, as shown in Figure 2. The motion of the blade at time t can be represented by two components, the rotating angular velocity Ω around the origin o and the translational velocity U of the origin o. The azimuthal angle θ indicates the position of the blade in the track circle, θ = θ 0 + ωt, in which ω is the angular velocity around the rotation center O and θ 0 is the initial azimuthal angle of the blade. At present, prediction methods for the hydrodynamic performance of vertical axis tidal turbines can broadly be divided into the following four types: Stream-tube model methods [4][5][6], vortex model methods [7][8][9], computational fluid dynamic (CFD) methods [2,10,11] and model test methods [12][13][14][15]. Stream-tube model methods are based on the momentum theorem and have difficulties predicting the flow field characteristics of turbines when the speed ratio increases to a point past which the momentum equation may diverge, leading to no solution. The classic vortex model uses a bound vortex filament to replace the rotor blade and changes its strength as a function of azimuthal position. However, the model cannot fulfill the Kutta condition, resulting in an inaccurate calculation of the forces on blades. The free vortex method combined with finite element analysis (FEVDTM) of the flow surrounding the blades was developed by Ponta [8]. This improves the accuracy of vortex models, but the computational cost is expensive. In recent years, the application of CFD methods to analyze turbine performance has become more popular, but too many factors can affect the results of CFD, thus the results need to be verified by experimental tests. Further, the computational cost of CFD is very high and becomes increasingly prohibitive when used for a full-scale turbine or multiple turbines. Compared with actual conditions in engineering applications, the computational results from a small-scale model may have scaling effects, which are difficult to evaluate with accuracy. Model tests are comparatively more reliable but take a longer time and greater effort and thus have a higher cost.
In this paper, a numerical model is developed based on a combination of the boundary element method (BEM) and vortex theory to predict the hydrodynamic performance of vertical axis turbines. The dimension-reducing technique of the boundary element method is used to reduce the computation cost and vortex theory is applied to wake vortices, which are simulated by a system of discrete vortex cores. The Kutta condition is implemented by forcing zero pressure jumping at the trailing edge between the upper and lower surfaces. The time-stepping method is used to solve the unknown strength of sources and vortices and obtain the velocity distribution along the blades at any time. A model test is also carried out to verify our numerical methods.

Methods
The blade of a vertical axis turbine is a slender body; it can be treated simply as the two-dimensional cross section of the blade. Here, we consider unsteady, two-dimensional, inviscid, incompressible flows around a multi-blade vertical turbine, as shown in Figure 1.   The velocity of uniform incoming flow is V A , which at infinity is Φ(p, t) is the velocity potential at point p and time t, which meets the Laplace equation in the fluid domain τ e : ∇ 2 Φ(p, t) = 0 (p ∈ τ e ) (2) and the non-penetration boundary condition on the blade surface: where n b is the normal of the blade surface S b ; and r is the vector from the origin o to any surface point in the coordinate system xoy. The velocity of uniform incoming flow is A V , which at infinity is ( ) Φ τ ∇ = ∈ e p t p (2) and the non-penetration boundary condition on the blade surface: where nb is the normal of the blade surface Sb; and r is the vector from the origin o to any surface point in the coordinate system xoy.
At infinity, ( ) Introducing φ as the perturbation velocity potential, Φ can be decomposed into two parts: Assuming zero initial perturbation, φ should satisfy the following equation and boundary conditions: At infinity, Φ(p, t) satisfies the velocity potential for uniform incoming flow, Introducing φ as the perturbation velocity potential, Φ can be decomposed into two parts: Assuming zero initial perturbation, φ should satisfy the following equation and boundary conditions: The boundary element method is used to solve the time-discrete Laplace equations in combination with the Neumann boundary conditions as described by Equations (6)- (9). As shown in Figure  shed from the trailing edge and move at the local particle velocity, as shown in Figure 3. As a result, at an arbitrary point p in the fluid field, the induced velocity can be expressed as: , which are shed from the trailing edge and move at the local particle velocity, as shown in Figure 3. As a result, at an arbitrary point p in the fluid field, the induced velocity can be expressed as: To compute Equation (10)  To compute Equation (10), we divide each blade surface S b into N elements at which the sources σ j (j = 1, N) are distributed, and surface S f into M elements with a linear distribution. Assuming vorticity intensity γ(0) = γ f at the leading edge and γ(l f ) = 0 at the end of the arced centerline, the vorticity intensity on the surface S f can then be described as: at which l f = M ∑ j=1 S f j ds is the length of the arced centerline. The total intensity of the vortices placed on S f is: With the movement of each blade after a time step ∆t, a point vortex γ w is generated at the trailing edge. As a result, a discrete vortex street is formed.
According to the surface boundary condition in Equation (7), Equation (10) for each surface element i can be discretized as where Z is the number of blades, n i is the unit normal vector of the ith element on surface S b and k represents the number of time steps at current time t k (t k = k × ∆t).
In Equation (13) First, the Kutta condition of equal pressure is used on each blade to obtain Z equations, which can be expressed as: That is, the pressure difference between the top and bottom surface of the airfoil trailing edge (T.E.) is zero. Applying the unsteady Bernoulli equation at an arbitrary point in the fluid field, the dynamic pressure P can be expressed as: Applying Equation (15) to Equation (14) produces: where p u is the pressure on the upper surface of the airfoil trailing edge, p d is the pressure on the lower surface of the airfoil trailing edge, φ u is the perturbation velocity potential on the upper surface of the airfoil trailing edge, φ d is the perturbation velocity potential on the lower surface of the airfoil trailing edge, r u is the vector of the upper surface of the airfoil trailing edge in the local body-fixed coordinate system and r d is the vector of the lower surface of the airfoil trailing edge in the local body-fixed coordinate system. After obtaining the perturbation velocity on the blade surface unit numerically from Equation (10), the pressure distribution in each element can be calculated using the unsteady Bernoulli equation. Then, the hydrodynamic force on the airfoil blade can be obtained by numerical integration. The time derivative of the perturbation potential can be obtained by the finite difference method: At any time t, the perturbation velocity potential φ(t) is determined by the line integral of the perturbation velocity V(t). Theoretically, the starting point of the integral should be infinity. However, in actual calculation the starting point is 10 times the distance from the rotation center O of the water turbine. Integrating from the starting point to the rotation center O, then from O to the lower point q 1 of each trailing edge, the perturbation velocity potential of any element center q in the whole airfoil can be derived: where φ 0 = φ(q 1 ) is the perturbation velocity potential at q 1 . Second, the Kelvin condition is applied at each blade to obtain Z equations as follows: where Γ where β is a constant coefficient, β = 0.4 ∼ 0.6. The blade surface pressure coefficient is defined as: where ρ is the density of water. The resulting force on the blade is: The moment on the blade reference point o is: The tangential force coefficient of the turbine blade is defined as: and the normal force coefficient of the turbine blade as: where f t is the turbine blade tangential force, f n is the turbine blade normal force, C is the turbine blade chord and b is the turbine blade length.

Model Validation against Past Studies
Previously, many researchers have studied the hydrodynamic performance of vertical axis turbines. A representative study is Strickland's model test study [15]. The parameters for one of the turbine models (Turbine A) in his work are shown in Table 1. The normal force, tangential force and wake evolution were also recorded in his experiment. Strickland's test is often used as a benchmark for the validation of numerical models. In addition, there are a number of representative modeling studies, such as the FEVDTM proposed by Ponta [8] and the CFD method proposed by Li [2]. This paper first compares the predictions from different methods with the turbine model given in Table 1. In the numerical computation, the number of elements is chosen as N = 80, M = 40 and the time step is ∆t = T/360, where T is the rotation period of the turbine, T = 2π/ω; further fine mesh and smaller time step present very similar result. Figures 4 and 5 compare the predicted central path of the turbine wake vortices from different studies for blade number Z = 1 and speed ratio λ = ωR/V A = 5. In Figure 4, the light streak shows the wake's conformation obtained from Strickland's model test using the dye injection visualization technique [15] and the line with marked circles is the prediction result based on Ponda's FEVDTM model [8]. Figure 5 demonstrates the result predicted by the method proposed in this paper. Comparing the results in Figures 4 and 5, the turbine vortex wakes predicted by the proposed model are in good agreement with the experiment and the FEVDTM model.
In the numerical computation, the number of elements is chosen as N = 80, M = 40 and the time step is / 360 t T Δ = , where T is the rotation period of the turbine, T = 2π/ω; further fine mesh and smaller time step present very similar result. Figures 4 and 5 compare the predicted central path of the turbine wake vortices from different studies for blade number Z = 1 and speed ratio λ = ωR/VA = 5. In Figure 4, the light streak shows the wake's conformation obtained from Strickland's model test using the dye injection visualization technique [15] and the line with marked circles is the prediction result based on Ponda's FEVDTM model [8]. Figure 5 demonstrates the result predicted by the method proposed in this paper. Comparing the results in Figures 4 and 5, the turbine vortex wakes predicted by the proposed model are in good agreement with the experiment and the FEVDTM model.  Next, a quantitative comparison is carried out. The normal and tangential force coefficients against the azimuthal angle θ (as shown in Figure 2) of one blade of the turbine are calculated for the cases of Z = 2 and λ = 2.5, 5, 7.5. The results are presented in Figures 6-8. The results from the experiment and the CFD model in the literature are also shown for comparison purposes. The quantitative results predicted by the present method agree with those predicted by Li's CFD method [2] and Strickland's experiment [15]. Compared with the experimental test, the measured force at the peak and trough based on numerical methods is somewhat different. Overall, both the CFD and our method predict the normal force more accurately than the tangential force when compared with the experimental results. This may be attributed to either the tangled wake field of the turbines or ignoring the viscous effects of the fluid. However, the tangential force is one order lower than the normal force, so it has less impact on the performance analysis.
In the case of Z = 2, the total number of elements is 240, the time step is / 360 t T Δ = , the motion time is 4T including 1440 time steps and the computation time is only 15 min using a one-core 2.5 GHz CPU. The calculation time and computer memory needed by the proposed method are far less than the CFD method. Based on the above analyses and comparisons, it appears that the proposed method can satisfy the engineering requirements of fast prediction of turbine performance and hydrodynamic interference between multiple turbines in tidal power plants.  Next, a quantitative comparison is carried out. The normal and tangential force coefficients against the azimuthal angle θ (as shown in Figure 2) of one blade of the turbine are calculated for the cases of Z = 2 and λ = 2.5, 5, 7.5. The results are presented in Figures 6-8. The results from the experiment and the CFD model in the literature are also shown for comparison purposes. The quantitative results predicted by the present method agree with those predicted by Li's CFD method [2] and Strickland's experiment [15]. Compared with the experimental test, the measured force at the peak and trough based on numerical methods is somewhat different. Overall, both the CFD and our method predict the normal force more accurately than the tangential force when compared with the experimental results. This may be attributed to either the tangled wake field of the turbines or ignoring the viscous effects of the fluid. However, the tangential force is one order lower than the normal force, so it has less impact on the performance analysis.
In the case of Z = 2, the total number of elements is 240, the time step is ∆t = T/360, the motion time is 4T including 1440 time steps and the computation time is only 15 min using a one-core 2.5 GHz CPU. The calculation time and computer memory needed by the proposed method are far less than the CFD method. Based on the above analyses and comparisons, it appears that the proposed method can satisfy the engineering requirements of fast prediction of turbine performance and hydrodynamic interference between multiple turbines in tidal power plants.
GHz CPU. The calculation time and computer memory needed by the proposed method are far less than the CFD method. Based on the above analyses and comparisons, it appears that the proposed method can satisfy the engineering requirements of fast prediction of turbine performance and hydrodynamic interference between multiple turbines in tidal power plants.

Comparison with a Laboratory Test
To study the hydrodynamic performance of the turbine further and validate the proposed method, a model test was conducted in a towing tank. The towing tank was 130 m long, 6 m wide and 4 m deep. The maximum speed of the trailer was 6.5 m/s. A vertical-axis turbine model was fixed on the trailer and then driven by it. The turbine was driven by a motor through the reduction gear box. Due to limited test conditions, only the torque of the hydraulic turbine model was measured. Figure 9 shows the apparatus used to measure the torque of the vertical-axis turbine:

Comparison with a Laboratory Test
To study the hydrodynamic performance of the turbine further and validate the proposed method, a model test was conducted in a towing tank. The towing tank was 130 m long, 6 m wide and 4 m deep. The maximum speed of the trailer was 6.5 m/s. A vertical-axis turbine model was fixed on the trailer and then driven by it. The turbine was driven by a motor through the reduction gear box. Due to limited test conditions, only the torque of the hydraulic turbine model was measured. Figure 9 shows the apparatus used to measure the torque of the vertical-axis turbine: Figure 9A shows the test layout, Figure 9B the five-blade vertical-axis turbine model, Figure 9C the torque meter and Figure 9D a 32-channel dynamic data acquisition instrument. A torque meter was arranged between the turbine and the driving motor to measure the torque produced by the turbine at a given inflow speed and speed in real time. To avoid the influence of the free surface, the turbine was immersed below 0.4 m. The accuracy of the torque meter is 0.3%. The center of the torque meter is a strain shaft. The upper end of the strain shaft is the driving motor and the lower end is the turbine. When the turbine rotates and causes the strain shaft to be slightly deformed by torsion, the resistance value of the strain meter pasted on the strain shaft changes correspondingly. The change in strain resistance can be transformed into a change in voltage signal which is output to the dynamic data acquisition instrument, and the transient value of the hydraulic turbine torque can be obtained after the transformation.
where D is the turbine diameter. The torque coefficient against the azimuthal angle θ after the third revolution for two different turbine models and speed ratios predicted by the proposed method are shown in Figures 10-13. Figures 10 and 11 show that the torque coefficient curve for a turbine with three blades has three peaks in one cycle, while Figures 12 and 13 show that a five-blade turbine produces five peaks. The number of peaks in one cycle is related to the number of blades. Overall, the computed results are in good agreement with those obtained from experiments, which indicates that the model can capture the torque performance of the turbine well. However, a difference in the peaks and valleys is observed, which may be attributed to the following reasons. First, the proposed method is based on potential flow theory, in which the viscous effect is ignored. Second, it is a two-dimensional model and the experiment is carried out under three-dimensional conditions. These effects will be investigated in greater depth in future studies.   The model parameters used in the experiment are listed in Table 2. Two models were tested: One with three blades and the other with five blades. To match operating conditions of the vertical-axis tidal current turbine power station in Guanshan Island, Zhoushan, China, two typical speed ratios were selected for each turbine model: 1.6 and 2.2 for the three-blade model and 1.65 and 2.23 for the five-blade model. The torque test results Q can be transformed into a torque coefficient where D is the turbine diameter. The torque coefficient against the azimuthal angle θ after the third revolution for two different turbine models and speed ratios predicted by the proposed method are shown in Figures 10-13. Figures 10 and 11 show that the torque coefficient curve for a turbine with three blades has three peaks in one cycle, while Figures 12 and 13 show that a five-blade turbine produces five peaks. The number of peaks in one cycle is related to the number of blades. Overall, the computed results are in good agreement with those obtained from experiments, which indicates that the model can capture the torque performance of the turbine well. However, a difference in the peaks and valleys is observed, which may be attributed to the following reasons. First, the proposed method is based on potential flow theory, in which the viscous effect is ignored. Second, it is a two-dimensional model and the experiment is carried out under three-dimensional conditions. These effects will be investigated in greater depth in future studies.

Conclusions
In this paper, the flow induced by the motion of the vertical-axis turbine has been simplified as a two-dimensional problem and solved using the unsteady boundary element method. A code for the proposed method was developed specifically for performance analysis of vertical-axis turbines. In comparison with previous studies, the proposed model is simpler and more efficient. It can capture our experimental results well and satisfy the engineering application requirements for performance prediction of vertical-axis turbines. The results imply that the proposed method can be adopted as a tool for analyzing the hydrodynamic performance of vertical-axis tidal turbines. The code developed provides a useful design tool for optimizing the airfoil profile and arrangement of tidal power devices. In this study, we ignored the viscous effect of the fluid and an uncertainty analysis was not carried out. These will be incorporated in our future research and may lead to a better prediction. Another future research direction is to develop the method to solve three-dimensional problems.