Lidar-based wind turbine impeller plane online wind speed prediction method

Abstract

The application provides a wind turbine impeller plane wind speed online prediction method based on a laser radar, and comprises the following steps: measuring the wind speed in front of a wind turbine by a four-beam laser radar anemometer arranged in front of the wind turbine in the main direction of wind speed in a wind farm, obtaining a time series of discrete four-point wind speed information of a certain plane in front of the impeller, solving the wind parameters of the four-point wind speed information of the plane, obtaining relevant wind parameters, performing differential transformation on the wind speed engineering model, extracting one-dimensional wind speed time delay information of the dynamic evolution of the radar wind measurement plane to the impeller hub and the rear of the unit, constructing an ultra-short-term plane wind speed model according to the relevant wind parameters, and completing online prediction of the wind speed of each point in the plane and plane equivalent wind speed prediction on the basis of the model. Compared with the traditional Taylor frozen processing, the wind speed prediction accuracy at the impeller of the wind turbine is improved.

Description

TECHNICAL FIELD

[0001] The application belongs to the technical field of automatic control, and particularly relates to a wind turbine impeller plane online wind speed prediction method based on a laser radar. BACKGROUND

[0002] With the development of large-scale wind turbines, the impeller diameter is getting larger and larger, and the whole impeller plane is subjected to more and more obvious unbalanced load due to wind shear, turbulence and tower shadow effect, which seriously affects the stable operation of the unit. The traditional mechanical anemometer is installed on the nacelle, and the wind speed measurement accuracy is greatly affected by the wake, and only the single-point wind speed at the nacelle can be measured, which cannot sense the wind speed change in front of the wind field in advance, is not conducive to the unit to react to the gust in advance, and is easy to cause the unit to overload and run, which endangers the safety of the wind turbine.

[0003] The laser wind measuring radar emits high-speed laser pulses which are reflected by particles in the air, and the time required for the pulse to return is measured by a sensor. Based on the transmission wavelength and the reflection wavelength, the wind speed is determined by using the Doppler effect. The nacelle type wind measuring laser radar has the characteristics of measuring the wind speed in front of the wind turbine in advance, and can calculate the wind speed at the impeller plane in advance according to the wind speed information at the measurement point, thereby providing a basis for the feedforward control and independent variable pitch of the unit. SUMMARY

[0004] The application provides a plane wind speed time delay processing algorithm based on a laser radar, which aims to complete the ultra-short-term plane wind speed online prediction of the wind farm by installing as few nacelle type wind measuring laser radars as possible (i.e. only installing laser radars on the front row of units in the main flow direction), and provides wind speed basis for the feedforward control and independent variable pitch control of the unit. The algorithm takes the line-of-sight wind speed at a certain distance in front of the wind farm measured by the laser radar as the input, solves the wind parameters, and then extracts the one-dimensional wind speed time delay information from the radar wind measuring plane dynamic evolution to the hub and the rear of the unit by differential transformation of the wind speed engineering model. Then, the ultra-short-term plane wind speed model is constructed to complete the online prediction of the unit plane wind speed of the wind farm. Meanwhile, for the unified variable pitch, since the unit is large-scale, the wind speed at the hub height cannot accurately represent the equivalent wind speed received by the impeller plane, so the impeller plane equivalent wind speed calculation method based on the impeller plane wind speed prediction model is also proposed.

[0005] The technical scheme adopted by the application is as follows:

[0006] I. A wind turbine impeller plane wind speed online prediction method based on a laser radar, comprising the following steps:

[0007] S1: The laser radar measures to obtain the time series of line-of-sight wind speeds at four points on the measurement plane;

[0008] The laser transmitter of the four-beam laser radar emits four laser beams in turn to the front at a very high frequency, and the time series of the discrete line-of-sight wind speed of four points on the measurement plane can be obtained.

[0009] The measurement plane in S1 is the measurement plane of the laser radar, located in front of the impeller plane, at a distance of 50-200 m.

[0010] S2: Establish an inertial coordinate system and a wind speed coordinate system with the laser radar installation position as the origin, and obtain four wind parameters: the average flow direction wind speed v H , the vertical wind shear index α, the horizontal incident angle The vertical incident angle The coordinates of the four-beam measurement points in the inertial coordinate system are determined by the measurement plane distance and the beam angle.

[0011] S2 specifically is:

[0012] 2.1) Obtain the line-of-sight wind speed v i,los of the four-beam measurement points by laser radar measurement;

[0013]

[0014] In formula (1), (x i,L , y i,L , z i,L ) represents the coordinates of the i-th measurement point in the inertial coordinate system, i=0,1,2,3;

[0015]

[0016] Wherein, L represents the horizontal distance between the laser radar and the measurement plane; θ1 represents the horizontal angle of the beam; θ2 represents the vertical angle of the beam;

[0017] In formula (1), l i represents the distance between the i-th measurement point and the origin of the laser radar, u i,L , v i,L , w i,L represent the wind speed components in the inertial coordinate system; the measurement point is a point on the measurement plane of the laser radar.

[0018] 2.2) Establish the transformation relationship between the wind speed coordinate system and the inertial coordinate system:

[0019]

[0020] In the formula, u i,W , v i,W , w i,W represent the wind speed components in the wind speed coordinate system, the horizontal incident angle of wind; the vertical incident angle of wind;

[0021] 2.3) The wind speed expression at height h is obtained according to the vertical direction wind shear exponent law:

