Large wind turbine generator blade aerodynamic performance free vortex wake prediction method

Abstract

The invention discloses a method for predicting the aerodynamic performance free vortex wake of a blade of a large-scale wind turbine generator, and the method comprises the steps: estimating the initial circular rector distribution of each section of the blade through a momentum blade element method or a preset vortex wake method, and determining the space shape of a trailing vortex surface and the initial circular rector distribution; updating the circumferential angle of the blade and the space coordinates of the attached vortex nodes on each section in the inertial system; calculating the induced velocity on the trailing vortex surface at the new moment; rotating and deforming the trailing vortex; updating blade attachment vortex element circular rector distribution, and obtaining a blade aerodynamic characteristic result at the current moment; and along with the update of the inflow wind speed, returning to the step S3 to perform cyclic calculation until a preset maximum time step number is reached. According to the method, the number of iterations of the attached vortex circular rector distribution of each time step can be remarkably reduced, and the calculation efficiency of dynamic aerodynamic characteristic prediction of the wind turbine generator based on a free vortex wake numerical method is improved, so that rapid prediction of dynamic aerodynamic performance of a large wind turbine blade and safety diagnosis under a batch wind condition are facilitated.

Description

Technical Field

[0001] The present invention relates to the technical field of new energy wind power generation, and particularly to a free vortex wake prediction method for the aerodynamic performance of large wind turbine blades. Background Art

[0002] The free vortex wake numerical algorithm combines the inviscid potential flow lifting line / surface theory with the vortex dynamics method. By predicting the evolution law of the wake vortex surface of the wind turbine blade and the induction effect on the rotor cross-section, the dynamic aerodynamic load information of the blade can be obtained, which can take into account both the prediction accuracy and the computational economy.

[0003] Compared with the existing numerical calculation methods, the free vortex wake numerical algorithm has obvious advantages over the momentum blade element method in terms of the prediction accuracy of the vortex induction effect; compared with the Reynolds-averaged Navier-Stokes equation solution method or CFD methods such as large eddy simulation, the free vortex wake algorithm can effectively control the computational time required for engineering evaluation, providing an efficient and accurate evaluation method for the dynamic aerodynamic loads of large wind turbine blades under complex inflow conditions, providing aerodynamic data support for the aeroelastic response analysis of large wind turbine ultra-long flexible blades, and assisting in the prediction of extreme aerodynamic loads of large and extra-large blades and the evaluation of structural aeroelastic stability and operation safety.

[0004] At present, the existing free vortex wake numerical algorithm has many iterative steps, the overall computational efficiency and prediction accuracy of the unsteady aerodynamic performance of large wind turbine blades are poor, and safety problems are likely to occur due to inaccurate prediction of blade dynamic loads. Summary of the Invention

[0005] The purpose of the present invention is to provide a free vortex wake prediction method for the aerodynamic performance of large wind turbine blades, effectively utilizing the vortex dynamics theory to numerically simulate the strength and spatial distribution of the vortex structure of large wind turbine blades, accurately simulating the downwash effect of the wake vortex system on the blade, and further reducing the number of iterative steps and improving the overall computational efficiency and prediction accuracy of the unsteady aerodynamic performance of large wind turbine blades by proposing a new attached vortex iteration method and a time marching method, so as to solve the safety problems caused by inaccurate prediction of blade dynamic loads.

[0006] To achieve the above technical objectives, the technical solution adopted by the present invention is as follows:

[0007] A free vortex wake prediction method for the aerodynamic performance of large wind turbine blades, comprising the following steps:

[0008] S1, the program reads the operating parameters, reference airfoil aerodynamic data, and blade geometric feature data required for the calculation loop from the data file;

[0009] S2. Use the momentum blade element method or the prescribed vortex wake method to estimate the initial circulation distribution of each blade section, and determine the spatial shape of the wake surface and the initial circulation distribution.

[0010] S3. Update the circumferential angle of the blade and the spatial coordinates of the attached vortex nodes on each section in the inertial system.

[0011] S4. Calculate the induced velocity on the wake surface at the new time.

[0012] S5. Combine the wind speed data and the induced velocity to update the coordinates of the wake vorticity element nodes, and realize the rotation and deformation of the wake.

