A Method for Predicting Unsteady Characteristics of Horizontal Axis Tidal Current Turbine Blades

Abstract

The present invention discloses a method for predicting the unsteady characteristics of a horizontal axis tidal current turbine blade, comprising: 1) obtaining blade profile parameters, including chord length distribution, pitch angle distribution, and thickness distribution; 2) dividing the blade into multiple blade elements and calculating the static and dynamic parameters of each blade element; 3) on the basis of a uniform incoming flow, superimposing a shear flow synthesized by the power law, a wave synthesized by the Stokes wave theory, and a turbulence synthesized by the turbulence spectrum method to obtain a synthesized incoming flow; 4) coupling the L-B dynamic stall model, the synthesized incoming flow, and the steady-state blade element momentum theory model to form a transient blade element momentum theory model, and through iterative calculation of the unsteady forces and moments received by each blade element, and then integrating along the radial direction to obtain the unsteady characteristics received by the blade. The prediction method proposed by the present invention can comprehensively consider the actual incoming flow, blade thickness, dynamic stall, and three-dimensional rotation enhancement effect, and the prediction result can be closer to the actual situation.

Description

Technical Field

[0001] The present invention relates to the field of tidal current power generation, and particularly to a method for predicting the unsteady characteristics of a horizontal-axis tidal current turbine blade. Background Art

[0002] Currently, most countries mainly meet the growing energy demand by consuming fossil fuels. The resulting greenhouse gas emissions and global warming problems have attracted widespread attention worldwide, and it has gradually been realized that it is necessary to incorporate renewable energy into the alternative energy policy for fossil fuels. Tidal current energy, with its advantages of high power density, strong predictability, and small impact on the power grid, has become the most commercially promising marine new energy technology after hydropower, wind power, and photovoltaic power generation. The horizontal-axis type occupies approximately 76% of the global tidal current power generation devices with its high efficiency and advantages that have been effectively verified in the wind power industry, and this share is still increasing.

[0003] However, there are still many obstacles to overcome to realize the commercial potential in this field. The complexity of the service environment affecting the hydrodynamic performance of horizontal-axis tidal current turbines is mainly reflected in typical unsteady phenomena such as wave motion in the impeller inlet flow field, high turbulence intensity levels caused by irregular seabeds, and non-uniform oncoming flow. This results in large load fluctuations experienced by the rotor during operation. The blades of a rotor may experience stall delay, load lag, and dynamic stall, which may lead to dynamic and fatigue failures. Therefore, accurately predicting unsteady loads is crucial.

[0004] Currently, the problems existing in the methods for predicting the unsteady loads of horizontal-axis tidal current turbine blades at home and abroad are as follows: First, the simulation of the service environment of the turbine is not realistic enough. Most only consider one or two of the effects of waves, shear flow, and turbulence, which is not comprehensive enough. Second, the impeller is simplified to have a constant thickness from the root to the tip, which has a large deviation from the actual impeller. The relative thickness of the actual impeller gradually thins from the root to the tip. Therefore, the predicted results have a certain deviation from the actual situation. Summary of the Invention

[0005] In view of the deficiencies of the prior art, in order to predict the unsteady characteristics of the blades of a non-uniform-thickness horizontal-axis tidal current turbine operating in a real marine environment, the present invention proposes a method for predicting the unsteady characteristics of a non-uniform-thickness horizontal-axis tidal current turbine blade. This method comprehensively considers the effects of waves, shear flow, and turbulence, as well as dynamic stall and three-dimensional rotation enhancement effects, providing a reference for the strength check during the design process of horizontal-axis tidal current turbines.

[0006] The object of the present invention is achieved through the following technical solutions:

[0007] A method for predicting the unsteady characteristics of a horizontal-axis tidal current turbine blade, comprising the following steps:

[0008] Step 1: Obtain the blade profile parameters, including chord length distribution, pitch angle distribution, and thickness distribution;

[0009] Step 2: Divide the blade into multiple blade elements and calculate the static and dynamic parameters of each blade element;

[0010] Step 3: On the basis of a uniform oncoming flow, superimpose a shear flow synthesized by the power law, a wave synthesized by the Stokes wave theory, and a turbulence synthesized by the turbulence spectrum method to obtain a synthesized oncoming flow;

[0011] Step 4: Couple the L - B dynamic stall model, the synthesized oncoming flow, and the steady - state blade element momentum theory model to form a transient blade element momentum theory model, and obtain the unsteady characteristics of the blade through iterative calculations;

