Reduced order modeling method for on-line state monitoring of blades of wind turbine generator

Abstract

The invention discloses a reduced order modeling method for online state monitoring of a wind turbine generator blade, which comprises the following steps of: 1) constructing a kinetic equation of the wind turbine generator blade by using a central rigid body-flexible beam model, and dispersing the kinetic equation by using a finite element method to obtain a finite element model; (2) solving the finite element model obtained in the step (1) to obtain a displacement solution set of blade nodes, constructing a transformation matrix through the displacement solution set of the nodes, and performing order reduction on the finite element model obtained in the step (1) by using the transformation matrix to construct an order reduction model; 3) giving the geometric dimension of the wind turbine generator blade and the actual operation parameters of the wind turbine generator, respectively calculating the transverse displacement of the tail end of the blade by using the reduced order model and the finite element model, comparing the precision of the reduced order model with different reduced order numbers, and selecting a proper reduced order number; a calculation result shows that the calculation cost of transient response analysis is effectively reduced on the premise that the blade transient response calculation precision is guaranteed.

Description

Technical Field

[0001] The present invention relates to the field of on-line monitoring of the operating state of wind turbines, and is a reduced-order modeling method for on-line state monitoring of wind turbine blades. Background Art

[0002] As a large-scale electromechanical equipment that converts wind energy into electrical energy, a wind turbine has the characteristics of complex structure, high manufacturing cost, and harsh service environment. Once a failure occurs, it will cause significant economic losses. Therefore, the on-line monitoring of the operating state of wind turbines has become a research hotspot.

[0003] The blade is a key component for a wind turbine to convert wind energy into kinetic energy, and its health state directly affects the power generation efficiency and operating safety of the entire wind turbine. The working environment of the blade is harsh, and it bears various loads during service, and it is prone to icing, sand erosion, etc., which may lead to faults such as cracking and fracture of the blade. Therefore, the state monitoring of wind turbine blades is of great significance in ensuring the safe and efficient operation of wind turbines and reducing the power generation cost during the life cycle of the unit. When the blade freezes or cracks, etc., the basic properties of the blade will change, such as the blade mass, stiffness, etc., and these changes will cause the transient response of the wind turbine blade to change.

[0004] Traditional numerical calculation methods for the transient response of wind turbine blades (such as the finite element method, the finite volume method, the finite difference method, etc.) have a high number of degrees of freedom in the equilibrium equation after discretization when solving the transient response of the blade, resulting in high calculation costs and unable to achieve the purpose of real-time calculation of the transient response of the blade.

[0005] In order to quickly obtain the transient response of wind turbine blades and master their real-time operating state, the present invention proposes a reduced-order modeling method for on-line state monitoring of wind turbine blades. Summary of the Invention

[0006] The purpose of the present invention is to provide a calculation method for reduced-order modeling of wind turbine blades with high accuracy, fast calculation speed, and good timeliness; the technical solution is as Figure 1 shown.

[0007] The technical solution of the present invention is as follows:

[0008] A reduced-order modeling method for on-line state monitoring of wind turbine blades includes the following steps:

[0009] Step 1: Use the central rigid body-flexible beam model to describe the rotational motion of the hub and blade of the wind turbine, obtain the dynamic equation characterizing the motion of the hub-blade system, and use the finite element method to discretize the dynamic equation of the hub-blade system to construct a finite element model;

[0010] Step 2: Solve the finite element model obtained in Step 1 to obtain the set of displacement solutions of the nodes. Construct a transformation matrix through the set of displacement solutions of the nodes, and use the transformation matrix to reduce the order of the finite element model obtained in Step 1 to construct a reduced-order model;

[0011] Step 3: For an actual wind turbine blade, construct a finite element model and a reduced-order model respectively to calculate the transient response of the blade, compare the accuracy of the reduced-order model at different orders, and determine the final reduced-order order. The final reduced-order order should meet the engineering actual requirements for the reduced-order model at this order. In this patent, it is required that the relative error between the calculated values of the final reduced-order model and the finite element model is less than 1%.

[0012] Further, the process described in Step 1 is as follows: Simplify the blade of the wind turbine into a variable cross-section beam that rotates together with the hub. Since the hub is a curved surface structure and its geometric size is much smaller than the blade length, the hub is simplified as a rigid body, and the deformation of the hub is not considered. Only the influence of the rotation of the hub on the transient response of the blade is considered; Use the central rigid-flexible beam model to describe the motion of the hub-blade system, obtain the dynamic equation of the hub-blade system, and use the finite element method for discretization to construct a finite element model.