[0013] S6. Use an efficient iterative method to update the circulation distribution of the blade attached vorticity elements, and obtain the blade aerodynamic characteristic results at the current time.

[0014] S7. If it is in the full unsteady dynamic aerodynamic performance prediction mode, as the inflow wind speed is updated, return to S3 to enter the loop calculation until the predetermined maximum number of time steps is reached.

[0015] In step S1, read the operating parameters required for the calculation loop, the aerodynamic data of the reference airfoil, and the blade geometric characteristic data from the data file by executing the program. The operating parameters read in step S1 include: the calculation parameters of the free vortex wake model and the inflow wind speed V at the hub height of the wind turbine in the wind field environment. ∞ ; The aerodynamic data of the reference airfoil includes: the lift coefficient C at different angles of attack of the reference airfoil. l and the drag coefficient C. d The blade geometric characteristics include: the chord length c on blade section j. j and the twist angle θ. j .

[0016] In step S2, using the momentum blade element method or the prescribed vortex wake method to estimate the initial circulation distribution of each blade section further includes:

[0017] The steps of using the momentum blade element method to estimate the initial circulation distribution of each blade section are:

[0018] S2-1 Divide the entire blade into N segments along the radial direction. On the j-th section, obtain the effective angle of attack α of the section according to the inflow angle of the air flow calculated by the velocity triangle and the twist angle θ of the section. j Based on the aerodynamic coefficient sequence of the airfoil of each reference section in the spanwise direction of the blade, interpolate to obtain the lift coefficient C per unit span on the j-th section. j , the drag coefficient C. l,j , d,j ,

[0019]

[0020] Cl,j = interp(α, C l , α j )

[0021] C d,j = interp(α, C d , α j )

[0022] In the formula: V axis is the axial velocity component of the wind turbine rotor; ω is the rotational angular velocity of the wind turbine rotor; r is the radial coordinate of the cross-section; a is the axial induction factor; a' is the tangential induction factor on the rotating plane of the wind turbine rotor. The determination methods of the axial and tangential induction factors are: determined by fixed-point iteration in combination with the Prandtl tip correction formula and the Glauert correction method; interp is a cubic spline interpolation program, α, C l , and C d are respectively the geometric angle of attack of the reference airfoil and the corresponding lift and drag coefficient sequences,

[0023] S2-2 converts the lift coefficient C l,j of the cross-section j calculated above into the two-dimensional vortex segment circulation Γ j at the 1 / 4 chord length of the cross-section,

[0024]

[0025] In the formula, V j is the absolute value of the relative wind speed at the blade cross-section j,

[0026] Use the predetermined vortex wake method to estimate the initial circulation distribution Γ j of each blade cross-section, and directly read the circulation distribution data calculated by the previous predetermined vortex wake method to assign values to the circulation Γ of the blade attached vortex and the vortex elements of the wake vortex.

[0027] In step S3, according to the input number of circumferential angle steps N azim update the circumferential angle ψ n+1 of the blade at the n+1 time step and the spatial coordinates x bound，

[0028]

[0029] ψ n+1 = ψ n + dψ

[0030] ψ n is the circumferential angle of the blade at the n time step.

[0031] In step S4, the induced velocity V θ on the wake vortex surface is calculated as follows:

[0032] S41. Divide the wake vortex structure into several straight vortex segments, and obtain the induced velocity component V at any vortex surface node in space according to the inviscid Biot-Savart law. inv ,

[0033]

[0034] In the formula: and are the length vectors from the two end nodes of the straight vortex segment to any vortex surface node x n ; r1 and r2 are their length values respectively.

[0035] S42. Regularize the inviscid induced velocity value V inv to obtain the induced velocity V θ on the wake vortex surface.

[0036]

[0037] In the formula: H is the perpendicular distance from the vortex surface node to the vortex segment, h0 is the initial vortex core radius, and the value of the exponent e is calibrated by the hot-wire anemometry test in the near wake.

[0038] In step S5, the method for updating the coordinates of the wake vortex surface vorticity element nodes by combining the wind speed data and the induced velocity is as follows: The coordinates x n+1 of the wake vortex surface vorticity element nodes at the (n + 1)-th time step are obtained by accumulating the oncoming flow wind speed V n of each node, the total induced velocity V ∞ caused by the wake vortex, and integrating on the basis of the coordinates x induced at the n-th time step; The time integration adopts the second-order predictor-corrector format or the third-order three-step Runge-Kutta format.