[0012] The specific steps for obtaining the unsteady characteristics of the blade through iterative calculations are as follows:

[0013] (1) Calculate the angle of attack of the blade elements at different radii through the synthesized oncoming flow; calculate the lift coefficient and drag coefficient of the blade elements through the L - B dynamic stall model;

[0014] (2) Consider the three - dimensional rotation enhancement effect to obtain the actual lift coefficient and drag coefficient of the blade elements;

[0015] (3) Calculate the forces and torques received by each blade element through the actual lift coefficient and drag coefficient, and then integrate along the span direction to obtain the unsteady forces and torques received by the blade.

[0016] Furthermore, the static parameters include the lift coefficient CL, drag coefficient CD, normal force coefficient CN, and static stall onset angle α of the blade element at different angles of attack ss , zero - lift drag coefficient CD0, and the slope CN of the linear region of the normal force coefficient α ; the dynamic parameters include the critical dynamic stall onset angle α ds0 , constant time - delay constant T α and the time T of vortex transfer vL .

[0017] Furthermore, the synthesized shear inflow preferably adopts the 1 / 7 power law

[0018]

[0019] where U xshear is the synthesized shear inflow, U hubis the average velocity at the impeller hub, z is the distance in the z - direction of the Cartesian coordinate system established at the hub center, with the coordinate system pointing positive towards the horizontal plane, h is the installation water depth of the hub, γ = 1 / 7; H is the water depth at the location of the water turbine.

[0020] Further, the synthesis of waves adopts the Stokes second - order wave theory

[0021]

[0022]

[0023] where H S is the wave height, g is the acceleration due to gravity, K is the wave number, ω a is the angular velocity wave frequency,

[0024]

[0025] where, T a is the wave period; H S is the wave height, g is the acceleration due to gravity, J is the wave number, ω a is the angular velocity wave frequency, u xwave represents the flow - direction velocity caused by the waves; u zwave represents the vertical velocity caused by the waves.

[0026] Further, the synthesis of turbulence adopts the von Kármán spectrum method, where the flow - direction is the x - direction, and its velocity spectrum S x is

[0027]

[0028] where, L x is the length scale in the x - direction, σ x is the standard deviation in the x - direction, σ x = I x U x where, I x is the flow - direction turbulence intensity; U x is the flow - direction average velocity; n represents the dimensionless frequency, n = L x f t / U x f t is the turbulence frequency component;

[0029] The velocity spectrum S in the y - direction y is

[0030]

[0031] where, σ y = R t σx , L y is the length scale in the y direction; L y = R t L x , R t the anisotropy ratio, where R t = 1; σ y is the standard deviation in the y direction;

[0032] The velocity spectrum S in the z direction z = S y ;

[0033] Then the turbulent synthetic velocity in each direction is

[0034]

[0035] where i represents the three directions of x, y, and z, Δf t is the frequency interval size, N is the number of f t , Φ j is the random phase angle, and its range is from 0 to 2π.

[0036] Furthermore, the static parameters of the blade element are calculated using Xfoil, and the calculation steps are as follows:

[0037] (1) Use Xfoil to calculate the lift coefficient CL and drag coefficient CD of the blade element at the local angle of attack; the local angle of attack is selected from -20° to 30°, and the calculation is performed every 1°;

[0038] (2) Calculate the lift coefficient and drag coefficient at angles of attack from -180° to 180° through Viterna extrapolation;

[0039] (3) Convert or extract other static parameters through the lift coefficient CL and drag coefficient CD.

[0040] Furthermore, the relationship between the critical dynamic stall onset angle of attack α ds0 and the constant time delay constant T α is:

[0041] α ds = α ds0 + T α * r

[0042] where r is the reduced frequency; α ds represents the dynamic stall onset angle of attack;

[0043] The dynamic parameters of the blade element are calculated using Fluent, and the calculation steps are as follows:

[0044] (1) Calculate the dynamic stall onset angle of attack α of the blade element at different reduced frequencies ds ;

[0045] (2) By linearly fitting the curves of different reduced frequencies and the onset angle of dynamic stall α ds , the critical onset angle of dynamic stall α ds0 and the constant time delay constant T α are obtained;

[0046] (3) Calculate the time T of vortex transfer through the following formula vL

[0047]

[0048] where V is the inflow velocity, c is the chord length of the blade element, t f is the time of dynamic stall vortex shedding, and t0 is the start time of dynamic stall. These two times are obtained by extracting the dynamic stall characteristics of the airfoil calculated by Fluent.