[0022]

[0023] wherein, v H represents the average flow direction wind speed at the hub height of the impeller; H represents the hub height of the impeller; a represents the vertical direction wind shear exponent; z represents the vertical distance between the target position on the impeller plane and the center of the hub; h represents the height of the target position on the impeller plane from the horizontal ground;

[0024] The Taylor expansion is performed on equation (3), and the terms of 3 orders and above are ignored to obtain:

[0025]

[0026] 2.4) In the wind speed coordinate system, the average wind direction wind speed v H at the hub height is taken as the reference, the vertical direction wind shear is considered, the horizontal direction wind shear is ignored, and the wind speed vector is represented as:

[0027]

[0028] 2.5) Equation (5) is substituted into equation (2) to obtain:

[0029]

[0030] The above equation is substituted into equation (1) to obtain four wind parameters, which are respectively: the average flow direction wind speed v H at the hub height, the vertical direction wind shear exponent a, the horizontal incident angle the vertical incident angle

[0031] wherein,

[0032] S3: Time delay processing algorithm: time delay processing is performed on the wind speed engineering model to extract the wind speed size and time delay information of the first unit and the units behind the unit in the stringer unit;

[0033] The time delay processing algorithm of the step S3 is:

[0034] The variable time t is introduced on the basis of the traditional wind speed engineering model, that is,

[0035] The traditional wind speed engineering model is:

[0036]

[0037] where x represents the position coordinate along the wind flow direction; v x represents the size of the wind speed at position x along the x direction of the inertial coordinate system, which is the direction perpendicular to the plane of the turbine unit;

[0038] v in represents the axial wind speed perpendicular to the plane of the impeller; for the wind evolution model, v in is obtained by the laser radar, and for the wake model, v in is the upstream unit inflow wind speed; f(x) is an expression of the relationship between the wind speed and the position, and the form of f(x) is different for different engineering models;

[0039] Substituting the traditional wind speed engineering model, we can obtain:

[0040]

[0041] Integrating both sides of equation (7) with respect to t, we obtain:

[0042] v in t = F(x) + C (8)

[0043] where: F(x) is the primitive function of f(x); C is a constant, which can be determined by the boundary conditions in specific engineering models; let the target position on the plane of the impeller be x t , then the delay time of the wind from the measuring point to the target position is where x t = 0 at the impeller of the unit; thus the wind speed v t at the target position x t and the time delay information t t are obtained:

[0044]

[0045] where v t is the wind speed at the target position x t after time t t .

[0046] S4: processing the average flow direction wind speed of the measuring plane according to the time delay processing algorithm, extracting the one-dimensional wind speed time delay information from the dynamic evolution of the radar measuring plane to the hub of the impeller and the rear of the unit, and further obtaining the wind speed time sequence;

[0047] The step S4 is specifically:

[0048] At the position of the hub of the unit corresponding to the measuring plane, the axial wind speed perpendicular to the plane of the impeller is represented as

[0049] ​

[0050] The axial wind speed v H,r is processed by using the time delay processing algorithm of step S3 to obtain one-dimensional time delay information of the wind in the propagation process, and further obtain a wind speed time sequence.

[0051] S5: Determine whether the one-dimensional wind speed time delay information has wind speed aliasing:

[0052] When the wind speed time sequence processed by the time delay method has two or more wind speeds at the same time, it is determined that wind speed aliasing occurs, such as when a gust occurs, which is a sudden increase in wind speed, and the target position will have multiple wind speed data at the same time.

[0053] If wind speed aliasing does not occur, output the wind speed time sequence of step S4;

[0054] If wind speed aliasing occurs, use the energy unification method to unify the multiple axial wind speeds and wind shear exponents that occur at the same time to obtain a one-dimensional wind speed time sequence at the hub height of the unit, specifically:

[0055] 5.1) Unify multiple axial wind speeds occurring at the same time

[0056] When wind speed aliasing occurs, after time delay processing, different axial velocities of atmospheric particles v1, v2,..., v n occur at the same time on the target plane; assuming that the number of particles of each velocity is equal, it is counted as unit 1, and the unified wind speed is v e ; thus obtaining the energy equation:

[0057]

[0058] In the formula, n represents different sizes of axial velocities;

[0059] The unified wind speed v e is expressed as:

[0060]

[0061] 5.2) Unify the vertical wind shear exponents of the measurement plane and the target plane

[0062] According to the exponential law wind shear formula, the wind shear speeds v z1 , v z2 ……v zn corresponding to different times at the position corresponding to the measurement plane height z of the impeller plane are obtained as:

[0063]

[0064] In the formula, α nTo measure the vertical wind shear index of the plane at different time;

[0065] According to formula (11), (12), the wind speed v of the impeller plane height z at time t after energy unification method processing is obtained zt The central wind speed v at the hub is t The ratio is:

[0066]

[0067] In the formula, α e Indicates the unified wind shear index;

[0068] Taylor expansion is performed on formula (13), and according to the first order coefficient, the following is obtained:

[0069]

[0070] Therefore, when the wind speed aliasing phenomenon occurs, the unified axial wind speed and the wind shear index are:

[0071]

[0072] S6, constructing an ultra-short-term plane wind speed model according to the related wind parameters, and completing online prediction of the wind speed of each point of the unit plane based on the constructed model;

[0073] The process of constructing the ultra-short-term plane wind speed model in step S6 is specifically:

[0074] When the wind speed aliasing phenomenon occurs, α appearing in the process of constructing the ultra-short-term plane wind speed model is the unified wind shear index α e , and v h is the unified axial wind speed v e ;

[0075] For the vertical wind shear, let z = rcosθ, and formula (4) is expressed in polar coordinate form as:

[0076]

[0077] In the formula, v h The axial wind speed v H,r of the hub height of the measurement plane is obtained by time delay processing algorithm; r is the distance from the target point on the impeller plane to the hub center; and θ is the azimuth angle;

[0078] Then the wind shear component is expressed as:

[0079]

[0080] For the tower shadow effect component:

[0081]

[0082] wherein Δx represents the distance from the impeller plane to the tower centerline, in m; b represents the hub radius of the machine, in m; R represents the impeller radius, in m;

[0083] The wake component of the 2D_K Jensen wake model is:

[0084]

[0085] wherein, a represents the axial induction factor of the machine.

[0086] wherein v0 represents the incoming wind speed (axial wind speed perpendicular to the impeller plane) of the previous machine; r x represents the wake radius at x downstream of the machine, r x = k wake x + R, k wake represents the wake expansion coefficient; for the first machine, v r = 0.

[0087] Considering the wind shear component, the tower shadow component, and the wake effect component, the ultra-short-term planar wind speed model is obtained:

[0088]

[0089] wherein, is the result after time delay processing of the wind evolution model, is the result after time delay processing of the wake model;

[0090] wherein v r,θ represents the wind speed at the corresponding polar coordinate on the impeller plane; C1 and C2 are constants.

[0091] S7. Calculate the equivalent wind speed on the impeller plane according to the predicted planar wind speed at each point.

[0092] The step S7 is specifically:

[0093] The equivalent wind speed calculation method is derived from the perspective of energy conservation:

[0094] Divide the impeller plane into n equal parts from top to bottom, and take the wind speed v i at the center of each part as the representative wind speed:

[0095]

[0096] wherein v eq represents the equivalent wind speed, in m / s; m eq represents the equivalent mass of air, m eq = ρAveq where p is the air density, unit kg / m 3 , A is the impeller area, unit m 2 ); m i represents the mass of the i-th part of air, m i = pA i v i , where A i is the i-th part area;

[0097] The impeller plane equivalent wind speed is represented as:

[0098]

[0099] In order to make the calculated equivalent wind speed more accurately represent the entire plane, the impeller plane is divided into an infinite number of small blocks, and the representative wind speed v i is obtained from the plane wind speed prediction model. Since the wind speed mathematical model obtained by S7 is relatively complex, under the condition of ensuring the same dimension, the integral of the simplified formula (22) is carried out:

[0100]

[0101] Thus, the equivalent wind speed of the impeller plane is obtained.

[0102] The beneficial effects of the present application are:

[0103] The present application provides a method for accurately sensing the wind speed information of the impeller plane of the unit in the future period of time in advance, which provides wind speed basis for unit feedforward control and independent variable pitch control. The plane time delay processing algorithm based on laser radar can predict the wind speed of the unit plane. The equivalent wind speed of the impeller plane obtained based on this algorithm has high correlation with the equivalent wind speed obtained by the unit operation data. Compared with the traditional Taylor frozen processing, the time delay processing method improves the wind speed prediction accuracy of the wind turbine impeller by considering the wind speed evolution process. BRIEF DESCRIPTION OF DRAWINGS

[0104] Figure 1 A time delay processing algorithm flowchart provided for examples of the present application;

[0105] Figure 2 A wind speed aliasing processing schematic diagram based on the energy unification method;

[0106] Figure 3 a, 3b are the results analysis of the time delay processing method and the Taylor frozen method;

[0107] Figure 4 The comparison of the predicted equivalent wind speed of the impeller plane, the predicted wind speed at the hub height and the unit backstepping equivalent wind speed;

[0108] Figure 5A schematic diagram of a line-of-sight beam of a four-beam lidar and a related coordinate system;

[0109] Figure 6 A wind speed propagation route map;

[0110] Figure 7 A wind speed component schematic diagram;

[0111] Figure 8 A lidar field measurement schematic diagram;

[0112] Figure 9 A prediction error distribution of the time delay processing method in the high wind speed section and the low wind speed section. Embodiments

[0113] The specific embodiments of the application will be further described in conjunction with the accompanying drawings, which are used to provide a further understanding of the application, form a part of this application, and do not constitute a limitation of the application.

[0114] The prediction method of the application is a plane wind speed time delay processing algorithm based on a four-beam lidar, which considers the aerodynamic interference between the tandem unit, establishes the plane wind speed propagation relationship from the upstream of the first wind turbine unit to the downstream, and extracts the quantitative relationship between the wind speed and the propagation time on the basis of the original engineering model. In addition, the equivalent wind speed calculation formula corresponding to the plane wind speed model is also proposed, which aims to replace the wind speed at the hub height with the equivalent wind speed of the impeller, to provide a more reliable wind speed basis for unified variable pitch, to reduce the blade and hub load, and to prolong the service life of the wind turbine.

[0115] The method of the application specifically comprises the following steps:

[0116] S1, the lidar measures the time series of line-of-sight wind speed of four points on the same plane at a certain distance in front of the impeller plane.

[0117] The laser emitter of the four-beam lidar emits four laser beams in turn at a very high frequency to the front, and the discrete wind speed time series of four points on the set distance plane can be obtained.