[0039] Among them, the second-order predictor-corrector format is calculated according to the following formula.

[0040] x pred = x n + Δt(V ∞ + V induced (x n ))

[0041]

[0042] In the formula: x pred is the coordinate prediction value; is the coordinate of the wake vortex surface vorticity element node at the (n + 1)-th time step obtained by the second-order predictor-corrector format; Δt is the time step size.

[0043] The third-order three-step Runge-Kutta format is calculated according to the following formula.

[0044] V(x) = V ∞ + V induced (x)

[0045] x (1) = x n + ΔtV(x n )

[0046]

[0047] where: x (1) 、x (2) are the intermediate values at each step of Runge - Kutta; is the coordinate of the wake surface vortex element node at the (n + 1)-th time step obtained by the third - order three - step Runge - Kutta format,

[0048] The coordinate x of the wake surface vortex element node at the (n + 1)-th time step calculated by the above second - order predictor - corrector format or third - order three - step Runge - Kutta format n+1 is used to represent the spatial position of the vortex surface and the update calculation of the induced velocity in the iteration of the next time step S4.

[0049] In step S6, the specific method for updating the circulation distribution of the blade attached vortex elements using an efficient iteration method is as follows:

[0050] The circulation distribution of the blade attached vortex elements at the (n + 1)-th time step The efficient iterative calculation includes the following sub - steps:

[0051] S61: Given the wake vortex element circulation , solve the first predicted value of the attached vortex element circulation and the updated value dΓ * ,

[0052]

[0053] where: the input parameters and of the lifting line solver LiftingLineSolver represent the circulation value of the blade attached vortex elements and the circulation value of the wake vortex elements at the n - th time step, respectively;

[0054] S62, use the above - mentioned circulation predicted value to solve the second time with the lifting line model, and obtain the second predicted value of the attached vortex element circulation and the updated value dΓ^,

[0055]

[0056] S63, use the attached vortex element circulation and the above two circulation update values to obtain the bound vortex element circulation at the new time step n+1 for the corrected distribution, as shown in the following formula

[0057]

[0058] where; w is the relaxation factor, and its value range is between 0.6 and 1.0. It can be adjusted according to the aerodynamic performance in the blade span direction through numerical experiments to promote iterative convergence

[0059] The bound vortex element circulation at the new time step n+1 obtained in step S6 and the wake vortex element node coordinates x obtained in step S5 at the new time step n+1 are used for calculating the induced velocity at the vortex element nodes on the wake surface in step S4 at the next moment, and based on the bound vortex element circulation and the wake vortex element circulation the local downwash induced velocity V at the 1 / 4 chord length of each blade section is calculated, and the magnitudes of the axial induction factor a and the tangential induction factor a' are obtained, so as to obtain the rising drag coefficient distribution of the section in the way of step S2, and the blade torque M, axial thrust Th and power P values at this moment are calculated by numerical integration using the following formula induced where: C

[0060]

[0061] and C n,j and C t,j are the normal force and axial force coefficients at section j in turn; q j is the local dynamic pressure at section j; r j is the radial coordinate of the midpoint of the jth blade section interval; Δr j is the radial width of the jth blade section interval, and ρ is the air density

[0062] Using the correction formula shown in step 6 for the bound circulation value at the new time step for the iterative method can significantly improve the convergence of the bound vortex solution iterative cycle and significantly reduce the total number of iterative steps

[0063] In step S7, when performing full unsteady dynamic aerodynamic performance prediction, as the inflow wind speed V ∞ is updated at each time step, return to S3 to enter the loop calculation until the predetermined maximum number of time steps N t,max is reached, and the operation is stopped

[0064] Compared with the prior art, the beneficial effects of the present invention are as follows:

[0065] First, the method for predicting the aerodynamic performance of a wind turbine blade based on a free vortex wake according to the present invention can comprehensively predict the circulation and three-dimensional spatial distribution of the wake vortices of the wind turbine blade, taking into account the viscous effects near the vortex core.