[0013] Establish a fixed coordinate system XOY with the center of the wind turbine hub as the origin, and establish a floating coordinate system x1o1y1 with the connection between the wind turbine blade and the hub as the origin, the axial direction of the blade as the x-axis, and the transverse direction as the y-axis. Use the floating coordinate system to describe the displacement of the blade; Due to the influence of the transverse bending of the blade, the axial displacement of a point P0 on the blade is represented by the displacement vector of this point as follows:

[0014]

[0015] In the formula, v1 and v2 are the actual axial displacement and transverse displacement of point P0, respectively, in m; w1, w2, w c are the axial elongation, transverse bending deformation, and axial elongation caused by the transverse bending deformation of point P0, respectively, in m;

[0016] After the blade deforms, point P0 reaches the position of point P. The coordinate vector r2 of point P is expressed in the coordinate system XOY as:

[0017] r2 = Θ(r0 + r1) + R (2)

[0018] In the formula, R is the coordinate vector of the origin o1 of the floating coordinate system x1o1y1 in the coordinate system XOY; Θ is the direction cosine matrix of the floating coordinate system x1o1y1 relative to the coordinate system XOY, r0 is the coordinate vector of point P0 in the floating coordinate system x1o1y1;

[0019] Calculate the kinetic energy T, potential energy H, and external work W of the hub - blade system as follows:

[0020]

[0021] where ρ is the blade density, kg / m 3 ; A(x) is the cross - sectional area of the wind turbine blade at the floating coordinate x, m 2 ; I(x) is the sectional moment of inertia of the blade at the floating coordinate x, m 4 ; J H is the moment of inertia of the hub, kg / m 2 ; τ is the torque applied to the hub, N·m; w′1 is the first - order partial derivative of w1 with respect to x, w″2 is the second - order partial derivative of w2 with respect to x; E is the elastic modulus of the blade material, Pa; θ is the hub rotation angle, rad; is the hub rotational angular velocity, rad / s; is the first - order partial derivative of the vector r2 with respect to time, m / s; L is the total length of the blade, m;

[0022] Based on Equation (3) and using the Hamilton variational principle, as shown in Equation (4); obtain the equations shown in Equations (5), (6), and (7):

[0023]

[0024]

[0025]

[0026]

[0027]

[0028] where:

[0029]

[0030] where, are the first - order partial derivatives of the variables w1, w2, w c with respect to time, m / s; are the second - order partial derivatives of the variables w1, w2 with respect to time, m / s 2 ; is the angular acceleration of the hub rotation, rad / s 2 ; r A is the hub radius, m; w′1 is the first - order partial derivative of w1 with respect to x, w″2, w″″2 are the second - order and fourth - order partial derivatives of w2 with respect to x respectively;

[0031] Discretize equations (5), (6), and (7) using the finite element method. Divide the wind turbine blade along the axial direction into u elements, with the number of nodes being u + 1. The axial elongation w of a point P within blade element i p1 and the lateral displacement w p2 are written in the following form:

[0032]

[0033] where, is the internal coordinate of the beam element; is the shape function matrix of the beam element; q i (t) is the displacement vector of the element nodes at time t. q i (t) is specifically expressed as:

[0034]

[0035] where, respectively represent the axial elongation, lateral displacement, and rotation angle of the I node of the element in the local coordinate system of the element; respectively represent the axial elongation, lateral displacement, and rotation angle of the J node of the element in the local coordinate system of the element; L i is the length of the i-th element of the blade, in m;

[0036] In the floating coordinate system, the displacement vector of point P is expressed as:

[0037]

[0038] where, is the coupling shape function matrix of the i-th element, and its specific form is as follows:

[0039]

[0040] where, L j is the length of the j-th element of the blade, in m;

[0041] The equilibrium equation of the i-th element of the discretized blade is expressed as:

[0042]

[0043] where, is the moment of inertia of the i-th element; and are the non-linear coupling inertia vectors caused by the rotational motion and elastic deformation of the i-th element, and are in a transpose relationship; is the generalized mass matrix of the i-th element; Derived from the gyroscopic effect; is the stiffness matrix of the i-th element of the blade; is the inertial force term generated by the rotational motion of the i-th element; F i is the node generalized load vector of the i-th element; q i 、 are respectively the node displacement vector, node velocity vector, and node acceleration vector of the i-th element;

[0044] The calculation formulas of the relevant variable matrices in the formula are:

[0045]

[0046] In the formula, A i is the cross-sectional area of the i-th element, m 2 ; I i is the cross-sectional moment of inertia of the i-th element, m 4 .

[0047] Assembling all the elements, the overall balance equation of the hub-blade system is obtained:

[0048]