[0118] The line-of-sight wind speed v i,los of the four-beam measurement point is measured by the lidar

[0119]

[0120] In the formula, (x i,L ,y i,L ,z i,L ) represents the coordinates of the i-th measurement point in the inertial coordinate system, i=0,1,2,3

[0121] Where, l irepresents the distance of the measuring point from the origin of the lidar, u i,L , v i,L , w i,L represents the wind speed component in the inertial coordinate system;

[0122] S2, an inertial coordinate system and a wind speed coordinate system are established with the lidar installation position as the origin, and based on the transformation relationship between the inertial coordinate system and the wind speed coordinate system, four wind parameters are obtained: the average flow direction wind speed v H at the hub height, the vertical direction wind shear index a, the horizontal incidence angle the vertical incidence angle

[0123] In the step S2:

[0124] An inertial coordinate system and a wind speed coordinate system are established and related, and the coordinates of the four-beam measuring points in the inertial coordinate system are determined by the measuring plane distance and the beam angle, Figure 5 a is the schematic diagram of the line-of-sight beam of the four-beam lidar. As shown in Figure 5 b, an inertial coordinate system (x L , y L , z L ) is established with the lidar installation position as the origin, wherein the position coordinates of the four measuring points can be represented as:

[0125]

[0126] In the formula: i represents the four measuring points on the measuring plane, i = 0, 1, 2, 3; j represents the ten measuring distances of the lidar, j = 0-9; L j represents the horizontal distance of the lidar from the measuring plane; θ1 represents the horizontal angle of the beam; θ2 represents the vertical angle of the beam. For the lidar, the wind speed vector along the direction of the emitted beam is called the line-of-sight wind speed (v LOS ).

[0127] In order to better represent the related wind parameters, a wind speed coordinate system (x Figure 5 , y W , z W , z W ) is established as shown in

[0128]

[0129] In the formula: u i,L , v i,L , w i,L represents the wind speed component in the inertial coordinate system; u i,W、v i,W 、w i,W denotes the wind speed component in the wind speed coordinate system; denotes the horizontal incidence angle; denotes the vertical incidence angle.

[0130] Let:

[0131]

[0132] Thus:

[0133]

[0134] For the natural wind in the near-surface layer, the wind speed is affected by the low-layer atmospheric stratification and the ground roughness, and will increase with the increase of height. Generally, the change law can be represented by the exponential law, that is:

[0135]

[0136] wherein v H denotes the average flow direction wind speed at the hub height, m / s; H denotes the hub height, m; and a denotes the vertical direction wind shear exponent. Taylor expansion is performed on equation (6), and terms of 3 orders and above are ignored to obtain:

[0137]

[0138] In the wind speed coordinate system, taking the average wind direction wind speed v H at the hub height as the reference, the vertical direction wind shear is considered, and the horizontal direction wind shear is ignored. The wind speed vector can be expressed as:

[0139]

[0140] For the convenience of writing, let:

[0141]

[0142] Simultaneous equations (4), (5), (8), and (9) can be obtained:

[0143]

[0144] The line-of-sight wind speed of the laser radar can be written as:

[0145]

[0146] wherein l i denotes the distance between the measurement point and the origin of the laser radar, unit m,

[0147] Simultaneous equations (10) and (11) can be obtained:

[0148]

[0149] Four-beam lidar returns four line-of-sight wind speeds at four measuring points, and four equations in the form of equation (12) can be solved simultaneously for four unknown variables The equation group has a solution, and the relevant wind parameters required to modify the wind speed model of the impeller plane can be obtained.

[0150] S3, time delay processing algorithm: time delay processing is performed on the wind speed engineering model to extract the wind speed and time delay information of the first unit and the units behind it in the string unit;

[0151] S3 specifically includes: in a wind farm, the natural wind is affected by wind turbines during its forward movement, and a region where the wind speed gradually decays appears in front of the units. In this region, the wind evolution model can be used to represent the wind speed information; after passing through the units, a wake is formed behind the impeller, and the wake model can be used to represent the wind speed information in the wake region, as shown in Figure 6 The commonly used wind speed engineering model only describes the relationship between the wind speed and the position during the propagation of the wind in the wind farm, and lacks the description of the relationship between the wind speed and the propagation time:

[0152]

[0153] In the formula: x represents the position coordinate of the wind flow in the upward direction; v represents the size of the wind speed at position x in the x direction, m / s; v in represents the size of the incoming flow wind speed, m / s; for the wind evolution model, v in can be obtained by lidar (v H,r ), and for the wake model, v in takes the incoming flow wind speed of the upstream unit, and f(x) has different forms for different engineering models.

[0154] In the actual propagation process, the wind will continuously evolve, and the traditional Taylor freezing method ignores this evolution process, considering that the wind moves forward in the Taylor frozen form, and the propagation time is calculated by dividing the propagation distance by the average wind speed. The time delay processing method is proposed, which considers the wind speed evolution process and replaces the average wind speed with the real-time wind speed. Through differential transformation of the wind speed engineering model, the time information of the wind on the propagation path is calculated. The time delay processing algorithm introduces the variable t on the basis of the traditional wind speed engineering model, that is Substituting equation (13) into equation (13) gives

[0155]

[0156] Then integrate both sides of the above equation with respect to t to get

[0157] v in t=F(x)+C (15)