[0066] Second, the method for predicting the aerodynamic performance of a wind turbine blade based on a free vortex wake according to the present invention has obvious advantages in terms of the iteration efficiency of calculating the attached vortex circulation, and is improved in terms of the stability of time advancement in the unsteady calculation condition, which is beneficial to the popularization and application of the present invention.

[0067] Third, the time data series of the wind turbine torque, power, and axial thrust calculated by the method of the present invention can be used for the dynamic aerodynamic performance evaluation of the blades of a large wind turbine under given inflow wind conditions, which is of great significance for predicting the extreme loads and safety performance of the wind turbine blades. Description of the Drawings

[0068] Figure 1 is a flowchart of the free vortex wake prediction method for the aerodynamic performance of the present invention;

[0069] Figure 2 is a basic layout diagram of the wake vortex surface generated by the vortex wake method;

[0070] Figure 3 is a comparison diagram of the convergence step accuracy of the attached vortex iteration step S6 of the present invention and the conventional fixed-point method;

[0071] Figure 4 is a comparison diagram of the power coefficient of the example of the present invention. Detailed Embodiment

[0072] The following further describes the embodiments of the present invention in detail with reference to the drawings.

[0073] See Figure 1 , the present invention discloses a free vortex wake prediction method for the aerodynamic performance of the blades of a large wind turbine, including the following steps:

[0074] S1, Read in the operating parameters required for the calculation loop, the aerodynamic data of the reference airfoil, and the blade geometric characteristics;

[0075] S2, Use the momentum blade element method or the predetermined vortex wake method to estimate the initial circulation distribution of each section of the blade, and determine the spatial shape of the wake vortex surface and the initial circulation distribution;

[0076] S3, Update the circumferential angle of the blade and the spatial coordinates of the attached vortex nodes on each section in the inertial system;

[0077] S4, Calculate the induced velocity field on the wake vortex surface at the new moment;

[0078] S5. Combine the wind speed data and the induced velocity to update the coordinates of the vorticity element nodes on the wake vortex surface, realizing the rotation and deformation of the wake vortex;

[0079] S6. Use an efficient iterative method to update the circulation distribution of the blade attached vortex system, and obtain the blade aerodynamic characteristic results at the current moment;

[0080] S7. If it is in the full unsteady dynamic aerodynamic performance prediction mode, as the incoming wind speed is updated, return to S3 to enter the loop calculation until the predetermined maximum number of time steps.

[0081] In step S1, the input operation parameters include: the wind turbine diameter, the blade profile geometric parameters (radial coordinates, chord lengths, twist angles, etc. of each reference section); the lift coefficient and drag coefficient sequences at different angles of attack of the reference airfoil; the main calculation parameters of the free vortex wake model, and the incoming wind speed V at the hub height of the wind turbine required for the calculation in the wind field environment ∞ 。

[0082] In step S2, use the momentum blade element method or the predetermined vortex wake method to estimate the initial circulation distribution of each blade section, which further includes:

[0083] The steps to estimate the initial circulation distribution of each blade section using the momentum blade element method are:

[0084] S2-1 Divide the entire blade into N segments along the radial direction, and respectively obtain the effective angle of attack α of the section according to the inflow angle of the air flow calculated by the velocity triangle and the twist angle θ of the section and the twist angle θ of the section j at each section, and obtain the effective angle of attack α of the section j , and interpolate based on the aerodynamic coefficient sequences of the airfoils of each reference section in the spanwise direction of the blade to obtain the lift coefficient C per unit span at the j-th section l,j 、drag coefficient C d,j ,

[0085]

[0086] C l,j =interp(α,C l ,α j )

[0087] C d,j =interp(α,C d ,α j )

[0088] where: V axis is the axial velocity component of the wind turbine; ω is the rotational angular velocity of the wind turbine; r is the radial coordinate of the section; a is the axial induction factor; a' is the tangential induction factor on the wind turbine rotation plane. The axial and tangential induction factors need to be determined by fixed-point iteration in combination with the Prandtl tip correction formula and the Glauert correction method.

[0089] S2-2 Convert the lift coefficient C of section j calculated above l,j into the circulation of the two-dimensional vortex segment at the 1 / 4 chord length of the section by using the theorem

[0090]

[0091] Use the predetermined vortex wake method to estimate the initial circulation distribution Γ of each section of the blade j , then directly read in the circulation distribution calculated by the previous predetermined vortex wake method to assign values to the blade attached vortex and wake vortex elements.