[0049] The beneficial effects of the present invention are as follows:

[0050] The unsteady characteristic prediction method for the blades of a horizontal-axis tidal current turbine proposed by the present invention comprehensively considers the effects of waves, shear flow, and turbulence, as well as dynamic stall and three-dimensional rotation enhancement effects. The prediction results are more consistent with the actual situation and can provide a reference for the strength check during the design process of horizontal-axis tidal current turbines. Description of the Drawings

[0051] The present invention will be further described in detail below in conjunction with the drawings and specific embodiments:

[0052] Figure 1 is the flow chart of the present invention;

[0053] Figure 2 is the flow chart of the embodiment of the present invention;

[0054] Figure 3 is the airfoil parameter diagram of a certain horizontal-axis tidal current turbine blade;

[0055] Figure 4 is the lift and drag coefficient distribution of 150 blade elements of the blade; among them, Figure (a) is the lift coefficient of 150 blade elements; Figure (b) is the drag coefficient of 150 blade elements;

[0056] Figure 5 is the schematic diagram of the airfoil doing sinusoidal oscillatory motion;

[0057] Figure 6 is the verification of the steady-state BEM model;

[0058] Figure 7 is the comparison diagram of the unsteady power coefficient CP model prediction and CFD calculation of the impeller;

[0059] Figure 8 It is a comparison diagram of the predicted value and CFD calculation of the unsteady thrust coefficient CT of the impeller. Specific implementation mode

[0060] The present invention will be described in detail below with reference to the drawings and preferred embodiments. The purpose and effect of the present invention will become more apparent. It should be understood that the specific embodiments described herein are only used to explain the present invention and are not used to limit the present invention.

[0061] As Figure 1 and 2 shown, the method for predicting the unsteady characteristics of the blades of a horizontal-axis tidal current turbine according to the present invention includes the following steps:

[0062] Step 1: Obtain the blade profile parameters, including chord length distribution, pitch angle distribution, and thickness distribution, as Figure 3 shown;

[0063] Step 2: Divide the blade into multiple blade elements, and calculate the static parameters and dynamic parameters of each blade element;

[0064] In this embodiment, the blade is divided into 150 blade elements, and 15 blade elements are selected at equal intervals as the main blade elements. The static parameters and dynamic parameters of the main blade elements are calculated by Xfoil and Fluent, and then the other blade elements that have not been calculated are obtained by interpolation.

[0065] Among them, the static parameters include the lift coefficient CL, drag coefficient CD, normal force coefficient CN, static stall starting angle of attack α ss , zero-lift drag coefficient CD0, and the slope CN of the linear region of the normal force coefficient α . The dynamic parameters include the critical dynamic stall starting angle of attack α ds0 , constant time delay constant T α , and the time T of vortex transfer vL . When calculating the static and dynamic parameters, only 15 main blade elements are calculated by Xfoil and Fluent. When obtaining the static parameters of each blade element, Xfoil needs to calculate the static lift coefficient CL and drag coefficient CD of the airfoil at the local angle of attack. The local angle of attack is preferably selected from -20° to 30°, and calculated every 1°. After obtaining the lift coefficient and drag coefficient at the local angle of attack, the lift coefficient and drag coefficient at the angle of attack from -180° to 180° are calculated by the Viterna extrapolation method, as Figure 4 shown. Other static parameters can be converted or extracted from the lift coefficient CL and drag coefficient CD. Using Fluent to calculate the dynamic parameters mainly controls the blade element to perform a sinusoidal oscillating motion to induce dynamic stall. The motion schematic diagram is as Figure 5As shown. Then, obtain the variation laws of the dynamic lift coefficient CL, drag coefficient CD, and moment coefficient CM of the blade element with the angle of attack, and extract the angle of attack α at the onset of dynamic stall and the time T of vortex transfer through dynamic characteristics ds and the time T of vortex transfer vL By calculating the dynamic characteristics at different reduced frequencies r, obtain the angle of attack α at the onset of dynamic stall and the time T of vortex transfer at different reduced frequencies ds and the time T of vortex transfer vL , through the relationship α ds = α ds0 + T α * r to obtain the critical angle of attack α at the onset of dynamic stall ds0 and the constant time delay constant T α . Generally, it is sufficient to calculate three to five points in the reduced frequency range of 0.01 - 0.05. Of course, the more calculations are performed, the more accurate the fitting will be. The blade elements that are not calculated are obtained by interpolation.