[0158] In the formula: C can be determined by boundary conditions in a specific engineering model. Let the target location be x. t The time delay for the wind to travel from the measurement point to the target location is... x at the impeller of the unit t =0.

[0159] Thus, information about wind speed with respect to time delay has been extracted based on the engineering model. Target location x t The wind speed and time delay information at that location can be represented as:

[0160]

[0161] In the formula: v t —After time t t Wind speed at the target location, m / s; f(x) t ) and F(x t The form varies depending on the engineering model.

[0162] S4. The average directional wind speed at the hub height of the measuring plane unit is processed according to the time delay processing algorithm to extract the one-dimensional wind speed time delay information of the radar measuring plane dynamic evolution to the impeller hub and the rear of the unit.

[0163] S4 specifically includes:

[0164] At the position corresponding to the hub height on the measuring plane, the axial wind speed perpendicular to the impeller plane is expressed as:

[0165]

[0166] For v H,r One-dimensional time delay information of wind propagation is obtained using a time delay processing algorithm.

[0167] S5. Determine whether wind speed aliasing occurs in the one-dimensional wind speed time delay information;

[0168] S5 specifically includes:

[0169] The wind speed time series processed by the time delay method has two or more wind speeds appearing at the same time, indicating that wind speed overlap has occurred. For example, when there is a sudden increase in wind speed due to gusts, the target location may have multiple wind speed data at the same time.

[0170] S6. If wind speed overlap occurs, the energy unification method is used to unify multiple axial wind speeds and wind shear indices that occur at the same time to obtain a one-dimensional wind speed time series at the hub height of the unit.

[0171] S6 specifically includes: if the wind speed aliasing phenomenon occurs, that is, after time delay processing, the target plane at the same time appears different axial velocity of atmospheric particles v1, v2…v n , and assuming that the number of particles of each speed is equal, it is 1, and the unified wind speed is v e

[0172] According to the energy equation

[0173]

[0174] n represents different sizes of axial velocity

[0175] Therefore, the unified wind speed can be represented as:

[0176]

[0177] Similarly, there are differences in the vertical wind shear index of the measurement plane and the target plane, and the following method is used for superposition:

[0178] Assuming that at a certain time t, n different sizes of axial velocity of atmospheric particles v h1 , v h2 …v hn appear at the hub center at the same time, then according to the exponential law wind shear formula, the wind shear wind speed v z1 , v z2 …v zn corresponding to each wind speed at the measurement plane height z of the impeller plane can be obtained.

[0179] At the measurement plane corresponding to the impeller plane height z, the wind shear wind speed v z1 , v z2 …v zn corresponding to different times is:

[0180]

[0181] α n is the vertical wind shear index of the measurement plane at different times

[0182] According to formulas (19) and (20), the ratio of the wind speed v zt at the impeller plane height z at time t to the hub center wind speed v t after energy unification processing can be obtained:

[0183]

[0184] Perform Taylor expansion on formula (21), and according to the first order coefficient, we can get:

[0185]

[0186] So when the wind speed aliasing phenomenon occurs, the unified axial wind speed and wind shear exponent are:

[0187]

[0188] S7, according to the relevant wind parameters to build a short-term plane wind speed model, on this basis to complete the unit plane each point wind speed online prediction;

[0189] S7 specifically includes: for the first row of units in the wind farm, the uneven distribution of the blade plane wind speed in space mainly needs to consider the vertical wind shear and tower shadow effect, for the rear row of units, the wind speed of the blade plane also needs to consider the wake effect generated by the unit, Figure 7 The schematic diagram of each wind speed component is shown in Fig. 1.

[0190] For vertical wind shear, let z = rcosθ, r is the distance from the target point on the blade plane to the hub center; θ is the azimuth angle. Equation (7) is expressed in polar form as:

[0191]

[0192] In the formula: v h The axial wind speed v H,r is obtained by time delay processing algorithm processing;

[0193] The axial wind speed at the hub center of the unit is obtained. For the first row of units, v h The axial wind speed v H,r is obtained by time delay processing algorithm processing after time delay processing of the wind evolution model; for the rear row of units, v h The wind shear component can be expressed as:

[0194]

[0195] For the tower shadow effect component:

[0196]

[0197] In the formula: Δx is the distance from the blade plane to the tower centerline, m; b is the hub radius of the unit, m; (where R is the radius of the blade, m).

[0198] 2D_K Jensen wake model, the plane expansion of the wake component is:

[0199]

[0200] wherein: v0 is the incoming wind speed of the previous unit; r x represents the wake radius at x downstream of the unit, r x = k wake x + R.

[0201] Considering the wind shear component, the tower shadow component and the wake effect component, the wind speed at the impeller plane of the unit can be written as:

[0202] v r,θ = v h + v ws + v ts + v r (28)

[0203] S8, calculating the equivalent wind speed at the impeller plane according to the predicted wind speed at each point of the plane.

[0204] S8 specifically comprises:

[0205] dividing the impeller plane into n equal parts from top to bottom, taking the wind speed v i at the center of each part as the representative wind speed:

[0206]

[0207] wherein: v eq — equivalent wind speed, m / s; m eq — equivalent mass of air, m eq = pAv eq (wherein p is the air density, kg / m 3 ; A is the impeller area, m 2 ); m i — the mass of air in the i-th part, m i = pA i v i (wherein A i is the area of the i-th part).

[0208] Then the equivalent wind speed at the impeller plane can be written as:

[0209]

[0210] In order to make the calculated equivalent wind speed more accurately represent the entire plane, the impeller plane is divided into an infinite number of small blocks, and the representative wind speed v i of each small block can be obtained from the plane wind speed prediction model, and finally calculated in the form of integration according to formula (30).

[0211] Since the wind speed mathematical model obtained by S7 is relatively complex, the simplified formula (30) is integrated under the condition of ensuring the same dimension:

[0212]

[0213] DETAILED DESCRIPTION: WIND FARM MEASURED WIND SPEED DATA VERIFICATION

[0214] 1) Time delay processing method verification