[0092] In step S3, according to the input number of circumferential angle steps N azim Update the circumferential angle ψ of the blade at the n+1 time step n+1 and the spatial coordinates x of the attached vortex nodes on each section in the inertial system bound .

[0093]

[0094] In step S4, it further includes:

[0095] S41, Divide the wake surface into several straight vortex segments, and obtain the induced velocity components V of any vortex surface node in space according to the inviscid Biot-Savart law inv ,

[0096]

[0097] In the formula: and are the length vectors from the endpoints of the straight vortex segment to the spatial point respectively, and r1 and r2 are their length values.

[0098] S42, Regularize the inviscid induced velocity V inv by using the viscous model to obtain the induced velocity V of the node θ ,

[0099]

[0100] In the formula: H is the perpendicular distance from the node to the vortex segment, h0 is the initial vortex core radius, and the value of the exponent e, both of which can be calibrated by the near-wake hot-wire anemometry test.

[0101] Step S5 further includes:

[0102] The coordinates x of the vortex element nodes on the wake surface at the n+1 time step n+1 are based on the coordinates x at the n time step n by accumulating the oncoming flow velocity V of each node ∞ , the vortex-induced velocity V inducedAnd integrating to obtain; in addition to using the second-order predictor-corrector format (PC2) for time integration, the third-order three-step Runge-Kutta format (RK33) can be used for calculation starting from the stability requirements.

[0103] Among them, the second-order predictor-corrector format is

[0104] x pred = x n + Δt(V ∞ + V induced (x n ))

[0105]

[0106] In the formula: x pred is the coordinate predicted value; is the new time-step coordinate obtained by the PC2 format; Δt is the time step size.

[0107] The third-order three-step Runge-Kutta format (RK33) is calculated according to the following formula

[0108] V(x) = V ∞ + V induced (x))

[0109] x (1) = x n + ΔtV(x n )

[0110]

[0111] In the formula: x (1) , x (2) are the intermediate values of each step of Runge-Kutta; is the new time-step coordinate obtained by the RK33 format; Δt is the time step size. For large wind turbine blades with different aerodynamic layouts, the initial circulation of the free-vortex wake method may have a large difference from the convergence result after multiple cycles; the RK33 format can be selected to improve the stability of time advancement under the condition of constant time step size.

[0112] Step S6 further includes the following sub-steps:

[0113] S61: On the premise of a given wake vortex element strength , solve the first predicted value of the attached vortex element circulation and the updated value dΓ * ,

[0114]

[0115] In the formula: the input parameters of the LiftingLineSolver and respectively and successively represent the blade bound circulation value and the circulation distribution of the wake vortex elements at the nth time step;

[0116] S62, using the above circulation prediction value Solve for the second time with the lifting line model to obtain the secondary prediction value of the bound vortex element circulation and the updated value dΓ^,

[0117]

[0118] S63, using the bound vortex element circulation at the nth step and the above two circulation updated values to obtain the corrected distribution of the bound vortex element circulation at the new time step n+1, as shown in the following formula In the formula; w is the relaxation factor, and its value range is between 0.6 and 1.0. It can be adjusted according to the blade spanwise aerodynamic performance through numerical experiments to promote iterative convergence.

[0119]

[0120] In the formula; w is the relaxation factor, and its value range is between 0.6 and 1.0. It can be adjusted according to the blade spanwise aerodynamic performance through numerical experiments to promote iterative convergence.

[0121] Adopting the correction formula shown in step 6 for the bound circulation value at the new time step The iterative method can significantly improve the convergence of the bound vortex solution loop iteration and significantly reduce the total number of iterative steps.

[0122] In step S7, when performing full unsteady dynamic aerodynamic performance prediction, as the inflow wind speed V ∞ is updated at each time step, return to step S3 to enter the loop calculation until the predetermined maximum number of time steps N t,max is reached, that is, stop the operation. The time data series such as wind turbine torque, power, and axial thrust calculated by the above method can be used for the dynamic aerodynamic performance evaluation of large wind turbine blades under given inflow wind conditions.

[0123] The present invention will be further described below with reference to the accompanying drawings. The following embodiments are only used to more clearly illustrate the technical solutions of the present invention and cannot be used to limit the protection scope of the present invention.

