A large wind turbine blade aerodynamic performance free vortex wake prediction method

Abstract

The application discloses a large wind turbine blade aerodynamic performance free vortex wake prediction method, which uses a momentum blade element method or a predetermined vortex wake method to predict initial circulation distribution of each section of the blade, determine the spatial shape of the wake surface and the initial circulation distribution; update the circumferential angle of the blade and the spatial coordinates of the attached vortex nodes on each section in the inertial system; calculate the induced velocity on the wake surface at the new time; rotate and deform the wake; update the attached vortex element circulation distribution of the blade to obtain the current time blade aerodynamic characteristic result; with the inflow wind speed update, return to S3 to enter the cycle calculation until the predetermined maximum time step is reached. The application can significantly reduce the iteration number of the attached vortex circulation distribution at each time step, which is beneficial to improve the calculation efficiency of the wind turbine dynamic aerodynamic characteristic prediction based on the free vortex wake numerical method, thereby helping to quickly predict the dynamic aerodynamic performance of the large wind turbine blade and diagnose the safety under batch wind conditions.

Description

TECHNICAL FIELD

[0001] The present application relates to the field of new energy wind power generation technology, in particular to a large wind turbine blade aerodynamic performance free vortex wake prediction method. BACKGROUND

[0002] The free vortex wake numerical algorithm combines the inviscid potential lifting line / surface theory with the vortex dynamics method, obtains the blade dynamic aerodynamic load information by predicting the evolution law of the wind turbine blade wake surface and the induced effect on the wind turbine section, and can balance the prediction accuracy and calculation economy.

[0003] Compared with the existing numerical calculation means, the free vortex wake numerical algorithm has obvious advantages in the prediction accuracy of vortex induction compared with the momentum blade element method; compared with the Reynolds average Navier-Srokes equation solution or the large eddy simulation CFD method, the free vortex wake algorithm can effectively control the calculation time required for engineering evaluation, and provides a high-efficiency and accurate evaluation method for the dynamic aerodynamic load of large wind turbine blades under complex inflow conditions, provides aerodynamic data support for the aeroelastic response analysis of large wind turbine blades, and helps the extreme aerodynamic load prediction of large and super large blades and the evaluation of structural aeroelastic stability and operation safety.

[0004] At present, the existing free vortex wake numerical algorithm has the problems of large number of iteration steps, poor overall calculation efficiency and prediction accuracy of the unsteady aerodynamic performance of large wind turbine blades, and safety problems caused by inaccurate blade dynamic load prediction. SUMMARY

[0005] The purpose of the present application is to provide a large wind turbine blade aerodynamic performance free vortex wake prediction method, effectively utilize the vortex dynamics theory to numerically simulate the strength and spatial distribution of the vortex structure of large wind turbine blades, accurately simulate the downwash effect of the wake vortex system on the blade, further reduce the number of iteration steps by proposing a new attached vortex iteration method and time marching method, improve the overall calculation efficiency and prediction accuracy of the unsteady aerodynamic performance of large wind turbine blades, and solve the safety problems caused by inaccurate blade dynamic load prediction.

[0006] To achieve the above technical purpose, the technical scheme adopted by the present application is as follows:

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

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

[0009] S2, estimate the initial circulation distribution of each section of the blade by using the momentum blade element method or the prescribed vortex wake method, 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, update the coordinates of the vortex nodes on the wake surface by combining the wind speed data and the induced velocity, realize the rotation and deformation of the wake;

[0013] S6, update the circulation distribution of the attached vortex elements of the blade by using the high-efficiency iteration method, obtain the aerodynamic characteristic results of the blade at the current time;

[0014] S7, if it is a full unsteady dynamic aerodynamic performance prediction mode, return to S3 for loop calculation until the maximum time step is reached.

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

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

[0017] The step of estimating the initial circulation distribution of each section of the blade by using the momentum blade element method is:

[0018] S2-1 divide the entire blade into N sections along the radial direction, obtain the effective angle of attack α j on the j-th section according to the calculated flow inflow angle and the twist angle θ j of the section, respectively, interpolate the lift coefficient C l,j and the drag coefficient C d,j of the unit span on the j-th section based on the aerodynamic force coefficient sequence of the blade spanwise reference section airfoil,