[0215] In order to verify the one-dimensional wind speed time delay processing method based on laser radar, in the measured data of the laser wind measuring radar in a certain wind farm, the laser radar measurement data at 250m and 160m in front of the unit were selected, Figure 8 The schematic diagram for laser radar field measurement.

[0216] In order to test the prediction performance of the model at different wind speed segments, this paper takes 250m as the measurement point and 160m as the target point, and divides the wind speed data at 250m into high and low wind speed segments. Table 1 shows part of the information of the test unit.

[0217] Table 1 Part of the information of the test unit

[0218]

[0219] For the wind evolution model, this paper selects the wind evolution model proposed by Medici et al. in 2011, which describes the attenuation process of the wind in front of the wind turbine during propagation, and the relationship between wind speed and position is expressed as:

[0220]

[0221] In the formula: v ∞ represents the wind speed at infinity, m / s; a represents the axial induction coefficient of the wind turbine.

[0222] Therefore:

[0223]

[0224] Since f(x t ) itself is relatively complex, there is no analytical solution for its original function, so this paper selects the numerical solution of the original function of f(x t ) instead of F(x t ). For the 160m in front of the unit, according to the time delay processing algorithm, we have:

[0225]

[0226] After calculation, the wind speed curve at 160m in front of the unit predicted by the time delay processing algorithm is obtained, and compared with the wind speed curve at 160m in front of the unit measured by the laser radar, Figure 3 a are the prediction wind speed curve comparison graphs of high and low wind speed segments respectively, Figure 3b are the correlation analysis of the predicted wind speed by the time delay processing method and the Taylor frozen method and the measured wind speed by the laser radar in the high wind speed section and the low wind speed section.

[0227] It is calculated that the average absolute error of the wind evolution model corrected by the time delay processing method in the high and low wind speed sections is 0.19 m / s and 0.18 m / s, and the average absolute percentage error is 2.4% and 3.6%; the average absolute error of the wind evolution model processed by the Taylor frozen method in the high and low wind speed sections is 0.23 m / s and 0.21 m / s, and the average absolute percentage error is 2.8% and 2.6% respectively. Compared with the Taylor frozen method, the time delay processing method improves the average absolute error by 17% and 14% respectively in the prediction of high and low wind speed at 160 m in front of the impeller. From the prediction trend, the calculated wind speed at 160 m by the wind evolution model corrected by the time delay processing method and the wind evolution model using the Taylor frozen method basically conforms to the measured wind speed by the laser radar anemometer; but from the prediction error, the wind evolution model corrected by the time delay processing method has smaller error range and smaller average error than the wind evolution model using the Taylor frozen method. Figure 3 It can be seen that the wind evolution model processed by the Taylor frozen method lags behind the measured wind speed by the laser radar, which is due to the fact that the wind speed used by the Taylor frozen method in the calculation of time is the average wind speed in a certain period of time, which has higher similarity requirement for wind speed data and cannot completely represent the wind speed at a certain moment, so compared with the real wind speed propagation, there will be the phenomenon of advance or lag in the propagation time. From the prediction trend, the wind evolution model processed by the Taylor frozen method basically conforms to the measured wind speed by the laser radar anemometer, but from the prediction error, the wind evolution model processed by the Taylor frozen method has larger error range and larger average error than the wind evolution model using the Taylor frozen method. Figure 3 It can be seen that compared with the Taylor frozen method, the predicted wind speed by the time delay processing method shows higher similarity with the measured wind speed by the laser radar.

[0228] Figure 9 The prediction error distribution of the time delay processing method and the Taylor frozen method at 160 m in the high wind speed section and the low wind speed section can be seen from the prediction error. For the wind evolution model corrected by the time delay processing method, the error range of 90% of the data in the high and low wind speed sections is-0.53-0.24 m / s and-0.63-0.21 m / s, and the average error is-0.12 m / s and 0.16 m / s respectively. The numerical value of the wind evolution model processed by the Taylor frozen method is slightly larger than that of the time delay processing method. From the above analysis, it can be seen that compared with the Taylor frozen method, the wind evolution model corrected by the time delay processing method has improved in the prediction of wind arrival time, wind speed error and error range.

[0229] 2) Plane equivalent wind speed verification

[0230] In order to compare with the real equivalent wind speed, the real equivalent wind speed is back calculated by using the operation data of the unit. For the mechanical anemometer, the wind speed can be back calculated according to the rotation speed of the wind cup, and the wind turbine also has similar aerodynamic characteristics with the mechanical anemometer, so the algorithm based on aerodynamic torque, pitch angle and tip speed ratio can also be used to estimate the equivalent wind speed of the wind wheel.

[0231] The measured data of a four-beam laser radar in a certain wind field are used for verification. Two groups of wind speed data at different time periods on a certain day are selected, and the predicted wind speed at the hub height, the predicted equivalent wind speed at the impeller plane, and the equivalent wind speed of the unit back-propagation are verified, Figure 4 The verification results of time periods 1 and 2 are shown in Figs. 6 and 7, respectively.