[0124] Embodiment 1: The number of blades of the wind turbine is 3, the wind turbine diameter is 63.0 m, the blade root radius is 1.5 m, the blade chord length and twist angle distribution are read from the original data file, the calculated wind speed range is 5.5 - 16.5 m / s, and the wind turbine tip speed ratio range is 4.6 - 13.7.

[0125] The blade aerodynamic model is replaced by straight vortex segments arranged at the 1 / 4 chord length position of each section, and forms a quadrilateral vortex ring form with the trailing vortex segments and shedding vortex segments downstream (see Figure 2 ).

[0126] Referring to Figure 3 , the numerical simulation results obtained by using the free vortex wake wind turbine blade aerodynamic performance evaluation method proposed by the present invention are compared with the reference method data. On the premise of not reducing the calculation accuracy, the number of iterations required for attached vortex calculation is significantly reduced, the overall calculation efficiency is improved, and the stability of the program time advancement is significantly improved (the calculation results are shown in Figure 4 ).

[0127] In summary, the method proposed by the present invention can effectively simulate the wake vortex circulation and spatial distribution of the wind turbine blade wake flow field, can accurately simulate the downwash effect of the wake vortex structure on the blade, and is sufficient to meet the engineering requirements in terms of the prediction accuracy of the wind turbine blade aerodynamic performance.

[0128] Those skilled in the art should understand that the embodiments of the present application can be provided as methods, systems, or computer program products. Therefore, the present application can adopt the form of a complete hardware embodiment, a complete software embodiment, or an embodiment combining software and hardware aspects. Moreover, the present application can adopt the form of a computer program product implemented on one or more computer-usable storage media (including but not limited to disk storage, CD-ROM, optical storage, etc.) containing computer-usable program code. The solutions in the embodiments of the present application can be implemented in various computer languages, for example, object-oriented compiled programming language C++ and scripting language Python, etc.

[0129] The present application is described with reference to the flowcharts and / or block diagrams of methods, devices (systems), and computer program products according to the embodiments of the present application. It should be understood that each process and / or block in the flowchart and / or block diagram, and the combination of processes and / or blocks in the flowchart and / or block diagram, can be implemented by computer program instructions. These computer program instructions can be provided to the processor of a general-purpose computer, a special-purpose computer, an embedded processor, or other programmable data processing devices to generate a machine, so that the instructions run by the processor of the computer or other programmable data processing devices generate means for realizing the functions specified in Figure 1 one process or multiple processes and / or blocks Figure 1 one block or multiple blocks.

[0130] These computer program instructions can also be stored in a computer-readable memory that can direct a computer or other programmable data processing device to work in a specific manner, so that the instructions stored in the computer-readable memory generate a manufactured article including an instruction device, and the instruction device realizes the functions in the process Figure 1One or more processes and / or blocks Figure 1 The functions specified in one or more blocks.

[0131] These computer program instructions can also be loaded onto a computer or other programmable data processing device, so that a series of operation steps are performed on the computer or other programmable device to generate a computer-implemented process. Thus, the instructions executed on the computer or other programmable device provide steps for implementing the functions specified in one Figure 1 One or more processes and / or blocks Figure 1 The steps of the functions specified in one or more blocks.

[0132] Although the preferred embodiments of the present application have been described, those skilled in the art can make additional changes and modifications to these embodiments once they learn the basic creative concepts. Therefore, the appended claims are intended to be construed to include the preferred embodiments as well as all changes and modifications falling within the scope of the present application.

[0133] Obviously, those skilled in the art can make various changes and modifications to the present application without departing from the spirit and scope of the present application. Thus, if these modifications and variations of the present application fall within the scope of the claims of the present application and their equivalent technologies, the present application is also intended to include these modifications and variations.

Claims

