DIGITAL TWINNING METHOD FOR MONITORING OPERATION STATE OF TOWER OF WIND TURBINE GENERATOR SYSTEM ONLINE

20250165673 ยท 2025-05-22

Assignee

Inventors

Cpc classification

International classification

Abstract

Disclosed is a digital twinning method for monitoring an operation state of a tower of a wind turbine generator system online. The method includes: 1) constructing a simplified model of the tower of the wind turbine generator system, and discretizing the simplified model according to a finite element method to obtain a finite element model; 2) reducing an order of the finite element model of the tower according to proper orthogonal decomposition, analyzing precision of a reduced-order model under different orders, and selecting the reduced-order model having a smallest reduced order as a final reduced-order model on the premise that the precision satisfies actual engineering requirements; and 3) programming upper computer software in a computer, deploying the reduced-order model to the upper computer software, further building a physical entity of the tower, and monitoring a stress and a strain of the physical entity online through the reduced-order model.

Claims

1. A digital twinning method for monitoring an operation state of a tower of a wind turbine generator system online, comprising: step 1, obtaining a simplified physical model of the tower of the wind turbine generator system according to geometrical characteristics of the tower of the wind turbine generator system, and analyzing an actual stress condition of the tower to construct a finite element model of the tower; step 2, carrying out order reduction analysis on the finite element model constructed in step 1 according to proper orthogonal decomposition (POD), comparing precision of a reduced-order model under different orders, and determining a final reduced order, wherein the final reduced order satisfies engineering requirements of the reduced model under the order; and step 3, designing and constructing upper computer software in a computer, and deploying the reduced-order model of the tower to the upper computer software; building a physical entity of the tower, installing a strain sensor and an inclination sensor, reading sensor data and uploading the data to an upper computer through a single chip microcomputer, so as to convert the data into a displacement boundary condition, and then inputting the displacement boundary condition into the reduced-order model for calculation; and determining accuracy of the digital twinning method for monitoring the operation state of the tower by calculating an error between a measured strain value of the strain sensor and a calculated strain value of the reduced-order model in the upper computer.