[0232] From Figure 4 It can be seen that the single-point predicted wind speed at the hub height and the equivalent wind speed at the impeller plane of the modified wind speed model have strong correlation with the equivalent wind speed of the unit back-propagation. However, compared with the equivalent wind speed of the unit back-propagation, the single-point predicted wind speed curve at the hub height is shifted upward as a whole, because the single-point wind speed at the hub height does not consider the influence of tower shadow effect and vertical wind shear. The equivalent wind speed at the impeller plane calculated by the modified wind speed model can better represent the actual equivalent wind speed during the operation of the unit, because the modified wind speed model considers the influence of tower shadow effect and vertical wind shear.

[0233] Compared with the equivalent wind speed of the unit back-propagation, the average absolute error of the single-point predicted wind speed at the hub height and the equivalent wind speed at the impeller plane in time period 1 is 0.92 m / s and 0.68 m / s, respectively, and the average absolute percentage error is 12% and 9%, respectively. The average absolute error of the wind speed prediction in time period 2 is 1.05 m / s and 0.76 m / s, respectively, and the average absolute percentage error is 15% and 11%, respectively. For time periods 1 and 2, compared with the single-point predicted wind speed at the hub height, the average absolute percentage error of the equivalent wind speed at the impeller plane of the modified wind speed model is reduced by 25% and 27%, respectively. Therefore, compared with the single-point wind speed at the hub height, the predicted equivalent wind speed proposed in the present application can better represent the actual equivalent wind speed during the operation of the unit, which also verifies the effectiveness of the impeller plane wind speed prediction method proposed in the present application.

Claims

1. A wind turbine impeller plane wind speed online prediction method based on laser radar, characterized in that, Comprise the following steps: S1: the laser radar is obtained by measuring the time series of the line-of-sight wind speed of four point positions on the measurement plane; S2: Establish an inertial coordinate system and a wind speed coordinate system with the laser radar installation position as the origin, and based on the transformation relationship of the inertial coordinate system and the wind speed coordinate system, four wind parameters are obtained: the average flow direction wind speed v H at the hub height, the vertical direction wind shear index a, the horizontal incidence angle the vertical incidence angle S3: time delay processing algorithm: the wind speed engineering model is processed by time delay, and the wind speed and time delay information of the first unit and the unit behind the stringer unit are extracted; S4: according to the time delay processing algorithm, the average flow direction wind speed of the measurement plane is processed, the one-dimensional wind speed time delay information of the dynamic evolution of the radar measurement plane to the hub of the impeller and the unit behind is extracted, and then the wind speed time sequence is obtained; S5: judge whether the one-dimensional wind speed time delay information occurs wind speed aliasing: If wind speed aliasing does not occur, output the wind speed time sequence of step S4; If wind speed aliasing occurs, use energy uniform method to unify the multiple axial wind speed and wind shear index appearing at the same time to obtain the one-dimensional wind speed time sequence at the hub height of the unit; S6, according to the related wind parameters, the ultra short term plane wind speed model is constructed, and the wind speed of each point on the unit plane is predicted on line based on the constructed model; S7, according to the predicted plane wind speed of each point, the equivalent wind speed of the impeller plane is calculated.

2. The LIDAR based wind turbine impeller plane wind speed online prediction method according to claim 1, characterized in that, The measurement plane in S1 is the measurement plane of the laser radar, which is located in front of the impeller plane, and the distance is 50-200m.

3. The LIDAR based wind turbine impeller plane wind speed online prediction method according to claim 1, characterized in that, S2 is specifically: 2.1) The line-of-sight wind speed v is measured by the laser radar at four measuring points i,los ; In formula (1): (x i,L ,y i,L ,z i,L ) represents the coordinates of the i-th measuring point in the inertial coordinate system, i = 0, 1, 2, 3; Wherein, L represents the horizontal distance between the laser radar and the measurement plane; θ1 represents the horizontal angle of the beam; θ2 represents the vertical angle of the beam; In formula (1), l i represents the distance between the i-th measuring point and the origin of the laser radar, u i,L , v i,L , w i,L represents the wind speed component in the inertial coordinate system; 2.2) construct the transformation relationship between the wind speed coordinate system and the inertial coordinate system: wherein u i,W , v i,W , w i,W denote the wind speed components in the wind speed coordinate system, denotes the horizontal angle of incidence of the wind; denotes the vertical angle of incidence of the wind; 2.3) according to the vertical direction wind shear index law, the wind speed expression at height h is obtained: where v H represents the average flow direction wind speed at the hub height of the impeller; H represents the hub height of the impeller; a represents the vertical direction wind shear index; z represents the vertical distance from the center of the hub to the target position on the impeller plane; h represents the height from the ground to the target position on the impeller plane; Taylor expansion is carried out on formula (3), and items higher than 3 are ignored: 2.4) in the wind speed coordinate system, with the average wind direction wind speed v at hub height H For reference, considering the vertical wind shear and ignoring the horizontal wind shear, the wind speed vector is represented as: 2.5) Simultaneous equations (5), (2), (1) are solved to obtain four wind parameters, respectively: the average flow direction wind speed v at the hub height H , the vertical direction wind shear exponent α, the horizontal incident angle the vertical incident angle 4. The LIDAR based wind turbine impeller plane wind speed online prediction method according to claim 1, characterized in that, The time delay processing algorithm of S3 is: On the basis of traditional wind speed engineering model, variable time t is introduced, i.e. The traditional wind speed engineering model is: where x represents the position coordinate in the flow direction; v x represents the magnitude of the wind speed at position x in the x direction of the inertial coordinate system, which is the direction perpendicular to the impeller set plane; v in indicates the axial wind speed perpendicular to the plane of the rotor: for the wind evolution model, v in obtained by the lidar, for the wake model, v in takes the upstream turbine inflow wind speed; f(x) is the expression of the wind speed and position, and the form of f(x) is different for different engineering models; Substituting into the traditional wind speed engineering model, we have: Substituting into the traditional wind speed engineering model, we have: Integrate both sides of formula (7) with respect to t, and get: v in t = F(x) + C (8) In the formula: F(x) is the original function of f(x); C is a constant; Let the target position on the impeller plane be x t The delay time of the wind from the measuring point to the target position is Where x t is the target position on the impeller plane t = 0; so as to obtain the wind speed size v t and the time delay information t t at the target position x where v t is the wind speed at the target position x t after time t t .