[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; ω is the angular velocity of the wind turbine; r is the radial coordinate of the cross section; a is the axial induction factor; a' is the tangential induction factor on the wind turbine's rotational surface. The axial and tangential induction factors are determined by combining the Prandtl tip correction formula and the Glauber's correction method through fixed-point iteration; interp is the cubic spline interpolation program, α, C l , and C d These are the geometric angles of attack of the reference airfoil and the corresponding lift and drag coefficient sequences.

[0023] S2-2 The lift coefficient C of section j calculated above. l,j Converted to two-dimensional vortex circulation Γ at 1 / 4 chord length of cross section j ,

[0024]

[0025] In the formula V j Let be the absolute value of the relative wind speed at blade section j.

[0026] Predicting the initial circulation distribution Γ of each blade cross section using the predetermined vortex wake method j The circulation distribution data calculated by the previously predetermined vortex wake method is directly read in to assign values ​​to the vortex segment circulation Γ of the blade attached vortex and wake vortex element.

[0027] In step S3, the number of circumferential angle steps N is calculated based on the input. azim Update the blade circumferential angle ψ at time step n+1. n+1 And the spatial coordinates x of the vortex nodes attached to each cross section in the inertial frame. bound，

[0028]

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

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

[0031] In step S4, the induced velocity V on the wake surface θ Calculate using the following steps:

[0032] S41, the wake vortex structure is divided into several straight vortex segments, and the induced velocity component V of any vortex surface node in space is obtained according to the inviscid Biot-Savart law inv ,

[0033]

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

[0035] S42, the inviscid induced velocity value V inv is regularized by using a viscous model to obtain the induced velocity V θ on the wake vortex surface,

[0036]

[0037] In the formula: H is the vertical distance of the vortex surface node to the vortex segment, h0 is the initial vortex core radius, and the value of the index e is calibrated by a near wake hot wire speed test.

[0038] In step S5, the method for updating the coordinates of the wake vortex surface vortex element nodes in combination with the wind speed data and the induced velocity is: the coordinates x n+1 of the wake vortex surface vortex element nodes at the n+1 time step are obtained by accumulating the incoming flow wind speed V ∞ , the total induced velocity V induced caused by the wake vortex on the basis of the coordinates x n at the n time step and integrating; the time integration adopts a second-order prediction correction format or a third-order three-step Runge-Kutta format,

[0039] wherein the second-order prediction correction 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 coordinates of the wake vortex surface vortex element nodes at the n+1 time step obtained by the second-order prediction correction format; and Δt is the time step length;

[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] In the formula: x (1) x (2) The intermediate values ​​of each step in the Runge-Kutta series; These are the coordinates of the vortex element nodes at the (n+1)th time step of the wake surface, obtained using the third-order, three-step Runge-Kutta scheme.

[0048] The coordinates x of the vortex element nodes at the (n+1)th time step, calculated using either the second-order predictor-corrector scheme or the third-order three-step Runge-Kutta scheme. n+1 Used to represent the spatial position of the vortex surface in the next time step iteration and the update calculation of the induced velocity in step S4.

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

[0050] Circulation distribution of blade-attached vortex elements at time step n+1 Efficient iterative computation includes the following sub-steps:

[0051] S61: Given the wake vortex circulation Under the premise of solving the first estimated value of the circulation of the attached vortex element according to the lift line model. and update value dΓ * ,

[0052]

[0053] Where: Input parameters of the LiftingLineSolver solver and These represent the circulation values ​​of the blade-attached vortex element and the wake vortex element at the nth time step, respectively.

[0054] S62, using the above-mentioned circulation estimate The second solution using the lift line model yields a secondary estimate of the circulation of the attached vortex element. and update value dΓ^,

[0055]

[0056] S63, utilizing the nth step attached vortex circulation The new time step n+1 of the attached vortex element circulation is obtained by updating the circulation value twice above The modified distribution of the circulation is shown in the following formula

[0057]

[0058] In the formula, w is a relaxation factor, which can be selected from the range of 0.6-1.0, and can be adjusted according to the blade spanwise aerodynamic performance through numerical experiments to promote the convergence of iteration,

[0059] The new time step n+1 of the attached vortex element circulation obtained in step S6 The new time step n+1 of the attached vortex element circulation obtained in step S6 n+1 The calculation of the induced velocity at the vortex element node on the trailing vortex surface in the next time step S4, and the induced velocity is calculated based on the attached vortex element circulation The new time step n+1 of the attached vortex element circulation obtained in step S6 The induced velocity V induced The size of the axial induction factor a and the tangential induction factor a', so that the lift-drag coefficient distribution on the section is obtained in step S2, and the torque M, the axial thrust Th and the power P value of the blade at this time are calculated by numerical integration using the following formula,

[0060]

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

[0062] The iteration method shown in step 6 is used to modify the attached circulation value of the new time step The convergence of the attached vortex iteration cycle can be obviously improved, and the total iteration steps can be significantly reduced.

[0063] In step S7, when the full unsteady dynamic aerodynamic performance is predicted, the inflow wind speed V ∞ Is updated at each time step, and the calculation is returned to S3 to enter the loop until the predetermined maximum time step N t,max Is reached, i.e. the operation is stopped.

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

[0065] First, the free vortex wake based wind turbine blade aerodynamic performance method can comprehensively predict the circulation and three-dimensional spatial distribution of wind turbine blade wake vortex, and considers the viscous effect near the vortex core.

[0066] Second, the free vortex wake based wind turbine blade aerodynamic performance method has obvious advantages in iteration efficiency for attached vortex circulation calculation, and improves the stability of time marching in unsteady calculation conditions, which is beneficial to the popularization and application of the method.

[0067] Third, the time data sequence of wind turbine torque, power and axial thrust calculated by the method can be used for dynamic aerodynamic performance evaluation of large wind turbine blades under given inflow wind conditions, which has important significance for the prediction of extreme load and safety performance of wind turbine blades. BRIEF DESCRIPTION OF DRAWINGS

[0068] Figure 1 The flow chart of the free vortex wake prediction method of the aerodynamic performance of the application;

[0069] Figure 2 The basic layout diagram of the vortex wake method to generate the wake surface;

[0070] Figure 3 The convergence step number accuracy of the attached vortex iteration step S6 of the application compared with the conventional fixed point method;

[0071] Figure 4 The power coefficient comparison chart of the application example. DETAILED DESCRIPTION

[0072] The embodiments of the application will be further described in detail below with reference to the accompanying drawings.

[0073] Referring to Figure 1 The application discloses a large wind turbine blade aerodynamic performance free vortex wake prediction method, which comprises the following steps:

[0074] S1, read in the running parameters required for calculation cycle, reference airfoil aerodynamic data and blade geometric characteristics;

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

[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 surface at the new time;

[0078] S5, combining wind speed data and induced velocity, updating the coordinates of the vortex element nodes of the tail vortex surface, realizing the rotation and deformation of the tail vortex;

[0079] S6, updating the circulation distribution of the blade attached vortex system using an efficient iteration method to obtain the current time blade aerodynamic characteristic results;

[0080] S7, if it is a full unsteady dynamic aerodynamic performance prediction mode, return to S3 to enter the loop calculation with the inflow wind speed update until the predetermined maximum time step.

[0081] In step S1, the running parameters are read in, including: wind turbine diameter, blade shape geometric parameters (radial coordinates of each reference section, chord length, twist angle, etc.); lift coefficient and drag coefficient sequences of the reference airfoil at different angles of attack; main calculation parameters of the free vortex wake model and the inflow wind speed V ∞ .

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

[0083] The step of estimating the initial circulation distribution of each section of the blade using the momentum blade element method is:

[0084] S2-1 divides the entire blade into N sections along the radial direction, and calculates the airflow inflow angle and the twist angle θ j of each section according to the velocity triangle, j obtains the effective angle of attack α l,j of the section based on the aerodynamic force coefficient sequence of the blade spanwise reference airfoil, d,j ,

[0085]

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

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

[0088] V axis is the wind turbine axial velocity component; ω is the wind turbine angular velocity; 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 surface. The axial and tangential induction factors need to be determined by fixed point iteration combined with Prandtl tip correction formula and Glauert correction method.

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

[0090]

[0091] Predicts the initial circulation distribution Γ of each section of the blade by using the predetermined vortex wake method j , and directly reads in the circulation distribution of the blade attachment vortex and the vortex element of the trailing vortex calculated by the previous predetermined vortex wake method to assign values.

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

[0093]

[0094] In step S4, further includes:

[0095] S41, divides the trailing vortex surface into several linear vortex segments, and obtains the induced velocity component V inv of any vortex surface node in space according to the inviscid Biot-Savart law

[0096]

[0097] In the formula: and are the length vectors of the endpoints of the linear vortex segment to the space point, and r1 and r2 are the length values thereof.

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

[0099]

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

[0101] Step S5 further includes:

[0102] The coordinates x of the vortex element nodes of the trailing vortex surface at the n+1 time step n+1 are obtained by accumulating the incoming flow velocity V n , the vortex-induced velocity V ∞ , and the like of each node on the basis of the coordinates x at the n time step inducedAnd integral to get; time integral except using second order prediction correction format (PC2), can from stability requirements, using three order three step Runge-Kutta format (RK33) to calculate

[0103] Wherein second order prediction correction format is,

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

[0105]

[0106] In the formula: x pred For coordinate prediction value; For PC2 format to obtain new time step coordinate;Δt is time step length.

[0107] Three 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 + Δt V (x n )

[0110]

[0111] In the formula: x (1) , x (2) For Runge-Kutta intermediate value; For RK33 format to obtain new time step coordinate;Δt is time step length.For different aerodynamic arrangement Large wind turbine blade, free vortex wake method The initial circulation of the tail may be different from the convergence result after several cycles;In the case of constant time step, RK33 format can be selected to improve the stability of time advancement.

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

[0113] S61: under the premise of given tail vortex element intensity According to the lift line model, the first estimated value And update value dΓ * Of attached vortex element circulation are solved,

[0114]

[0115] Where: Input parameters of the LiftingLineSolver solver and The values ​​of blade attachment circulation and the circulation distribution of the wake element at time step n are respectively represented.

[0116] S62, using the above-mentioned circulation estimate The second solution using the lift line model yields a secondary estimate of the circulation of the attached vortex element. and update value dΓ^,

[0117]

[0118] S63, utilizing the nth step attached vortex circulation The above two circulation updates yield the attached vortex circulation at time n+1. The corrected distribution is shown in the following formula.

[0119]

[0120] In the formula, w is the relaxation factor, which can take values ​​between 0.6 and 1.0. It can be adjusted through numerical experiments according to the aerodynamic performance of the blade spanwise to promote iterative convergence.

[0121] The attachment loop value for the new time step is modified as shown in step 6. The iterative method can significantly improve the convergence of the iterative solution of attached vortices and significantly reduce the total number of iterations.

[0122] In step S7, when performing fully unsteady dynamic aerodynamic performance prediction, as the inflow velocity V... ∞ At each time step update, return to step S3 to enter the loop calculation until the predetermined maximum number of time steps N is reached. t,max The calculation is then stopped. The time data sequences of wind turbine torque, power, and axial thrust obtained by the above method can be used to evaluate the dynamic aerodynamic performance 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 solution of the present invention, and should not be used to limit the scope of protection of the present invention.

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

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

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

[0127] In summary, the method proposed in the present application can effectively simulate the wake vortex ring quantity 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 wind turbine blade aerodynamic performance prediction accuracy.

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

[0129] The present application is described with reference to flowcharts and / or block diagrams according to the methods, devices (systems), and computer program products of the embodiments of the present application. It should be understood that each flow and / or block in the flowcharts and / or block diagrams, and the combination of flows and / or blocks in the flowcharts and / or block diagrams 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 apparatus to produce a machine, so that the instructions executed by the processor of the computer or other programmable data processing apparatus produce a device that implements the functions specified in the flow Figure 1 The device that implements the functions specified in the flow Figure 1 one or more flows and / or blocks.

[0130] These computer program instructions can also be stored in a computer-readable memory that can cause the computer or other programmable data processing apparatus to work in a specific manner, so that the instructions stored in the computer-readable memory produce a manufactured product including instruction devices that implement the functions specified in the flow Figure 1one or more processes and / or blocks Figure 1 the function specified in the one or more blocks.

[0131] These computer program instructions can also be loaded into computer or other programmable data processing devices, so that a series of operation steps are performed on the computer or other programmable data processing devices to generate computer-implemented processes, thus the instructions running on the computer or other programmable data processing devices provide a process for implementing the flow Figure 1 one or more processes and / or blocks Figure 1 the function specified in the one or more blocks.

[0132] Although preferred embodiments of the application have been described, those skilled in the art will recognize that additional modifications and variations may be possible in light of the above disclosure. It is therefore intended that the appended claims cover all such modifications and variations as fall within the scope of the application.

[0133] Obviously, numerous modifications and variations of the present application are possible in light of the above teachings. It is therefore intended that the application be covered within the scope of the claims, and their equivalents.

Claims

1. A method for predicting the free vortex wake of aerodynamic performance of large wind turbine blades, characterized in that, Includes the following steps: S1, The program reads the running parameters, reference airfoil aerodynamic data, and blade geometric characteristic data required for the calculation cycle from the data file; 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 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 frame; S4, Calculate the induced velocity on the wake surface at the new moment; S5, combining wind speed data and induced velocity, updates the coordinates of the vortex element nodes on the wake vortex surface, realizing the rotation and deformation of the wake vortex; S6. Update the circulation distribution of the attached vortex elements on the blade using an efficient iterative method to obtain the blade aerodynamic characteristics at the current moment. S7, If it is a fully unsteady dynamic aerodynamic performance prediction mode, as the inflow velocity is updated, return to S3 to enter the loop calculation until the predetermined maximum number of time steps is reached; In step S5, the method for updating the coordinates of the vortex element nodes of the wake surface by combining wind speed data and induced velocity is as follows: the coordinates of the vortex element nodes of the wake surface at time step n+1. Coordinates at time step n Based on the sum of the incoming wind speeds at each node Total induced velocity caused by wake vortex The integral is obtained by integration; the time integral is obtained using a second-order predictor-corrector scheme or a third-order three-step Runge-Kutta scheme. The second-order prediction correction scheme is calculated using the following formula. ； ； In the formula: Estimated coordinate values; The coordinates of the vortex element nodes at the (n+1)th time step are obtained from the second-order prediction correction scheme. For time step; The inflow wind speed at the height of the wind turbine hub in a wind field environment; The third-order, three-step Runge-Kutta formula is calculated using the following formula. ； ； ； ； In the formula: , The intermediate values ​​of each step in the Runge-Kutta series; These are the coordinates of the vortex element nodes at the (n+1)th time step of the wake surface, obtained using the third-order, three-step Runge-Kutta scheme. The coordinates of the vortex element nodes at the (n+1)th time step, calculated using either the second-order predictor-corrector scheme or the third-order three-step Runge-Kutta scheme. Used to represent the spatial position of the vortex surface in the next time step iteration and the update calculation of the induced velocity in step S4.

2. The method for predicting the free vortex wake of aerodynamic performance of large wind turbine blades according to claim 1, characterized in that, In step S1, the program reads the operating parameters required for the calculation cycle, the aerodynamic data of the reference airfoil, and the geometric characteristic data of the blade from the data file. The operating parameters read in step S1 include: the calculation parameters of the free vortex wake model and the inflow wind speed at the height of the wind turbine hub under the wind field environment. The reference airfoil aerodynamic data includes: the lift coefficient of the reference airfoil at different angles of attack. and drag coefficient The geometric features of the blade include: chord length on the section of the blade and twist angle .

3. The method for predicting the free vortex wake of aerodynamic performance of large wind turbine blades according to claim 2, characterized in that, In step S2, the initial circulation distribution of each section of the blade is estimated using the momentum blade element method or the predetermined vortex wake method, further including: The steps for predicting the initial circulation distribution of different sections of a blade using the momentum leaf element method are as follows: S2-1 divides the entire blade radially into N segments, and in the... The airflow angles calculated from the velocity triangle on the cross section. and the angle of twist of the cross section Obtain the effective angle of attack of the cross section The first value is obtained by interpolation based on the aerodynamic coefficient sequence of each reference section airfoil along the blade span. Lift coefficient per unit span of cross section drag coefficient , ； ； ； ； In the formula: This refers to the axial velocity component of the wind turbine. The rotational angular velocity of the wind turbine; The radial coordinates of the cross section; It is the axial induction factor; The tangential induction factor on the rotor's rotating surface is determined by combining the Prandtl tip correction formula and the Glauber's correction method through fixed-point iteration. This is a cubic spline interpolation procedure. , and These are the geometric angles of attack of the reference airfoil and the corresponding lift and drag coefficient sequences. S2-2 The cross-section calculated above lift coefficient Converted to two-dimensional vortex circulation at 1 / 4 chord length of cross section , ； ； In the formula For the first The absolute value of the relative wind speed at the section of the blade. Predicting the initial circulation distribution of each section of the blade using the predetermined vortex wake method The circulation distribution data calculated using the previously predetermined vortex wake method is directly read into the vortex segment circulation of the blade attached vortex and wake vortex elements. Perform the assignment.

4. The method for predicting the free vortex wake of aerodynamic performance of large wind turbine blades according to claim 3, characterized in that, In step S3, the number of circumferential angle steps is determined based on the input. Update the blade circumferential angle at time step n+1. And the spatial coordinates of the vortex nodes attached to each cross section in the inertial frame. , ； ； Let be the circumferential angle of the blade at time step n.

5. The method for predicting the free vortex wake of aerodynamic performance of large wind turbine blades according to claim 4, characterized in that, In step S4, the induced velocity on the wake surface Calculate using the following steps: S41, the wake structure is divided into several straight vortex segments, and the induced velocity components at any node of the vortex surface in space are obtained according to the inviscid Biot-Savart law. , ； In the formula: and These represent the distances from the two ends of the straight vortex segment to any vortex surface node. The length vector; and These are their length values; S42, Using a viscous model to determine the inviscid-induced velocity values By regularizing, the induced velocity on the wake surface is obtained. , ； In the formula: This represents the vertical distance from the vortex surface node to the vortex segment. The initial vortex core radius is the exponent. The value is calibrated by near-wake hot-wire velocity measurement test.

6. The method for predicting the free vortex wake of aerodynamic performance of large wind turbine blades according to claim 5, characterized in that, In step S6, the specific method for updating the circulation distribution of the blade's attached vortex elements using an efficient iterative approach is as follows: Circulation distribution of blade-attached vortex elements at time step n+1 Efficient iterative computation includes the following sub-steps: S61: Given the wake vortex circulation Under the premise of solving the first estimated value of the circulation of the attached vortex element according to the lift line model. and update value , ； ； Where: Lift line solver Input parameters and These represent the circulation values ​​of the blade-attached vortex element and the wake vortex element at the nth time step, respectively. S62, using the first estimate The second solution using the lift line model yields a quadratic estimate of the circulation of the attached vortex element. and update value , ； ； S63, utilizing the nth step attached vortex circulation The two circulation updates yield the attached vortex circulation at the new time n+1 step. The corrected distribution is shown in the following formula. ； In the formula; The relaxation factor, ranging from 0.6 to 1.0, is adjusted through numerical experiments according to the blade's spanwise aerodynamic performance to promote iterative convergence. The attached vortex circulation at the new time n+1 obtained in step S6 The new time-step wake element node coordinates obtained in step S5 Used for calculating the induced velocity at the vortex element node on the wake surface in step S4 at the next time step, and based on the attached vortex element circulation. With the circulation of the wake vortex The distribution of the local ground wash-induced velocity at 1 / 4 chord length of each blade section was calculated. and axial inducing factor and tangential inducing factor The magnitude of the drag coefficient is determined by step S2, and the blade torque M and axial thrust at the new moment are calculated using numerical integration. And the power P value, ； ； ； ； ； ； In the formula: and The coefficients of the normal force and axial force at the j-th cross section are, in order. Let J be the local dynamic pressure at the j-th cross section; Let J be the radial coordinate of the midpoint of the j-th blade segment; Let j be the radial width of the blade segment. This refers to air density.

7. The method for predicting the free vortex wake of aerodynamic performance of large wind turbine blades according to claim 6, characterized in that, In step S7, when performing fully unsteady dynamic aerodynamic performance prediction, the inflow velocity... At each time step update, return to S3 to enter a loop calculation until the predetermined maximum number of time steps is reached. That is, stop the operation.

Images