2. The digital twinning method for monitoring an operation state of a tower of a wind turbine generator system online according to claim 1, wherein the step 1 further comprises: simplifying a flange connection structure between the tower of the wind turbine generator system, and simplifying the tower into a beam structure, so as to construct the simplified physical model of the tower; a point of intersection of an axis of the tower before bending deformation and a lower end surface of the tower is an origin O of a coordinate system XYZ; and normal vectors of an end surface at a top end of the tower before bending deformation and after bending deformation are n.sub.0 and n.sub.1 respectively, and an included angle between the two vectors is an inclination angle of the end surface at the top end of the tower, and an expression of deflection w and the inclination angle of the end surface at the top end of the tower is as follows: { w = F 6 EI L 3 - FL 2 EI L 2 = F 2 EI L 2 - FL EI L .Math. w = 2 3 L ( 1 ) wherein F is a concentrated load applied at the top end of the tower, and has a unit of N; E is an elastic modulus of a tower material, and has a unit of Pa; I is cross sectional moment of inertia of the tower, and has a unit of m.sup.4; and L is a length of the tower, and has a unit of m; the deflection w and an azimuth need to be determined to describe the bending deformation of the tower in a three-dimensional space, wherein the azimuth is configured to express a bending orientation of the tower; and therefore, the displacement boundary condition of the finite element model of the tower is determined by the two parameters of w and , and according to a principle of linear superposition and decomposition of force, the displacement boundary condition of the tower is determined to obtain: q b = w x q bx + w y q by = ( w cos ) q bx + ( w sin ) q by ( 2 ) wherein q.sub.bx=[1, 0, 1, 0, . . . , 1, 0].sup.T, q.sub.by=[0, 1, 0, 1, . . . , 0, 1].sup.T, and a vector length is determined by a specific displacement boundary condition; q.sub.b is a joint displacement vector determined by the displacement boundary condition of the tower, i.e. a joint displacement vector corresponding to a joint having a displacement constraint of the tower; w.sub.x is a component of total deflection in a direction X; and w.sub.y is a component of the total deflection in a direction Y; and the deflection w at the top end of the tower is converted into the displacement boundary condition, the finite element model of the tower is constructed, and an equilibrium equation of the finite element model of the tower under the displacement boundary condition is shown in formula (3): [ K a a K a b K b a K b b ] [ q a q b ] = [ g a g b ] + [ 0 R r ] ( 3 ) wherein K.sub.aa is a stiffness matrix corresponding to a joint node having no displacement constraint of the tower; K.sub.ab and K.sub.ba are coupling matrices between the joint having the displacement constraint and the node having no displacement constraint of the tower; K.sub.bb is a stiffness matrix corresponding to the joint having the displacement constraint of the tower; q.sub.a is a joint displacement vector of the joint having no displacement constraint of the tower; g.sub.a and g.sub.b are joint gravity load vectors corresponding to the joint having no displacement constraint and the joint having the displacement constraint of the tower respectively; and R.sub.r is a joint support reaction force vector corresponding to the joint having the displacement constraint of the tower.

3. The digital twinning method for monitoring an operation state of a tower of a wind turbine generator system online according to claim 1, wherein the step 2 further comprises: constructing the reduced-order model of a displacement field of the tower by combining the POD with a finite element method as follows: solving the finite element model of step 1 to obtain a displacement solution set of a joint, and selecting m samples from the set to form a matrix X=[q.sub.1, q.sub.2, . . . , q.sub.m], wherein q.sub.i=[u.sub.1, u.sub.2, . . . , u.sub.z].sup.T is a displacement solution of the joint, and z is the number of degrees of freedom of the finite element model; {.sub.1, .sub.2, . . . , .sub.n} is set as a group of n orthonormal bases of X, and is written as a matrix form denoted as ; and for the finite element model of the tower, a joint displacement value q.sub.t at any moment is represented by {.sub.1, .sub.2, . . . , .sub.n}: q t = .Math. i - 1 n b ti i ( 4 ) wherein b.sub.ti=.sub.i.sup.Tq.sub.t; selecting first k orthonormal basis vectors from {.sub.1, .sub.2, . . . , .sub.n} to represent q.sub.t: q t ( k ) = .Math. i = 1 k b ti i = k b t ( k ) ( 5 ) wherein .sub.k=[.sub.1, .sub.2, . . . , .sub.k]; k<n; and b.sub.t(k)=[b.sub.t1, b.sub.t1, . . . , b.sub.tk].sup.T; establishing an error function as follows: 2 = .Math. t = 1 n .Math. q t - q t ( k ) .Math. 2 2 = .Math. j = k + 1 n j T XX T j ( 6 ) in order to obtain an optimal solution of the error function under a constraint condition {.sub.1, .sub.2, . . . , .sub.n} of the orthogonal bases, introducing Lagrange factors u.sub.ij (i, j=k+1, k+2, . . . , n) to construct a Lagrange function, so as to find the orthonormal bases corresponding to the optimal solution: L = .Math. j = k + 1 n j T XX T j - .Math. i = k + 1 n .Math. j = k + 1 n [ u ij ( i T j - ij ) ] ( 7 ) wherein .sub.ij is a Kronecker function; solving a partial derivative of .sub.j from two ends of formula (7) to obtain: L j = 2 ( XX T j - .Math. i = k + 1 n u ij j ) = 2 XX T j - 2 n - k u j ( 8 ) .sub.nk=[.sub.k+1, .sub.k+2, . . . , .sub.n].sup.T, and u.sub.j=[u.sub.k+1,j, u.sub.k+2,j, . . . , u.sub.n,j].sup.T; writing the formula into a matrix form to obtain the following formula: L n - k = 2 ( XX T n - k - n - k U n - k ) ( 9 ) wherein U.sub.nk=[u.sub.k+1, u.sub.k+2, . . . , u.sub.n].sup.T; in order to solve the optimal solution, making equation (9) equal to 0, and multiplying the equation left by .sub.nk.sup.T, wherein internal vectors of a matrix .sub.n are orthonormal, and therefore U.sub.nk=.sub.nk.sup.TXX.sup.T.sub.nk=(X.sup.T.sub.nk).sup.TX.sup.T.sub.nk is obtained, and it is easy to know that U.sub.nk is a positive semi-definite matrix, and therefore an orthogonal matrix P exists to diagonalize U.sub.nk, so as to obtain P.sup.T.sub.nk.sup.TXX.sup.T.sub.nkP=P.sup.TU.sub.nkP=; and multiply two ends of the formula left by .sub.nkP to obtain the following formula: n - k PP T n - k T XX T n - k P = n - k PP T U n - k P = n - k P ( 10 ) XX T n - k P = n - k U n - k P = n - k P wherein is a diagonal matrix of U.sub.nk; arranging diagonal elements of in descending order, wherein it is seen from formula (10) that an i-th diagonal element of is an eigenvalue .sub.i of XX.sup.T, and an i-th column of a matrix .sub.nkP is an eigenvector of XX.sup.T corresponding to .sub.i; and considering that the orthogonal matrix is norm-preserving under a Frobenius norm .sub.F.sup.2, formula (6) is rewritten as: 2 = .Math. X T n - k .Math. F 2 = .Math. X T n - k P .Math. F 2 = trace { ( X T n - k P ) T X T n - k P } = trace { } ( 11 ) forming a transformation matrix .sub.n from left singular value vectors corresponding to first n singular values of X, wherein .sub.n is an optimal solution of the error function under the condition that the constraint condition {.sub.1, .sub.2, . . . , .sub.n} is the orthonormal bases, and a minimum value of the error function is the sum of last n-k eigenvalues of the matrix XX.sup.T; defining the following formula: , I ( k ) = .Math. i = 1 k i / .Math. i = 1 n i ( 12 ) wherein .sub.i is an eigenvalue of XX.sup.T arranged in descending order; and retaining, by a k-dimensional base vector, d % of characteristic information of an original sample under the condition of I(k)d %; and constructing the reduced-order model under different orders, and calculating precision of the reduced-order model under different reduced orders separately to determine the appropriate reduced order and obtain the reduced-order model, wherein the finally determined reduced order satisfies actual engineering requirements of the reduced-order model under the order.

4. The digital twinning method for monitoring an operation state of a tower of a wind turbine generator system online according to claim 1, wherein in the step 3, the upper computer reads the sensor data uploaded by the single chip microcomputer, converts the data into the displacement boundary condition, and inputs the displacement boundary condition into the reduced-order model to achieve rapid calculation, so as to obtain calculated strain and stress values of the current physical entity of the tower and a displacement solution of a joint; after calculation of the reduced-order model is completed, a calculation result is processed in the upper computer, that is, deformation of the physical entity of the tower is rapidly quantified on the basis of the displacement solution of the joint and displayed in the upper computer; an error between a strain calculation result of the reduced-order model and an actual strain of the physical entity is obtained by means of a strain value calculated by the reduced-order model and a strain value measured by a sensor; and stress nephogram is visualized by combining the displacement solution of the joint with the calculated stress value to intuitively and clearly express stress distribution of the physical entity of the tower.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

[0033] FIG. 1 is a road map of a technical implementation of the present disclosure.

[0034] FIG. 2 is a schematic diagram of deformation of a tower of a wind turbine generator system.

[0035] FIG. 3 is a graph of logarithmic distribution of singular values.

[0036] FIG. 4 is a graph of characteristic information of each order of transformation matrix without retaining original samples.

[0037] FIG. 5 is a diagram of relative error logarithmic values of a 5-order reduced-order model at each degree of freedom.

[0038] FIG. 6 is a schematic diagram of a physical entity of a simplified tower, a diagram of the physical entity, and a diagram of a digital twinning model.

[0039] FIG. 7 is an interface diagram of upper computer software.

[0040] FIG. 8 is a diagram of calculation of online monitoring of a maximum von-mises stress of a simplified tower.

[0041] FIG. 9 is a diagram of calculation of online monitoring of a strain of a simplified tower.

DETAILED DESCRIPTIONS OF THE EMBODIMENTS

[0042] The present disclosure will be further described below in combination with the accompanying drawings of the description and the examples.

[0043] A digital twinning method for monitoring an operation state of a tower of a wind turbine generator system online specifically includes:

[0044] Step 1: a simplified physical model of the tower of the wind turbine generator system is established, the tower is approximated to a slender beam having a flange structure at one end, a beam cross section is in a shape of a circular ring, and a material is polyvinyl chloride (PVC); and the simplified model is divided into meshes through a 3D8 joint entity structure unit to obtain a mesh model including 7080 joints and 3780 units.

[0045] In order to calculate a joint displacement vector q.sub.a of a joint having no displacement constraint of the tower, formula (3) is deformed to obtain the following formula:

[00013] { K aa q a = g a - K ab q b K bb q b = g b + R r - K ba q a .Math. K aa q a = g a - K ab q b ( 13 )

[0046] Deflection w at a top end of the tower and a bending azimuth of the tower are sampled according to a Latin hypercube sampling method to determine a sampling range of the deflection as [0, 0.15], a sampling range of an azimuth as [0, 2), and the number of sampling points as 20. Results are shown in Table 1:

TABLE-US-00001 TABLE 1 Sampling results obtained according to Latin hypercube sampling method for deflection at top end of tower and bending azimuth of tower Sampling point sequence number Deflection w/m Azimuth /rad 1 0.123 0.317 2 0.072 0.888 3 0.054 1.489 4 0.081 2.074 5 0.024 4.060 6 0.036 3.554 7 0.139 2.701 8 0.145 2.832 9 0.091 5.842 10 0.021 5.145 11 0.045 4.349 12 0.061 1.765 13 0.131 2.462 14 0.004 0.075 15 0.011 6.009 16 0.042 5.544 17 0.108 4.497 18 0.116 3.363 19 0.098 1.089 20 0.084 4.765

[0047] Data in Table 1 is substituted into formula (2) to obtain a displacement boundary condition, and then substituted into formula (13) for calculation to obtain a displacement solution set of joints, and the displacement solution set of a joint is used in step 2 to construct a reduced-order model. Physical quantities finally monitored online in the present disclosure are a stress and a strain, and the stress and the strain at any unit joint may be calculated according to a displacement solution of the joint with reference to formula (14):

[00014] ( x 0 , y 0 , z 0 ) = [ xx yy zz xy yz zx ] = D ( x 0 , y 0 , z 0 ) = D [ xx yy zz xy yz zx ] = D [ ] N ( x 0 , y 0 , z 0 ) q ( 14 ) [0048] (x.sub.0, y.sub.0, z.sub.0) is the stress at the unit joint at a coordinate (x.sub.0, y.sub.0, z.sub.0); D is a material elastic coefficient matrix; (x.sub.0, y.sub.0, z.sub.0) is the strain at the unit joint at the coordinate (x.sub.0, y.sub.0, z.sub.0); [] is an operator matrix of a geometric equation; N(x.sub.0, y.sub.0, z.sub.0) is a shape function matrix at the coordinate (x.sub.0, y.sub.0, z.sub.0); and .sub.xx, .sub.yy, .sub.zz, .sub.xy, .sub.yz and .sub.zx are stress components at the unit joint separately, and .sub.xx, .sub.yy, .sub.zz, .sub.xy, .sub.yz, and .sub.zx are strain components at the unit joint separately.

[0049] Step 2: according to Table 1, the displacement boundary condition is determined, a finite element model is solved, the calculated displacement solution of the joint is sorted corresponding to sampling point sequence numbers to form a matrix V=[v.sub.1, v.sub.1, . . . , v.sub.20], and the matrix V is decomposed through singular values to obtain a singular value matrix and a left singular value vector matrix W.

[0050] First 20 singular values on a diagonal of the singular value matrix are drawn in FIG. 3. Due to a large difference in numerical values, logarithmic values are drawn after a base 10 is taken as a logarithm, sequence numbers are given according to positions in which the singular values are located in the matrix , and the sequence number of the first singular value on the diagonal is given as 1, and so on. It may be seen through analysis of results obtained in FIG. 3 that the singular values having sequence numbers 6 to 20 are close in size. Left singular value vectors corresponding to the first k singular values are selected to construct a transformation matrix, the reduced-order model is constructed through the transformation matrix, and characteristic information U(k) of the transformation matrix without retaining original samples is calculated according to the following formula:

[00015] U ( k ) = log 10 ( ( 1 - I ( k ) ) 100 % ) ( 15 )

[0051] When k is less than 10, the characteristic information of the transformation matrix without retaining the original samples is shown in FIG. 4. It may be seen from FIG. 4 that when a value of k is greater than or equal to 5, the characteristic information of the transformation matrix without retaining the original samples is changed slightly. It may be seen through calculation that when the reduced order is 5, the transformation matrix does not retain 0.0001% of the characteristic information of the original samples, and therefore it is preliminarily determined that the reduced order is 5.

[0052] First five columns of a left singular value vector matrix W are selected to construct a transformation matrix .sub.5, the order of the finite element model is reduced through the transformation matrix to obtain a 5-order reduced-order model, formula (5) is substituted into formula (13), and two ends of formula (13) are multiplied left by .sup.T.sub.5 to obtain the following formula:

[00016] 5 T K aa 5 b = 5 T ( g a - K ab q b ) ( 16 ) [0053] b is a coefficient vector, and the length is the reduced order. A displacement vector may be reconstructed according to formula (5).

[0054] Under the same displacement boundary condition, the displacement solution of the joint is obtained by calculating the reduced-order model and the finite element model separately, and a relative error of each degree of freedom in the displacement solution of the joint of the reduced-order model is calculated. Since numerical values of most of relative errors are close, in order to clearly express the relative errors of the reduced-order model, relative error logarithmic values are used. The definition is as follows:

[00017] err i = log 1 0 ( .Math. "\[LeftBracketingBar]" q i - q i ( 5 ) q i .Math. "\[RightBracketingBar]" ) ( 17 )

[0055] err.sub.i is a relative error logarithmic value of the reduced-order model on an i-th degree of freedom; and q.sub.i is a displacement value corresponding to the i-th degree of freedom in the displacement solution of the joint calculated by the finite element model, and q.sub.i(5) is a displacement value corresponding to the i-th degree of freedom in the displacement solution of the joint calculated by the 5-order reduced-order model. All degrees of freedom in the displacement solution of the joint are calculated, and calculation results are drawn in FIG. 5. It may be seen from FIG. 5 that the relative error logarithmic value is concentrated around 10, and a maximum relative error logarithmic value does not exceed 7. It may be seen through calculation that a maximum relative error is 0.8210.sup.7, which satisfies the requirements of practical engineering application. Therefore, the reduced order of the tower is determined as a 5 order.

[0056] Step 3: a schematic structural diagram of the tower is shown in (a) of FIG. 6. A total length of the physical entity is 1 m. An inclination sensor WT61C is installed at the top end of the tower, and may detect an inclination angle of the end surface of the top end of the tower; a strain sensor BF350-3EB is installed at a position 0.5 m away from the ground on a tower body, and may measure an average strain at an installation position of the tower and output a voltage; and the data of the strain sensor and the inclination sensor is collected through the STM32 single chip microcomputer, the collected data is uploaded to the computer through a serial communication port of the single chip microcomputer, and the computer processes the data and calculates the stress and the strain of the tower. The finally established physical entity of the tower is shown in (b) of FIG. 6, and a digital twinning model is shown in (c) of FIG. 6.

[0057] An inclination angle of an end surface of a top end of the tower is detected through an inclination sensor, and the displacement boundary condition is obtained by combining formula (1) with formula (2); and a voltage value of a strain sensor is sampled and read through an analog-to-digital converter (ADC), and an average strain of the physical entity at an installation position of the strain sensor may be calculated according to formula (18).

[00018] e 0 = E 4 ( R 1 R 1 + R 3 R 3 - R 2 R 2 - R 4 R 4 ) = E K s ( 18 )

[0058] e.sub.0 is an output voltage of the strain sensor, and has a unit of V; E is a supply voltage of the strain sensor, and has a unit of V; R.sub.1, R.sub.2, R.sub.3 and R.sub.4 are bridge resistors of the strain sensor, and have a unit of ; K.sub.s is a sensitivity coefficient of the strain sensor, which is related to a material used by the strain sensor, and the sensitivity coefficient of the BF350-3EB strain sensor is 2.1; and is an average strain of the physical entity at the installation position of the strain sensor.

[0059] A von-mises stress at the joint may be calculated through six stress components at the unit joint on the basis of the following formula:

[00019] von - mises = 2 2 ( xx - yy ) 2 + ( yy - ) 2 + ( - xx ) 2 + 6 ( xy 2 + y 2 + x 2 ) ( 19 )

[0060] .sub.xx, .sub.yy, .sub.zz, .sub.xy, .sub.yz and .sub.zx are six stress components of the stress at the unit joint separately, related definitions of which have been given in formula (14).

[0061] When the physical entity of the tower and the reduced-order model on the computer are interactively analyzed, the data of the inclination sensor and the strain sensor is read through an STM32 single chip microcomputer; the data is obtained, and then the data is uploaded to the computer through the STM32 single chip microcomputer; inclination data is processed and converted into the displacement boundary condition in the upper computer software that has been compiled, the displacement boundary condition is substituted into the 5-order reduced-order model that has been deployed in the upper computer to obtain the displacement solution of the joint, an interface of the upper computer software is shown in FIG. 7, and the interface is divided into three blocks, i.e. a tower three-dimensional model display area, a real-time data display area, and a communication and data display control function area; the von-mises stress at each joint is calculated through the displacement solution of the joint, results are drawn in the tower three-dimensional model display area through three-dimensional nephogram, a maximum joint von-mises stress value is displayed in real time in the data display area, and results are shown in FIG. 8; an average strain of the physical entity of the tower at the installation position of the strain sensor is calculated through the displacement solution of the joint, a deformed tower model is drawn in the three-dimensional model display area of the tower through the displacement solution of the joint, a measured value of the strain sensor uploaded by the STM32 single chip microcomputer and a calculated value of the reduced-order model are drawn in the same data visualization control in the real-time data display area, and moreover, a relative error of the reduced-order model is calculated, the relative error is displayed in another data visualization control in the real-time data display area, and strain online monitoring results of the simplified tower are shown in FIG. 9.

[0062] It may be seen from FIG. 8 that a maximum von-mises stress curve is changed gently and data does not have discontinuity and jump, which conforms to physical reality; and it may be seen from FIG. 9 that although a certain error exists between a calculated strain value of the reduced-order model and a measured strain value of the strain sensor, a change rule of the calculated strain value is consistent with a change rule of the measured strain value.

[0063] The present disclosure has the beneficial effects as follows: [0064] the present disclosure provides the digital twinning method for monitoring an operation state of a tower of a wind turbine generator system online. Compared with the prior art, the method can rapidly obtain the current state of the tower of the wind turbine generator system and visualize the current state of the tower of the wind turbine generator system in a manner of stress nephogram and a strain curve, and an analysis result has high precision, thereby achieving the purpose of monitoring the operation state of the tower online.

[0065] To sum up, the digital twinning method based on the POD can effectively monitor the operation state of the tower of the wind turbine generator system online.