[0049] Furthermore, the process described in step 2 is as follows: Solve the finite element model in step 1 to obtain the set of node displacement solutions, and select m samples from the set to form the matrix X = [q1, q2,..., q m ; where q i = [u1, u2,..., u z T is the displacement solution of the node, z is the number of degrees of freedom of the finite element model; Let {ζ1, ζ2,..., ζ n} be a set of n orthonormal bases of X and write it in matrix form Φ; For the blade finite element model, the node displacement value q t at any time t can be represented by {ζ1, ζ2,..., ζ n}:

[0050]

[0051] In the formula,

[0052] Select the first k orthonormal basis vectors to represent q t There is:

[0053]

[0054] In the formula, Φ k = [ζ1, ζ2,..., ζ​k , k < n; b t(k) = [b t1 , b t1 ,..., b tk T . Let the error function be:

[0055]

[0056] To find the optimal solution of the error function under the constraint condition {ζ1, ζ2,..., ζ n}} being an orthonormal basis, introduce the Lagrange factor u ij (i, j = k + 1, k + 2,…, n), construct the Lagrange function, and find the orthonormal basis corresponding to the optimal solution.

[0057]

[0058] In the formula, δ ij is the Kronecker function.

[0059] Take the partial derivative of both sides of equation (19) with respect to ζ j to obtain:

[0060]

[0061] In the formula, Φ n-k = [ζ k+1 , ζ k+2 ,..., ζ n T , u j = [u k+1,j , u k+2,j ,..., u n,j T . Write the above formula in matrix form to obtain the following formula:

[0062]

[0063] In the formula, U n-k = [u k+1 , u k+2 ,..., u n T . To find the optimal solution, let equation (21) be equal to 0, and multiply both sides on the left by Since the column vectors of Φ n are orthonormal, we can obtain It is easy to know that U n-k is a positive semi - definite matrix, so there exists an orthogonal matrix P that can diagonalize U n-k to obtain Multiply both sides of the above formula on the left by Φ n-k P to obtain the following formula: ​​​​

[0064]

[0065] where Λ is the diagonal matrix of U n-k . By arranging the diagonal elements of Λ in descending order, it can be seen from Equation (22) that the i-th diagonal element of Λ is the eigenvalue λ T of XX i , and Φ n-k P is the matrix composed of the eigenvectors corresponding to λ T of XX i . Considering that the orthogonal matrix preserves the norm under the Frobenius norm , Equation (18) is rewritten as:

[0066]

[0067] Take the left singular value vectors corresponding to the first n singular values of X to form the matrix Φ n . At this time, Φ n is the optimal solution of the error function under the constraint condition {ζ1, ζ2,..., ζ n} being the standard orthogonal basis. The error function has a minimum value and the numerical value is the sum of the last n - k eigenvalues of the matrix XX T . Define the following formula:

[0068]

[0069]

[0070] where λ i is the eigenvalue of the matrix XX T arranged in descending order. If I(k) ≥ d%, it is said that the k-dimensional vector retains d% of the characteristic information of the original sample, where d is a real constant and 0 ≤ d ≤ 100.

[0070] Further, for Step 3, taking the blade of a certain onshore wind turbine as an example, the blade is simplified to a variable cross-section beam, and the hub of the wind turbine is simplified to a central rigid body; a finite element model is constructed according to the geometric dimensions and force conditions of the blade, and the displacement solution set of the nodes is obtained by solving and written as a matrix; the singular value decomposition (SVD) is used for this matrix; the characteristic information of the original samples not retained in the reduced-order models of different orders is calculated according to the singular values, and the reduced-order number is preliminarily determined based on the un-retained characteristic information; a reduced-order model is constructed according to the preliminarily determined reduced-order number, and the error between the reduced-order model and the finite element model is calculated. If the error does not meet the engineering actual requirements, a reduced-order model with a higher order is established; the above process is repeated until the error of the reduced-order model meets the engineering actual requirements. In this patent, when the reduced-order number is set to 10, the reduced-order model can obtain the transient response of the wind turbine blade in real time on the premise of ensuring the calculation accuracy. The calculation results show that the reduced-order modeling method of the wind turbine blade based on the central rigid body-flexible beam model and proper orthogonal decomposition can effectively reduce the calculation cost of obtaining the transient response of the wind turbine blade.

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

[0072] The present invention provides a reduced-order modeling method for online condition monitoring of wind turbine blades, which can quickly and accurately calculate the transient response of wind turbine blades to achieve the purpose of online monitoring of the blade operation state. Brief Description of the Drawings

[0073] Figure 1 It is the technical implementation route map;

[0074] Figure 2 It is a schematic diagram of the central rigid body-flexible beam model;

[0075] Figure 3 It is the curve graph of the angular velocity of the hub rotation;

[0076] Figure 4 It is the curve graph of the transverse displacement at the blade tip obtained by solving the finite element model under no-load conditions;

[0077] Figure 5 It is the logarithmic graph of the singular values;

[0078] Figure 6 It is the logarithmic graph of the characteristic information of the original samples not retained;

[0079] Figure 7 It is the comparison graph of the curve of the transverse displacement at the blade tip obtained by solving the 10th-order reduced-order model and the finite element model under load;

[0080] Figure 8It is a logarithmic value graph of the relative error of the lateral displacement at the blade tip predicted by the 10th-order reduced-order model when the hub speed is 1.27 rad / s;

[0081] Figure 9 It is a logarithmic value graph of the relative error of the lateral displacement at the blade tip predicted by reduced-order models of different orders when the hub speed is 1.27 rad / s;

[0082] Figure 10 It is a logarithmic value graph of the relative error of the lateral displacement at the blade tip predicted by the 10th-order reduced-order model at different hub speeds. Specific implementation manner

[0083] The present invention will be further described below in conjunction with the accompanying drawings of the specification and embodiments.

[0084] A reduced-order modeling method for online condition monitoring of wind turbine blades includes the following steps:

[0085] Step 1: Establish a fixed coordinate system XOY with the center of the wind turbine hub as the origin, and establish a floating coordinate system x1o1y1 with the connection between the wind turbine blade and the hub as the origin, the axial direction of the blade as the x-axis, and the lateral direction as the y-axis. Use the floating coordinate system to describe the displacement of the blade, as Figure 2 shown. Due to the influence of the lateral bending of the blade on the axial displacement of a point P0 on the blade, the displacement vector of this point is expressed as:

[0086]

[0087] In the formula, v1 and v2 are respectively the actual axial displacement and lateral displacement of point P0, in m; w1, w2, w c are respectively the axial elongation, lateral bending deformation, and axial elongation caused by the lateral bending deformation of point P0, in m. Since the axial dimension of the blade is much larger than the radial dimension, it can be simply considered that v2(x,t) = w2(x,t), that is, the lateral displacement is equal to the lateral bending deformation. In the present invention, the lateral displacement and the lateral bending deformation are collectively referred to as the lateral displacement.

[0088] After the blade deforms, point P0 reaches the position of point P. The coordinate vector r2 of point P is expressed in the coordinate system XOY as:

[0089] r2 = Θ(r0 + r1) + R (2)

[0090] In the formula, R is the coordinate vector of the origin o1 of the floating coordinate system x1o1y1 in the coordinate system XOY; Θ is the direction cosine matrix of the floating coordinate system x1o1y1 relative to the coordinate system XOY, r0 is the coordinate vector of point P0 in the floating coordinate system x1o1y1.

[0091] Calculate the kinetic energy T, potential energy H, and external work W of the hub - blade system as follows:

[0092]

[0093] where ρ is the blade density, kg / m 3 ; A(x) is the cross - sectional area of the wind turbine blade at the floating coordinate x, m 2 ; I(x) is the sectional moment of inertia of the blade at the floating coordinate x, m 4 ; J H is the moment of inertia of the hub, kg / m 2 ; τ is the torque applied to the hub, N·m; w′1 is the first - order partial derivative of w1 with respect to x, w″2 is the second - order partial derivative of w2 with respect to x; E is the elastic modulus of the blade material, Pa; θ is the hub rotation angle, rad; is the angular velocity of hub rotation, rad / s; is the first - order partial derivative of the vector r2 with respect to time, m / s; L is the total length of the blade, m.

[0094] Based on Equation (3) and using Hamilton's variational principle (as shown in Equation (4)), the equations shown in Equations (5), (6), and (7) can be obtained:

[0095]

[0096]

[0097]

[0098]

[0099] where:

[0100]

[0101] where, are the first - order partial derivatives of the variables w1, w2, w c with respect to time, m / s; are the second - order partial derivatives of the variables w1, w2 with respect to time, m / s 2 ; is the angular acceleration of hub rotation, rad / s 2 ; r A is the hub radius, m; w′1 is the first - order partial derivative of w1 with respect to x, w″2, w″″2 are the second - order and fourth - order partial derivatives of w2 with respect to x respectively.

[0102] The equations (5), (6), and (7) are discretized using the finite element method. The wind turbine blade is axially divided into u elements, and the number of nodes is u + 1. The axial elongation w p1 and the lateral displacement w p2 of a point P within the blade element i are written in the following form:

[0103]

[0104] where is the internal coordinate of the beam element; is the shape function matrix of the beam element; q i (t) is the displacement vector of the element nodes at time t. q i (t) is specifically expressed as:

[0105]

[0106] where respectively represent the axial elongation, lateral displacement, and rotation angle of the I node of the element in the local coordinate system of the element; respectively represent the axial elongation, lateral displacement, and rotation angle of the J node of the element in the local coordinate system of the element; L i is the length of the i-th element of the blade, in m.

[0107] Therefore, in the floating coordinate system, the displacement vector of point P can be expressed as:

[0108]

[0109] where is the coupling shape function matrix of the i-th element, and its specific form is as follows:

[0110]

[0111] where L j is the length of the j-th element of the blade, in m. The equilibrium equation of the i-th element of the discretized blade can be expressed as:

[0112]

[0113] where is the moment of inertia of the i-th element; and are the non-linear coupling inertia vectors caused by the rotational motion and elastic deformation of the i-th element, and are in a transpose relationship; is the generalized mass matrix of the i-th element; results from the gyroscopic effect; is the stiffness matrix of the i-th element of the blade; is the inertial force term generated by the rotational motion of the i-th element; F i is the generalized load vector of the nodes of the i-th element; q i 、 are the node displacement vector, node velocity vector, and node acceleration vector of the i-th element, respectively. The calculation formulas of the relevant variable matrices in the formula are:

[0114]

[0115] In the formula, A i is the cross-sectional area of the i-th element, m 2 ; I i is the cross-sectional moment of inertia of the i-th element, m 4 .

[0116] Assemble all the elements to obtain the overall equilibrium equation of the hub-blade system:

[0117]

[0118] For the blades of a wind turbine, the motion law of the hub is known. Therefore, a non-inertial system is established with the hub as the reference system, and the following formula is obtained:

[0119]

[0120] Use the Newmark method to construct an integral solution format for Equation (25) to obtain the following formula:

[0121] Z i+1 q i+1 = Q i+1 + F i+1 + R i+1,1 + R i+1,2 + R i+1,3 (26)

[0122] In the formula:

[0123]

[0124] In the formula, where Δt is the integration time step, and δ and α are two control parameters in the Newmark method.

[0125] Use Equation (14) to update matrices M, G, and K and solve Equation (26) to obtain the set of displacement solutions of the nodes.

[0126] Step 2: Select m node displacement vectors from the set of displacement solutions of the nodes to form matrix X m = [q1, q2,..., q m, perform singular value decomposition on matrix X, and select the left singular value vectors corresponding to the first k left singular values to form the transformation matrix Φ k =[ψ1,ψ2,...,ψ k .

[0127] Use the transformation matrix Φ k to reduce the order of the finite element model to obtain a reduced-order model, and use the direct integration method to construct an integration solution format to obtain the following formula:

[0128]

[0129] In the formula:

[0130]

[0131] Solve the above formula to obtain the coefficient vector b i+1 , and the displacement solution of the node can be reconstructed according to Equation (17).

[0132] Step 3: The operating parameters of an onshore wind turbine with a rated power of 5 MW are shown in Table 1:

[0133] Table 1: Operating parameters of 5 MW onshore wind turbine

[0134] Cut-in wind speed 3 m / s Rated wind speed 11.4 m / s Cut-out wind speed 25 m / s Rated blade rotational speed 12 r / min, 1.27 rad / s Hub height above ground 87.6m Distance from yaw bearing to hub center 1.75m

[0135] The structural information of the blades of this wind turbine is shown in Table 2:

[0136] Table 2: Blade structural properties

[0137]

[0138] The blade material is glass fiber reinforced composite material, and the elastic modulus of this material is 73.9 GPa and the density is 1510 kg / m 3 ; EI is the sectional flexural stiffness in the flapwise direction of the blade, N·m 2 ; EA is the tensile stiffness, N.

[0139] According to the actual operating conditions of the wind turbine, the function of the hub rotational angular velocity is given as follows, and the function is plotted in Figure 3 :

[0140]

[0141] In the formula, T = 15 s, which is the duration of the blade rotation acceleration stage; ω = 1.27 rad / s, which is the rotational speed of the blade during stable rotation; T1 = 20 s, which is the total duration of numerical calculation.

[0142] According to the data in Table 2, a finite element model of the blade was constructed, with a total of 169 elements and 170 nodes; fixed constraints were applied to the connection between the wind turbine blade and the hub, and the degree of freedom of the final blade finite element model was 507. When calculating, let Δt = 0.001 s, α = 0.25, and δ = 0.5. The transverse displacement at the blade tip calculated for the blade finite element model without applying a load is as Figure 4 shown.

[0143] The gravity load and aerodynamic load were applied to the blade finite element model for solution, and the displacement solution set of the nodes was obtained. The set was written in matrix form and subjected to singular value decomposition. The first twenty singular values were taken and numbered 1 - 20 in descending order; the base-10 logarithms of the first twenty singular values were taken and plotted in Figure 5 . It can be seen from the figure that the magnitude of the singular values does not stabilize around a certain value as the serial number increases, but the value of the 20th singular value is much smaller than the value of the 1st singular value. The transformation matrix containing different numbers of left singular value vectors was calculated according to the following formula to obtain the characteristic information U(k) of the original sample not retained, and the result is as Figure 6 shown:

[0144] U(k)=log 10 ((1 - I(k))×100%) (31)

[0145] where the definition of I(k) is given by Equation (24).

[0146] From Figure 6 , it can be seen that when the order reduction order is greater than or equal to 10, the change in the characteristic information value of the original sample not retained by the transformation matrix is small. Through calculation, it can be known that the information volume of the transformation matrix with an order reduction order of 10 not retaining the characteristics of the original sample is less than 0.01%. Therefore, the order reduction order is initially determined to be 10.

[0147] A 10th-order reduced-order model was constructed using the transformation matrix, and the error between the 10th-order reduced-order model and the finite element model was calculated. Whether a higher-order reduced-order model is required was judged according to the calculation results.

[0148] The transient responses of the blade under the condition of the stable rotational speed ω = 1.27 rad / s of the hub were calculated using the finite element model and the 10th-order reduced-order model respectively. The calculation time of the 10th-order reduced-order model was 0.1 s, and the calculation time of the finite element reduced-order model was 71 s. The transverse displacements at the blade tip in the calculation results of the finite element model and the 10th-order reduced-order model were plotted to obtain Figure 7 the curve shown. It can be known that the calculation results of the 10th-order reduced-order model are in good agreement with the calculation results of the finite element model.

[0149] Further quantitative calculations are made on the calculation errors of the 10th-order reduced-order model, and the relative errors of the reduced-order model at different times are calculated. Since the numerical values of most of the relative errors are relatively close, in order to clearly represent the relative errors of the reduced-order model, the logarithmic values of the relative errors are used, as defined in Equation (32), and the calculation results are as follows Figure 8 shown.

[0150]

[0151] where q t(10) is the calculated value of the 10th-order reduced-order model for the lateral displacement of the blade tip node at the t-th moment, and q t is the calculated value of the finite element model for the lateral displacement of the blade tip node at the t-th moment.

[0152] Under the condition that other conditions remain unchanged, reduced-order models with reduced-order numbers of 5, 10, 15, and 20 are used to predict the lateral displacement at the tip of the blade respectively. According to Equation (32), the logarithmic values of the relative errors of the reduced-order models with different orders for the lateral displacement at the tip of the blade are calculated, and the calculation results are plotted in Figure 9 .

[0153] Under the condition that other conditions remain unchanged, the 10th-order reduced-order model is used to predict when the stable rotational angular velocity ω of the hub is 0.1 rad / s, 0.4 rad / s, 0.7 rad / s, and 1.27 rad / s respectively. Similarly, according to Equation (32), the logarithmic values of the relative errors for the lateral displacement at the tip of the blade are calculated, and the calculation results are plotted in Figure 10 .

[0154] From Figure 8 it can be seen that the logarithmic values of the relative errors of the 10th-order reduced-order model for calculating the lateral displacement at the tip of the blade are concentrated around -4, and the maximum logarithmic value of the relative error does not exceed -2. Through calculation, it can be known that the maximum relative error is 0.74×10 -2 , and the error meets the engineering actual requirements; from Figure 9 it can be seen that the distribution regions of the logarithmic values of the relative errors of the reduced-order models with reduced-order numbers of 10, 15, and 20 mostly overlap, while when the reduced-order number is 5, the distribution region of the logarithmic values of the relative errors of the reduced-order model is above the distribution regions of the logarithmic values of the relative errors of the reduced-order models with orders of 10, 15, and 20, and does not overlap with the distribution regions of the logarithmic values of the relative errors of the reduced-order models with orders of 10, 15, and 20. It can be seen that it is reasonable to select the reduced-order number of 10; from Figure 10 it can be seen that although the angular velocities during the stable rotation of the hub are different, the logarithmic values of the relative errors predicted by the 10th-order reduced-order model are distributed between -8 and -2, and are concentrated around -4. It can be seen that the 10th-order reduced-order model can effectively calculate the transient response of the blade.

[0155] In summary, the 10th-order reduced-order model constructed by the blade reduced-order modeling method of wind turbines based on the central rigid body-flexible beam model and proper orthogonal decomposition can greatly reduce the computational cost of blade transient response analysis while meeting the computational accuracy requirements.

Claims

1. A reduced-order modeling method for on-line condition monitoring of wind turbine blades, characterized in that, It includes the following steps: Step 1: Use the central rigid body-flexible beam model to describe the rotational motion of the hub and blades of the wind turbine, obtain the dynamic equation characterizing the motion of the hub-blade system, and use the finite element method to discretize the dynamic equation of the hub-blade system to construct a finite element model; Step 2: Solve the finite element model obtained in Step 1 to obtain the set of displacement solutions of the nodes, construct a transformation matrix through the set of displacement solutions of the nodes, and use the transformation matrix to reduce the order of the finite element model obtained in Step 1 to construct a reduced-order model; Step 3: For the blades of the wind turbine, respectively construct a finite element model and a reduced-order model to calculate the transient response of the blades, compare the accuracy of the reduced-order models at different orders, and determine the final reduced-order order. The final reduced-order order should meet the engineering actual requirements for the reduced-order model at this order.

2. A reduced-order modeling method for on-line condition monitoring of wind turbine blades according to claim 1, characterized in that The process in Step 1 is as follows: Simplify the blade of the wind turbine into a variable cross-section beam rotating together with the hub, simplify the hub into a central rigid body, and only consider the influence of the hub rotation on the transient response of the blade; Establish a fixed coordinate system XOY with the center of the wind turbine hub as the origin, and establish a floating coordinate system x1o1y1 with the connection between the wind turbine blade and the hub as the origin, the axial direction of the blade as the x-axis, and the transverse direction as the y-axis, and use the floating coordinate system to describe the displacement of the blade; The axial displacement of a point P0 on the blade is affected by the transverse bending of the blade, and the displacement vector of this point is expressed as: Wherein, v1 and v2 are respectively the actual axial displacement and lateral displacement of point P0, in m; w1, w2, and w c are respectively the axial elongation, lateral bending deformation, and axial elongation caused by lateral bending deformation of point P0, in m; After the blade deforms, point P0 reaches the position where point P is located. The coordinate vector r2 of point P is expressed in the coordinate system XOY as: r2 = Θ(r0 + r1) + R (2) Wherein, R is the coordinate vector of the origin o1 of the floating coordinate system x1o1y1 in the coordinate system XOY; Θ is the direction cosine matrix of the floating coordinate system x1o1y1 relative to the coordinate system XOY, r0 is the coordinate vector of the point P0 in the floating coordinate system x1o1y1; Calculate the kinetic energy T, potential energy H, and external work W of the hub-blade system as follows: where ρ is the blade density, kg / m 3 ; A(x) is the cross-sectional area of the wind turbine blade at x in the floating coordinate system, m 2 ; I(x) is the sectional moment of inertia of the blade at x in the floating coordinate system, m 4 ; J H is the moment of inertia of the hub, kg / m 2 ; τ is the torque applied to the hub, N·m; w1′ is the first-order partial derivative of w1 with respect to x, w″2 is the second-order partial derivative of w2 with respect to x; E is the elastic modulus of the blade material, Pa; θ is the hub rotation angle, rad; is the hub rotational angular velocity, rad / s; is the first-order partial derivative of the vector r2 with respect to time, m / s; L is the total length of the blade, m; Based on Equation (3) and using Hamilton's variational principle as shown in Equation (4); obtain the equations shown in Equations (5), (6), and (7): Where: wherein, are the first-order partial derivatives of variables w1, w2, w c with respect to time, m / s; are the second-order partial derivatives of variables w1, w2 with respect to time, m / s 2 ; is the angular acceleration of the hub rotation, rad / s 2 ; r A is the hub radius, m; w1′ is the first-order partial derivative of w1 with respect to x, and w″2, w″″2 are the second-order and fourth-order partial derivatives of w2 with respect to x respectively; Discretize equations (5), (6), and (7) using the finite element method. Divide the wind turbine blade axially into u elements with u + 1 nodes. The axial elongation w p1 and the lateral displacement w p2 of a point P within blade element i are written in the following form: In the formula, is the internal coordinate of the beam element; is the shape function matrix of the beam element; q i (t) is the displacement vector of the element node at time t. q i (t) is specifically expressed as: In the formula, respectively represent the axial elongation, lateral displacement and rotation angle of the I node of the element in the local coordinate system of the element; respectively represent the axial elongation, lateral displacement and rotation angle of the J node of the element in the local coordinate system of the element; L i is the length of the i-th element of the blade, in m; In the floating coordinate system, the displacement vector of point P is expressed as: In the formula, is the coupling shape function matrix of the i-th unit, and its specific form is as follows: where L j is the length of the j-th element of the blade, in m; The equilibrium equation of the i-th element of the blade after discretization is expressed as: Wherein, is the moment of inertia of the i-th unit; and is the non-linear coupling inertia vector caused by the rotational motion and elastic deformation of the i-th unit, and are in a transpose relationship with each other; is the generalized mass matrix of the i-th unit; results from the gyroscopic effect; is the stiffness matrix of the i-th unit of the blade; is the inertial force term generated by the rotational motion of the i-th unit; F i is the nodal generalized load vector of the i-th unit; q i 、 are respectively the nodal displacement vector, nodal velocity vector, and nodal acceleration vector of the i-th unit; The calculation formulas of the relevant variable matrices in the formula are: Where, A i is the cross-sectional area of the i-th unit, m 2 ; I i is the sectional moment of inertia of the i-th unit, m 4 . Assemble all the elements to obtain the overall equilibrium equation of the hub-blade system:

3. A reduced-order modeling method for on-line condition monitoring of wind turbine blades according to claim 1, characterized in that The process described in Step 2 is as follows: Solve the finite element model in Step 1 to obtain the set of nodal displacement solutions, and select m samples from the set to form the matrix X = [q1, q2,..., q m ; where q i = [u1, u2,..., u z T is the displacement solution of the node, and z is the number of degrees of freedom of the finite element model; Let {ζ1, ζ2,..., ζ n} be a set of n orthonormal bases of X and write it in matrix form Φ; For the blade finite element model, the nodal displacement value q t at any time t can be represented by {ζ1, ζ2,..., ζ n}:​ In the formula, Select the first k orthonormal basis vectors to represent q t There is: where Φ k = [ζ1, ζ2,..., ζ k , k < n; b t(k) = [b t1 , b t1 ,..., b tk T ;​ Let the error function be: To find the optimal solution of the error function under the constraint conditions {ζ1, ζ2,..., ζ n}, Lagrange factor u ij (i, j = k + 1, k + 2,..., n) is introduced to construct the Lagrange function and find the orthonormal basis corresponding to the optimal solution; where δ ij is the Kronecker function; Take the partial derivative of both sides of Equation (19) with respect to ζ j to obtain: where, Φ n-k = [ζ k+1 , ζ k+2 ,..., ζ n T , u j = [u k+1,j , u k+2 , j,..., u n,j T ;​​ Write the above formula in matrix form to obtain the following formula: where U n-k = [u k+1 , u k+2 ,..., u n T ;​ To find the optimal solution, let Equation (21) equal to 0 and multiply both sides on the left by Since Φ n is a column vector of orthonormal vectors, we get It is easy to know that U n-k is a positive semi - definite matrix, so there exists an orthogonal matrix P that can diagonalize U n-k to get Multiply both sides of the above equation on the left by Φ n-k P to get the following equation: where Λ is the diagonal matrix of U n-k ; Arrange the diagonal elements of Λ in descending order. From Equation (22), the i-th diagonal element of Λ is XX T eigenvalue λ of i , while Φ n- k P is XX T corresponding to λ i matrix composed of eigenvectors; Considering that the orthogonal matrix preserves the norm under the Frobenius norm Therefore, rewrite Equation (18) as follows: Form a matrix Φ by taking the left singular value vectors corresponding to the first n singular values of X n , where Φ n is the optimal solution of the error function under the constraint that {ζ1, ζ2,..., ζ n} is an orthonormal basis. The error function has a minimum value, and the numerical value is the sum of the last n - k eigenvalues of the matrix XX T ; Define the following formula: where λ i is the eigenvalue of matrix XX T sorted in descending order; If I(k) ≥ d%, then it is said that the k-dimensional basis vector retains d% of the characteristic information of the original sample, where d is a real constant, 0 ≤ d ≤ 100.

4. A reduced-order modeling method for on-line condition monitoring of wind turbine blades according to claim 1, characterized in that In Step 3, simplify the blade into a variable cross-section beam and simplify the wind turbine hub into a central rigid body; construct a finite element model according to the geometric dimensions and force conditions of the blade, solve to obtain the set of displacement solutions of the nodes, and write the set as a matrix; Perform singular value decomposition (SVD) on the matrix; calculate the characteristic information of the original samples not retained in the reduced-order models of different orders according to the singular values, and preliminarily determine the reduced order based on the un-retained characteristic information; Construct a reduced-order model according to the preliminarily determined reduced order, calculate the error between the reduced-order model and the finite element model. If the error does not meet the engineering actual requirements, establish a reduced-order model with a higher order; repeat the above process until the error of the reduced-order model meets the engineering actual requirements.

Images