1. A free-vortex wake prediction method for the aerodynamic performance of large wind turbine blades, characterized in that, It includes the following steps: S1. The execution program reads the operating parameters required for the calculation loop, the aerodynamic data of the reference airfoil, and the blade geometric feature data from the data file; S2. Use the momentum blade element method or the prescribed vortex wake method to estimate the initial circulation distribution of each blade section, determine the spatial shape of the wake surface and the initial circulation distribution; S3. Update the circumferential angle of the blade and the spatial coordinates of the attached vortex nodes on each section in the inertial system; S4. Calculate the induced velocity on the wake surface at the new moment; S5. Combine the wind speed data and the induced velocity to update the coordinates of the wake surface vortex element nodes to realize the rotation and deformation of the wake; S6. Use an efficient iterative method to update the circulation distribution of the blade attached vortex elements to obtain the blade aerodynamic characteristic results at the current moment; S7. If it is a full unsteady dynamic aerodynamic performance prediction mode, as the inflow wind speed is updated, return to S3 to enter the loop calculation until the predetermined maximum number of time steps is reached.

2. A free-vortex wake prediction method for the aerodynamic performance of large wind turbine blades according to claim 1, characterized in that In step S1, the operating parameters required for the calculation loop, the aerodynamic data of the reference airfoil, and the blade geometric feature data are read from the data file through an execution program. The operating parameters read in step S1 include: the calculation parameters of the free-vortex wake model and the inflow wind speed V at the hub height of the wind turbine in the wind field environment ∞ ; the aerodynamic data of the reference airfoil includes: the lift coefficient C at different angles of attack of the reference airfoil l and the drag coefficient C d , the blade geometric features include: the chord length c at the blade section j j and the twist angle θ j .

3. A free-vortex wake prediction method for the aerodynamic performance of large wind turbine blades according to claim 2, characterized in that In step S2, using the momentum blade element method or the prescribed vortex wake method to estimate the initial circulation distribution of each blade section further includes: The steps of using the momentum blade element method to estimate the initial circulation distribution of each blade section are: S2-1 divides the entire blade into N segments radially, and calculates the inflow angle of the air flow obtained according to the velocity triangle on the j-th section and the twist angle θ of the section j to obtain the effective angle of attack α of the section j , and interpolates based on the aerodynamic coefficient sequence of the airfoil of each reference section in the spanwise direction of the blade to obtain the lift coefficient C per unit span on the j-th section l,j , drag coefficient C d,j , C l,j = interp(α, C l , α j ) C d,j = interp(α, C d , α j ) Where: V axis is the axial velocity component of the wind turbine rotor; ω is the rotational angular velocity of the wind turbine rotor; r is the radial coordinate of the cross-section; a is the axial induction factor; a' is the tangential induction factor on the rotating plane of the wind turbine rotor. The determination methods of the axial and tangential induction factors are as follows: they are determined by fixed-point iteration in combination with the Prandtl tip correction formula and the Glauert correction method; interp is a cubic spline interpolation program, α, C l , and C d are respectively the geometric angle of attack of the reference airfoil and the corresponding lift and drag coefficient sequences, S2-2 Convert the lift coefficient C of cross-section j calculated above l,j into the circulation Γ of the two-dimensional vortex segment at the 1 / 4 chord length of the cross-section j , where V j is the absolute value of the relative wind speed at blade section j, Estimate the initial circulation distribution Γ of each blade section using the predetermined wake method j , directly read in the circulation distribution data calculated by the previous predetermined wake method to assign the circulation Γ of the vortex segments of the blade bound vortex and wake vortices.

4. A free-vortex wake prediction method for the aerodynamic performance of large wind turbine blades according to claim 3, characterized in that, In step S3, according to the input number of circumferential angle steps N azim update the blade circumferential angle ψ at the n+1 time step n+1 and the spatial coordinates x of the bound vortex nodes on each cross-section in the inertial system bound , ψ n is the circumferential angle of the blade at n time steps.

5. A free-vortex wake prediction method for the aerodynamic performance of blades of a large wind turbine according to claim 4, characterized in that In step S4, the induced velocity V on the trailing vortex surface θ is calculated according to the following steps: S41. Divide the wake vortex structure into several straight vortex segments, and obtain the induced velocity component V of any vortex surface node in space according to the inviscid Biot-Savart law inv , Wherein: and are respectively the length vectors from the nodes at both ends of the straight vortex segment to any vortex surface node x n ; r1 and r2 are respectively their length values; S42. Regularize the inviscid induced velocity value V using a viscous model to obtain the induced velocity V on the wake surface inv θ ,​ Where: H is the perpendicular distance from the vortex surface node to the vortex segment, h0 is the initial vortex core radius, and the value of the exponent e is calibrated by the near wake hot wire anemometry test.