[0066] Step 3: On the basis of a uniform oncoming flow, superimpose a shear flow synthesized by the power law, a wave synthesized using Stokes wave theory, and a turbulence synthesized using the turbulence spectrum method to obtain a synthesized oncoming flow.

[0067] In this embodiment, the synthesized shear oncoming flow preferably adopts the 1 / 7 power law

[0068]

[0069] where U xshear is the synthesized shear oncoming flow, U hub is the average velocity at the impeller hub, z is the distance in the z - direction of the Cartesian coordinate system established at the hub center, with the coordinate system pointing upward in the horizontal plane as positive, h is the installation water depth of the hub, γ = 1 / 7; H is the water depth at the location of the water turbine.

[0070] The synthesis of the wave adopts Stokes second - order wave theory

[0071]

[0072]

[0073] where H S is the wave height, g is the acceleration due to gravity, K is the wave number, ω a is the angular velocity wave frequency, u xwave represents the flow - direction velocity caused by the wave; u zwave represents the vertical velocity caused by the wave.

[0074] The synthesis of the turbulence adopts the von Karman spectrum method, where the flow - direction is the x - direction, and its velocity spectrum S x is

[0075]

[0076] Among them, L x is the length scale in the x direction, and σ x is the standard deviation in the x direction. σ x = I x U x , where I x is the streamwise turbulence intensity; U x is the streamwise mean velocity; n represents the dimensionless frequency, and n = L x f t / U x , f t is the turbulence frequency component;

[0077] The velocity spectrum S y in the y direction is S y

[0078] Among them, σ y = R t σ x , L y is the length scale in the y direction; L y = R t L x , R t is the anisotropy ratio. Here, R t = 1; σ y is the standard deviation in the y direction.

[0079] The velocity spectrum S z in the z direction is S y .

[0080] Then the synthetic turbulence velocity in each direction is

[0081]

[0082] where i represents the three directions of x, y, and z, Δf t is the frequency interval, N is the number of f t , and Φ j is the random phase angle, and its range is from 0 to 2π;

[0083] Step 4: Couple the L - B dynamic stall model, the synthetic inflow, and the steady - state blade element momentum theory model to form a transient blade element momentum theory model, and obtain the unsteady characteristics of the blade through iterative calculations;

[0084] The specific steps for obtaining the unsteady characteristics of the blade through iterative calculations are as follows:

[0085] (1) Calculate the angle of attack of the blade element at different radii through the synthetic inflow calculation; calculate the lift coefficient and drag coefficient of the blade element through the L-B dynamic stall model;

[0086] (2) Considering the three-dimensional rotation enhancement effect, obtain the actual lift coefficient and drag coefficient of the blade element;

[0087] (3) Calculate the force and moment received by each blade element through the actual lift coefficient and drag coefficient, and then integrate along the span direction to obtain the unsteady force and moment received by the blade.

[0088] The method of the present invention is named HATT-ULPM. As Figure 6 shown, by comparing the power coefficient CP and thrust coefficient CT calculated by the steady blade element momentum theory (BEM) with the CFD simulation results, verify the correctness of the steady BEM model applied by this method. As Figure 7 and Figure 8 shown, under the same unsteady inflow conditions, compared with CFD, the CP and CT calculated by applying HATT-ULPM have a higher fitting degree with the CFD simulation results, thus verifying the accuracy of the unsteady characteristics of the blade predicted by this method.

[0089] Those of ordinary skill in the art can understand that the above are only preferred examples of the invention and are not used to limit the invention. Although the invention has been described in detail with reference to the foregoing examples, for those skilled in the art, they can still modify the technical solutions recorded in the foregoing examples, or perform equivalent replacements for some of the technical features. Any modifications, equivalent replacements, etc. made within the spirit and principle of the invention shall be included within the protection scope of the invention.

Claims