5. The LIDAR based online prediction of wind speed at the plane of the rotor of a wind turbine method according to claim 1, characterized in that, S4 is specifically: At the position of the unit hub corresponding to the measurement plane, the axial wind speed perpendicular to the impeller plane is represented as The time delay processing algorithm of step S3 is used to process the axial wind speed v H,r and obtain one-dimensional time delay information of the wind in the propagation process, and further obtain the wind speed time series.

6. The LIDAR based wind turbine impeller plane wind speed online prediction method according to claim 1, characterized in that, In S5, the method for judging whether the one-dimensional wind speed time delay information occurs wind speed aliasing is: When the wind speed time sequence processed by the time delay method exists two or more wind speeds at the same time, it is determined that wind speed aliasing occurs.

7. The LIDAR based wind turbine impeller plane wind speed online prediction method according to claim 1, characterized in that, In S5, the process of using energy uniform method to unify the multiple axial wind speed and wind shear index appearing at the same time is specifically: 5.1) unify the multiple axial wind speed appearing at the same time When the wind speed aliasing phenomenon occurs, after time delay processing, the target plane appears different axial velocity atmospheric particles v1, v2…v n ; assuming that the number of particles of each speed is equal, it is counted as unit 1, and the unified wind speed is v e ; thus the energy equation is obtained: In the formula, n represents the different size of axial velocity; Unified wind speed v e is expressed as: 5.2) unify the vertical wind shear index of the measurement plane and the target plane According to the exponential law wind shear formula, the wind shear speed v corresponding to different time at the position corresponding to the measuring plane of the impeller plane height z is obtained z1 、 z2 … zn v In the formula, α n is the vertical wind shear index of the measurement plane at different times; According to the formula (11), (12), the wind speed v at the height z of the impeller plane at the time t after the energy unification method is processed zt The ratio is: t The ratio is: In the formula, α e denotes the unified wind shear index; Taylor expansion is carried out on formula (13), and according to the first order coefficient, the following formula is obtained: Therefore, when wind speed aliasing occurs, the unified axial wind speed and wind shear index are:

8. The online prediction method of wind speed at the plane of the impeller of a wind turbine based on lidar according to claim 2, characterized in that, The process of constructing the ultra short term plane wind speed model in S6 is specifically: For vertical wind shear, let z=rcosθ, and formula (4) is expressed in polar coordinate form as: where: v h The axial wind speed v at the hub height of the plane is measured by the four-beam laser radar H,r The time delay processing algorithm is processed; r is the distance from the target point on the impeller plane to the hub center; θ is the azimuth angle; When wind speed aliasing occurs, where a is the unified wind shear exponent a e , where v h is the unified axial wind speed v e ; Then the wind shear component is represented as: For tower shadow effect component: wherein: Δx represents the distance from the impeller plane to the tower midline; b represents the machine hub radius; R is the impeller radius; The wake component of 2D_K Jensen wake model is: wherein, a represents the axial induction factor of the machine set; where v0represents the incoming wind speed of the previous turbine; r x represents the wake radius downstream of the turbine at x, r x = k wake x + R, k wake represents the wake expansion coefficient; for the first turbine, v r = 0; Considering the wind shear component, tower shadow effect component and wake effect component, the ultra short term plane wind speed model is obtained: wherein, is the wind evolution model time delay processed result, is the wake model time delay processed result; where v r,θ represents the wind speed at the corresponding polar coordinate on the impeller plane; C1 and C2 are constants.

9. The LIDAR based wind turbine impeller plane wind speed online prediction method according to claim 1, characterized in that, S7 is specifically: Divide the impeller plane into n equal parts from top to bottom, and take the wind speed v at the center of each part i Representative wind speed: where v eq represents the equivalent wind speed; m eq represents the equivalent mass of air, m eq = pAv eq where p is the air density and A is the impeller area; m i represents the mass of air in the i-th portion, m i = pA i v i where A i is the i-th portion area; Then the equivalent wind speed of the impeller plane is represented as: Divide the impeller plane into infinite small blocks, each block represents the wind speed v i From the plane wind speed prediction model, integrate the simplified equation (22) under the condition of ensuring the same dimension: Thus, the equivalent wind speed of the impeller plane is obtained. Thus, the equivalent wind speed of the impeller plane is

Images