6. A free-vortex wake prediction method for the aerodynamic performance of large wind turbine blades according to claim 5, characterized in that In step S5, the method for updating the coordinates of the wake surface vortex element nodes by combining the wind speed data and the induced velocity is as follows: the coordinate x of the wake surface vortex element nodes at the n+1 time step n+1 at the coordinate x at the n time step n is obtained by accumulating the oncoming wind speed V of each node ∞ and the total induced velocity V caused by the wake induced and integrating; the time integration adopts the second-order predictor-corrector format or the third-order three-step Runge-Kutta format. Among them, the second-order prediction-correction format is calculated according to the following formula, x pred = x n + Δt(V ∞ + V induced (x n )) where: x pred is the coordinate prediction value; is the coordinate of the wake surface vortex element node at the (n + 1)-th time step obtained by the second-order prediction-correction format; Δt is the time step size; The third-order three-step Runge-Kutta format is calculated according to the following formula, V(x) = V ∞ + V induced (x) x (1) = x n + ΔtV(x n ) where: x (1) , x (2) are the intermediate values at each step of Runge-Kutta; is the coordinate of the trailing vortex panel vorticity element node at the (n + 1)-th time step obtained by the third-order three-step Runge-Kutta scheme, The coordinates x of the trailing vortex surface vorticity element nodes at the (n + 1)-th time step calculated using the above second-order predictor-corrector format or third-order three-step Runge-Kutta format n+1 are used to represent the spatial position of the vortex surface and the updated calculation of the induced velocity in the iteration of the next time step in step S4.

7. A free-vortex wake prediction method for the aerodynamic performance of large wind turbine blades according to claim 6, characterized in that In step S6, the specific method of using an efficient iterative method to update the circulation distribution of the blade attached vortex elements is: The circulation distribution of the blade attached vorticity element at the (n + 1)-th time step The efficient iterative calculation thereof includes the following sub-steps: S61: Given the trailing vortex element annulus Under the premise of solving the first estimate of the attached vortex annulus according to the lift line model, and update value dΓ * , Where: the input parameters of the LiftingLineSolver and respectively and sequentially represent the circulation value of the blade attached vortex element and the circulation value of the wake vortex element at the n-th time step; S62, using the above-mentioned circulation pre-estimation value Solve for the second time using the lifting-line model to obtain the secondary pre-estimation value of the circulation of the attached vorticity element and the updated value dΓ ^ , S63. Utilize the circulation of the attached vortex element at the n-th step and the above two circulation update values to obtain the corrected distribution of the circulation of the attached vortex element at the (n + 1)-th step at the new moment, as shown in the following formula as follows Where; w is the relaxation factor, and its value range is between 0.6 - 1.

0. It can be adjusted according to the blade spanwise aerodynamic performance through numerical experiments to promote iterative convergence, The bound vortex element circulation at the new time step n+1 obtained in step S6 and the coordinates x of the trailing vortex element nodes at the new time step obtained in step S5 n+1 are used to calculate the induced velocity at the vortex element nodes on the trailing vortex surface in step S4 at the next time step, and based on the bound vortex element circulation and the trailing vortex element circulation the local downwash induced velocity V at the 1 / 4 chord length of each blade section is calculated from the distribution, induced as well as the magnitudes of the axial induction factor a and the tangential induction factor a', so as to obtain the distribution of the sectional lift drag coefficient in the manner of step S2, and the blade torque M, axial thrust Th and power P values at this time are calculated by numerical integration using the following formula Where: C n,j and C t,j are the normal force and axial force coefficients at section j, respectively; q j is the local dynamic pressure at section j; r j is the radial coordinate of the midpoint of the j-th blade section interval; Δr j is the radial width of the j-th blade section interval, and ρ is the air density.

8. A free-vortex wake prediction method for the aerodynamic performance of large wind turbine blades according to claim 7, characterized in that In step S7, when performing the full unsteady dynamic aerodynamic performance prediction, as the incoming flow velocity V ∞ is updated at each time step, return to S3 to enter the loop calculation until the predetermined maximum number of time steps N t,max is reached, i.e., stop the operation.

Images