1. A method for predicting the unsteady characteristics of blades of a horizontal-axis tidal current turbine, characterized in that, It includes the following steps: Step 1: Obtain the blade profile parameters, including chord length distribution, pitch angle distribution, and thickness distribution; Step 2: Divide the blade into multiple blade elements and calculate the static and dynamic parameters of each blade element; Step 3: On the basis of a uniform oncoming flow, superimpose the shear flow synthesized by the power-law, the wave synthesized by the Stokes wave theory, and the turbulence synthesized by the turbulence spectrum method to obtain the synthesized oncoming flow; Step 4: Couple the L-B dynamic stall model, the synthesized oncoming flow, and the steady-state blade element momentum theory model to form a transient blade element momentum theory model, and obtain the unsteady characteristics of the blade through iterative calculations; The specific steps for obtaining the unsteady characteristics of the blade through iterative calculations are as follows: (1) Calculate the angle of attack of the blade elements at different radii of the blade through the synthesized oncoming flow; calculate the lift coefficient and drag coefficient of the blade elements obtained through the L-B dynamic stall model; (2) Consider the three-dimensional rotation enhancement effect to obtain the actual lift coefficient and drag coefficient of the blade elements; (3) Calculate the forces and torques received by each blade element through the actual lift coefficient and drag coefficient, and then integrate along the span direction to obtain the unsteady forces and torques received by the blade; The synthesized shear oncoming flow adopts the 1 / 7 power-law Among them, U xshear is the synthesized shear inflow, and U hub is the average velocity at the impeller hub. z is the distance in the z-direction of the Cartesian coordinate system established at the hub center, with the coordinate system pointing to the horizontal plane being positive, h is the installation water depth of the hub, γ = 1 / 7; H is the water depth at the location of the water turbine; The synthesis of the wave adopts the Stokes second-order wave theory where T a is the wave period; H S is the wave height, g is the acceleration due to gravity, K is the wave number, ω a is the angular velocity wave frequency, u xwave represents the flow velocity induced by the wave; u zwave represents the vertical velocity induced by the wave; The synthesis of turbulence adopts the von Karman spectrum method, where the flow direction is the x direction, and its velocity spectrum S x is Among them, L x is the length scale in the x direction, σ x is the standard deviation in the x direction, σ x = I x U x where I x is the streamwise turbulence intensity; U x is the streamwise mean velocity; n represents the dimensionless frequency, n = L x f t / U x where f t is the turbulence frequency component; Velocity spectrum S in the y direction y is where σ y = R t σ x , L y is the length scale in the y direction; L y = R t L x , R t is the anisotropy ratio, where R t = 1; σ y is the standard deviation in the y direction; Velocity spectrum S in the z direction z = S y ; Then the synthesized turbulence velocity in each direction is where i represents the three directions of x, y, and z, and Δf t is the frequency interval size, N is the number of f t , and Φ j is the random phase angle, and its range is from 0 to 2π.

2. The unsteady characteristic prediction method for the blades of a horizontal-axis tidal current turbine according to claim 1, wherein: The static parameters include the lift coefficient CL, drag coefficient CD, normal force coefficient CN of the blade element at different angles of attack, and the static stall onset angle of attack α ss , zero-lift drag coefficient CD0, and the slope CN of the linear region of the normal force coefficient α ; the dynamic parameters include the critical dynamic stall onset angle of attack α ds0 , constant time delay constant T α and the time T for vortex transfer vL .

3. The unsteady characteristic prediction method for the blades of a horizontal-axis tidal current turbine according to claim 2, characterized in that: Use Xfoil to calculate the static parameters of the blade element. The calculation steps are as follows: (1) Use Xfoil to calculate the lift coefficient CL and drag coefficient CD of the blade element at the local angle of attack; the local angle of attack is selected from -20° to 30°, and the calculation is performed every 1°; (2) Calculate the lift coefficient and drag coefficient at angles of attack from -180° to 180° through the Viterna extrapolation method; (3) Convert or extract other static parameters through the lift coefficient CL and drag coefficient CD.

4. The unsteady characteristic prediction method for the blades of a horizontal axis tidal current turbine according to claim 2, wherein: Critical dynamic stall onset angle of attack α ds0 and constant time delay constant T α The relationship is as follows: α ds ＝α ds0 +T α *r where r is the reduced frequency; α ds represents the angle of attack at the onset of dynamic stall; Use Fluent to calculate the dynamic parameters of the blade element. The calculation steps are as follows: (1) Calculate the dynamic stall onset angle of attack α of the blade element at different reduced frequencies ds ; (2) By linearly fitting the curves of different reduced frequencies and the dynamic stall onset angle α ds to obtain the critical dynamic stall onset angle α ds0 and the constant time delay constant T α ; (3) Calculate the time T for vortex transfer by the following formula vL where V is the inflow velocity, c is the chord length of the blade element, t f is the time of dynamic stall vortex shedding, and t0 is the start time of dynamic stall. These two times are obtained by extracting the airfoil dynamic stall characteristics calculated by Fluent.

Images