Dynamic icing model construction method under three-dimensional rotating motion of fan blade

Abstract

The invention provides a method for constructing a dynamic icing model under three-dimensional rotating motion of a fan blade, which comprises the following steps of: constructing a three-dimensional blade model, and fitting a Lagrange point on the three-dimensional blade model by adopting an RBF (Radial Basis Function) to obtain a three-dimensional contour surface function; obtaining a plurality of two-dimensional section models of the three-dimensional blade model, and for each two-dimensional section model, fitting a Lagrange point of the edge of the two-dimensional section model by adopting RBF to obtain a corresponding two-dimensional contour surface function; and for each two-dimensional contour surface function, adopting Level-set to construct a corresponding dynamic section icing model, then based on the three-dimensional contour surface function, aggregating all the dynamic section icing models to obtain a three-dimensional icing model, and describing ice layer characteristics in the blade icing process according to the dynamic icing model. The method solves the problems that in the prior art, when a model test method is used for fitting the icing state of the blade, static experiment and empirical data are depended, so that experiment data are not coherent, the fitting effect is poor, and the characteristics of the icing process of the blade are difficult to dynamically restore.

Description

Technical Field

[0001] The present invention relates to the technical field of computational gas dynamics, and in particular to a method for constructing a dynamic icing model of a fan blade under three-dimensional rotational motion. Background Art

[0002] In recent years, high mountain wind farms in northwest my country have developed rapidly. Due to their special geographical location, these places are susceptible to icing weather due to factors such as humidity, altitude, and temperature, which brings about the problem of icing on wind turbine blades. Blade icing leads to frequent accidents in wind farms, which greatly affects the normal operation and power generation efficiency of wind farms. In mountain wind turbines, long-term icing on blades will change the aerodynamic shape of the blades, increase surface roughness, increase the drag coefficient, reduce the lift coefficient, and seriously damage the aerodynamic performance of the blades in extreme environments. Another thing worth noting is that different forms of icing will lead to uneven mass distribution of the entire impeller. The additional periodic centrifugal force generated by icing intensifies the vibration of the entire machine components, increases fatigue loads, and shortens the service life of wind power equipment. It can be seen that the problem of icing on wind turbine blades has become an important factor restricting the development and operation of wind power in icing areas.

[0003] Since the wind turbine icing problem is very complicated, the physical models involved have problems such as multi-phase coupling, strong discontinuity and difficult convergence. Traditional research methods mostly use model experiments. This method uses multiple sets of blade icing data under different static states to fit the blade icing state. It not only relies on experimental and empirical data, but also has high requirements for site and equipment. In addition, due to the non-dynamic fitting, the experimental data is incoherent, the fitting effect is poor, and it is difficult to dynamically restore the characteristics of the blade icing process. Summary of the invention

[0004] In view of the shortcomings existing in the prior art, the present invention provides a method for constructing a dynamic icing model of wind turbine blades under three-dimensional rotational motion, which solves the problem that the model test method in the prior art relies on static experimental and empirical data when fitting the blade icing state, resulting in incoherent experimental data, poor fitting effect, and difficulty in dynamically restoring the characteristics of the blade icing process.

[0005] According to an embodiment of the present invention, a method for constructing a dynamic icing model of a wind turbine blade under three-dimensional rotational motion comprises:

[0006] A three-dimensional blade model is constructed, and the Lagrange points on the three-dimensional blade model are fitted using RBF to obtain a three-dimensional isosurface function.

[0007] Acquire multiple two-dimensional cross-sectional models of the three-dimensional blade model, and for each two-dimensional cross-sectional model, use RBF to fit the Lagrange points at the edge of the two-dimensional cross-sectional model to obtain a corresponding two-dimensional isosurface function;

[0008] For each two-dimensional isosurface function, Level-set is used to construct the corresponding dynamic cross-section icing model. Then, based on the three-dimensional isosurface function, all dynamic cross-section icing models are aggregated to obtain a three-dimensional icing model. The characteristics of the ice layer during the blade icing process are described according to the dynamic icing model.

[0009] Preferably, the method of using RBF to fit the Lagrange points on the three-dimensional blade model to obtain the three-dimensional isosurface function includes:

[0010] The three-dimensional blade model is converted into a Lagrange point representation for representing the interface position information, and the field function is constructed using RBF based on all Lagrange points;

[0011] Based on the distance between two adjacent Lagrange points, a three-dimensional isosurface function is constructed.

[0012] Preferably, the field function is as follows:

[0013]

[0014] The three-dimensional isosurface function is as follows:

[0015] ψ(x)=-[d(x)-λ]

[0016] Where M is the total number of Lagrange points, (x i ,y,z) are the coordinates of any Lagrange point, (x j ,y j ,z j ) is the coordinate of the interpolation point of the RBF isosurface function, β j is the weight coefficient of RBF, R(||xx j ||) is the kernel function.

[0017] Preferably, the method of using RBF to fit the Lagrange points at the edge of the two-dimensional cross-sectional model to obtain the corresponding two-dimensional isosurface function is as follows:

[0018] The finite element method is used to construct a finite element mesh of a two-dimensional cross-section model, and the finite element equation is solved to obtain the node information of the finite element nodes on the finite element mesh;

[0019] Based on the rigid interface assumption of the beam model and the node information on the finite element nodes, the flow field variables of the Lagrangian points at the edge of the two-dimensional cross-sectional model are calculated;

[0020] Based on the flow field variables of the Lagrangian points at the edge of the two-dimensional cross-sectional model, RBF is used for fitting to obtain a two-dimensional isosurface function.

[0021] Preferably, after obtaining the two-dimensional isosurface function, the two-dimensional cross-sectional model is divided into a plurality of grid regions, and the grid regions are divided into gas grids, solid grids and virtual grids according to the two-dimensional isosurface function, and then the flow field variables of the virtual grids are reconstructed using the RBF difference format, and the flow field variables of the virtual grids of all the two-dimensional cross-sectional models are used as the boundary conditions of the two-dimensional isosurface function.

[0022] Preferably, the steps of constructing multiple dynamic cross-section icing models using the Level-set method include:

[0023] Construct the Level-set signed distance function Φ(x,y;t), and use the coordinates of all Lagrange points in the two-dimensional isosurface function as the initial values ​​of the signed distance function;

[0024] Solve the partial differential governing equation ▽Φ of the signed distance function Φ(x,y;t) to obtain the ice cross-section position corresponding to each Lagrangian point;

[0025] All ice-covered section positions are fitted to obtain a dynamic section ice-covered model.

[0026] Preferably, the interface region is divided according to the signed distance function Φ(x, y; t), and the partial differential control equation corresponding to each Lagrangian point is It is also necessary to determine whether it is satisfied If it is not satisfied, the signed distance function Φ(x, y; t) is reinitialized and the mass conservation correction is performed on the interface area.

[0027] Preferably, the reinitialization equation is as follows:

[0028]

[0029] Among them, τ is the virtual time step.

[0030] Compared with the prior art, the present invention has the following beneficial effects:

[0031] RBF is used to fit the Lagrangian points on the three-dimensional blade model of the wind turbine blade, and a three-dimensional isosurface function is generated to globally define the position of the blade during the dynamic rotation process. In addition, RBF is used to fit the Lagrangian points on the edge of the two-dimensional cross-sectional model of the blade, and a two-dimensional isosurface function is generated to globally define the deformation of the blade during the dynamic rotation process. Then, Level-set is used to construct the dynamic process of icing on the blade surface to describe the ice layer characteristics of the blade icing process. Since the entire fitting is performed under dynamic conditions, the data is continuous and not missing, which can ensure a better fitting effect and can better dynamically restore the characteristics of the blade icing process. BRIEF DESCRIPTION OF THE DRAWINGS

[0032] Figure 1 2 is a diagram of a prediction method according to an embodiment of the present invention.

[0033] Figure 2 It is a two-dimensional cross-sectional model diagram of an embodiment of the present invention.

[0034] Figure 3 This is a grid phase state identification diagram of an embodiment of the present invention.

[0035] Figure 4 Graph showing the positional relationship between the grid and the blade cross section according to an embodiment of the present invention.

[0036] Figure 5 This is a diagram of a dynamic cross-section ice coating model according to an embodiment of the present invention.

[0037] Figure 6 This is a parametric characterization diagram of airfoil icing according to an embodiment of the present invention.

[0038] Figure 7 It is a schematic diagram of the Level-set ice cover interface tracking process according to an embodiment of the present invention. DETAILED DESCRIPTION

[0039] The technical solution of the present invention is further described below in conjunction with the accompanying drawings and embodiments.

[0040] like Figure 1 As shown, the embodiment of the present invention proposes a method for constructing a dynamic icing model of a wind turbine blade under three-dimensional rotational motion, comprising:

[0041] A three-dimensional blade model is constructed, and the Lagrange points on the three-dimensional blade model are fitted using RBF to obtain a three-dimensional isosurface function.

[0042] A three-dimensional blade model is constructed for the wind turbine blades. For any immersed three-dimensional blade interface of the three-dimensional blade model, the RBF isosurface function is introduced based on the Lagrange point representing the interface position information to construct the following field function:

[0043]

[0044] Where M is the total number of Lagrange points, (x i ,y,z) are the coordinates of any Lagrange point, (x j ,y j ,z j ) is the coordinate of the interpolation point of the RBF isosurface function, β j is the weight coefficient of RBF.

[0045] Linear polynomial P(x) = p1 + p2x + p3y + p4z, this linear polynomial is not a necessary term. R(||xx j ||) is the kernel function, here we use the MQ kernel function:

[0046] R(||r ij ||)={r 2 +C 2} 1 / 2

[0047] in The shape parameter C is a constant, and here C=0.

[0048] RBF function weight coefficient β j The following orthogonality conditions are met:

[0049]

[0050] In the field function d(x) of the immersed interface, each Lagrangian point has an equal distance λ (λ>0), λ is an arbitrary constant greater than zero, and λ=1. Therefore, according to the Lagrangian point M and its isoline distance λ, the following solution to the (M+3)-order linear equation system about the coefficient β can be obtained:

[0051]

[0052] in Represents the distance between two points on the object surface, superscript b represents the object surface point, and subscript i and j represent the labels of arbitrary Lagrange points. In the coefficient matrix [A], subscripts 1, 2...M represent the numbers of the Lagrange points, and superscript bb represents that both the field point and the point to be interpolated are Lagrange points. The MQ function is conditionally positive definite. After adding the polynomial basis, the coefficient matrix of equation (4) becomes a positive definite matrix, ensuring that equation (4) has a unique solution. The implicit interface description method of the RBF isosurface function can obtain higher accuracy with limited Lagrange points, and the isosurface equation (4) can be quickly solved by the traditional Gaussian elimination method. After obtaining the coefficient matrix, the obtained field function d(x) is the three-dimensional contour function.

[0053] The following three-dimensional isosurface function ψ(x) is further constructed to facilitate the identification of the properties of the grid node ψ(x) in the field.

[0054] ψ(x)=-[d(x)-λ]

[0055] The resulting ψ(x) has the following signed distance property: for any coordinate node x i,j,k , if ψ i,j,k >0, indicating point x i,j,k In the object; i,j,k <0 indicates point x i,j,k Outside the object; i,j,k =0 represents point x i,j,k On the physical plane.

[0056] At the initial moment, select any Lagrangian point on the blade and substitute it into the three-dimensional isosurface function to obtain the result ψ(x, y; t=0). During the rotation of the blade, the coordinates of the Lagrangian point on the blade obtained at any time are substituted into the three-dimensional isosurface function. If the obtained result ψ(x, y; t) is the same as ψ(x, y; t=0), it can be considered that the Lagrangian point and the Lagrangian point at the initial moment belong to the same blade, so as to globally track the position of the blade during the dynamic rotation process.

[0057] Acquire multiple two-dimensional cross-sectional models of the three-dimensional blade model, and for each two-dimensional cross-sectional model, use RBF to fit the Lagrange points at the edge of the two-dimensional cross-sectional model to obtain a corresponding two-dimensional isosurface function;

[0058] First, during the rotation of the blade, it will be affected by external forces and produce slight deformation. Therefore, the shape of the two-dimensional cross-section of the three-dimensional blade model will also change continuously during the rotation of the blade. Therefore, the two-dimensional cross-sectional model of the blade is also dynamically changing. Therefore, RBF can also be used to track the cross-sectional shape of the blade in real time. At this time, it is necessary to consider the acceleration, velocity and displacement information of the Lagrange point of the two-dimensional cross-sectional model to be transmitted to the virtual grid to consider the effect of the flexible body on the flow field.

[0059] The finite element method is used to construct the finite element mesh of the two-dimensional cross-section model. The finite element mesh is usually constructed on the central axis of the blade section. By solving the finite element equation, the node information at the finite element nodes stored on the finite element mesh can be obtained. Subsequently, based on the rigid section assumption of the beam model, the finite element node information is transferred to the Lagrangian point on the blade section surface.

[0060] The acceleration, velocity and displacement of the finite element nodes at the current moment are known. For any finite element nodes 1, 2, 3 and 4, their node coordinates are (x1, y1), (x2, y2), (x3, y3), (x4, y4), respectively. Figure 2 The blade two-dimensional cross-section model in . Consider the cross section of the flexible beam with a thickness of h. Since the cross section of the central axis is always a plane before and after deformation, taking unit node 2 as an example, the interface points W2 and E2 on the beam boundary symmetrical about node 2 can be expressed as:

[0061]

[0062] (n x ,n y ) is the normal vector of point 2, and its expression is:

[0063]

[0064] Due to the rigid cross-section assumption, the velocity and acceleration of finite element node 2 are equal to the velocity and acceleration of the Lagrange points W2 and E2 at the edge of the blade section corresponding to it:

[0065] u W2 =u2,u E2 =u2

[0066] a W2 =a2,a E2 =a2

[0067] Then, referring to the above method, the acceleration, velocity and displacement information of a series of Lagrangian points constituting the boundary of the two-dimensional cross-sectional model are calculated.

[0068] Then, the above-mentioned RBF is used to fit the Lagrangian points of the flexible blade cross-sectional surface to construct a two-dimensional isosurface function of the two-dimensional cross-sectional model. In this way, during the rotation of the blade, no matter what kind of deformation the blade undergoes, it can be judged that the blade belongs to the same blade as the blade at the initial moment. Then, according to the two-dimensional isosurface function, the immersed boundary virtual grid method is used to divide the two-dimensional cross-sectional model into multiple grids, and the virtual grids, gas grids and solid grids are identified and located.

[0069] According to the two-dimensional isosurface function, it is necessary to distinguish between virtual grid (GC), gas grid (AC) and solid grid (SC), such as Figure 3 As shown, the air domain variable values ​​of the background virtual grid are then reconstructed to impose object surface boundary conditions.

[0070] After constructing the isosurface function of the object surface, the next step is to determine whether the Lagrange point is inside or outside the object, so as to identify the background grid properties, including gas grid (AC), solid grid (SC) and virtual grid (GC). Gas grids are defined as those grids whose grid cell centers are located in the gas domain; virtual grids refer to those grids whose grid cell centers are located in the fixed domain but adjacent to the gas domain; solid grids refer to those grids whose grid cell centers are located in the solid domain but exclude virtual grids. The traditional identification method is the ray method, but the numerical experiments mentioned in the patent show that the ray method may fail when dealing with irregular boundaries, resulting in the grid in the gas domain being mistakenly identified as a virtual grid. Therefore, the present invention proposes a direct identification method based on the sign characteristics of the distance function (RBF), which effectively avoids such misjudgment. This method is applicable to any geometric interface.

[0071] High-order convection discrete formats such as TVD-MUSCL generally occupy multiple grids. In order to be consistent with the discrete format, the present invention applies two layers of virtual grids to impose boundaries and divides the two-dimensional cross-sectional model into multiple grid areas. The identification algorithm of the virtual grid, gas grid, and solid grid is as follows:

[0072] (1) If ψ i,j,k ≥0 and (ψ i-2,j,k ≥0 and ψ i-1,j,k ≥0 and ψ i+1,j,k ≥0 and ψ i+2,j,k ≥0 and ψ i,j-2,k ≥0 and ψ i,j-1,k ≥0 and ψ i,j+1,k ≥0 and ψ i,j+2,k ≥0 and ψ i,j,k-2 ≥0 and ψ i,j,k-1 ≥0 and ψ i,j,k+1 ≥0 and ψ i,j,k+2 ≥0); then it is a solid grid (SC).

[0073] (2) If ψ i,j,k ≥0 and (ψ i-2,j,k <0 or ψ i-1,j,k <0 or ψ i+1,j,k <0 or ψ i+2,j,k <0 or ψ i,j-2,k <0 or ψ i,j-1,k <0 or ψ i,j+1,k <0 or ψ i,j+2,k <0 or ψ i,j,k-2 <0 or ψ i,j,k-1 <0 or ψ i,j,k+1 <0 or ψ i,j,k+2 <0); it is a virtual grid (GC).

[0074] (3) If ψ i,j,k <0 and (ψ i-2,j,k <0 and ψ i-1,j,k <0 and ψ i+1,j,k <0 and ψ i+2,j,k <0 and ψ i,j-2,k <0 and ψ i,j-1,k <0 and ψ i,j+1,k <0 and ψ i,j+2,k <0 and ψ i,j,k-2 <0 and ψ i,j,k-1 <0 and ψ i,j,k+1 <0 and ψ i,j,k+2 <0); then it is the gas grid (AC).

[0075] After that, the RBF interpolation format is used to reconstruct the virtual grid variable values. Since the projection points of the Lagrangian points at the edge of the two-dimensional section model and the Lagrangian points of the three-dimensional blade model on the edge of the two-dimensional section model are inconsistent, the linear interpolation format is used here to calculate the flow field variables (velocity or acceleration) at the interface intersection of the two-dimensional section model and the three-dimensional blade model. Figure 4 Take the interface intersection (BI) in as an example:

[0076]

[0077] is the velocity component or acceleration component of the interface projection point. Finally, the boundary conditions of the blade section are applied by satisfying the following mirror conditions:

[0078]

[0079] For the Dirichlet boundary conditions that impose velocity, ζ1 = -1 and is the x-, y-, or z-velocity component of the projection point. Neumann boundary conditions for pressure variables The coefficients are calculated as ξ1 = 1 and ξ2 = Δd GM ·n Γ ·a BI , a BI is the acceleration vector of the object surface, n Γ is the boundary normal vector, Δd GM is the distance from the virtual grid point to the mirror point. Unlike rigid boundaries, for elastic boundaries, the gas variables corresponding to each virtual grid projection point are different.

[0080] Finally, in order to suppress the problem of severe pressure oscillation, the present invention develops an area fraction representation method for flexible moving boundaries, and for the intersection of the flexible boundary and the background grid, a method similar to the structural finite element shape function interpolation is used to calculate the intersection velocity.

[0081] For each two-dimensional isosurface function, Level-set is used to construct the corresponding dynamic cross-section icing model. Then, based on the three-dimensional isosurface function, all dynamic cross-section icing models are aggregated to obtain a three-dimensional icing model. The characteristics of the ice layer during the blade icing process are described according to the three-dimensional icing model.

[0082] The dynamic icing of blades requires the establishment of a dynamic cross-sectional icing model of blades in a two-dimensional cross-sectional shape in the three-dimensional length direction ( Figure 5 As shown in Figure 2, based on the data provided by the wind turbine test site under snowy weather in the mountainous area, a modified surface icing model is used to construct the dynamic cross-section icing characteristic equation to capture the icing morphology under different cross-sections. The Lagrangian points on the icing surface are fitted by the radial basis isosurface function to form an icing change isosurface function to dynamically track and capture the changes on the icing surface.

[0083] Streamlined icing is usually formed under rime icing conditions with low ambient temperature and small changes in wind speed and direction. Its shape is regular and evenly distributed along the span of the blade. Angular icing often occurs under rime icing conditions with high ambient temperature and large changes in wind speed and direction. It contains more irregular structures and is unevenly distributed along the span of the blade.

[0084] like Figure 6As shown in Figure 1, under normal circumstances, the best way to obtain icing data is through field observation, where the icing shape at the blade tip (r / R=0.99) can be recorded more accurately. Therefore, the icing shape of the airfoil at this location is used to approximate the icing condition of the entire blade. At the same time, naturally formed icing usually contains many tiny protrusions and burr-like structures. These features are difficult to form a fitted surface geometry mesh when meshing, which easily leads to a decrease in the quality of the boundary layer mesh and may cause difficulties in numerical calculation convergence. Therefore, in the modeling process, only the more significant icing contour structures are usually retained. By selecting appropriate contour points on the upper and lower surfaces of the airfoil in turn (especially for encryption near the stagnation point with large curvature), and successively calculating the icing thickness along the external normal direction corresponding to each contour point, the parametric characterization of the icing outer contour can be achieved.

[0085] After obtaining certain blade ice data, it is necessary to track and capture the ice layer on the three-dimensional wind turbine blades at all times. This invention adopts the classic Level-set interface capture method, calibrates the tracking interface as the zero isosurface of the high-dimensional function, solves the partial differential equation with high-order finite differences, and implicitly captures the motion topology interface:

[0086] The first step is to construct the Level-set signed distance function.

[0087] In two-dimensional space, the scalar Φ is defined as the signed distance function from the grid node to the motion interface, that is, the shortest distance from any Euler grid point in the calculation area to the interface, and it satisfies Φ<0 inside the interface and Φ>0 outside the interface. Therefore, the Level-set function value Φ(x,y)=0 on the two-dimensional section isosurface satisfies {(x,y)|(x,y)∈interface point set}. The initial time is Φ(x,y;t=0)=d(x,y), where d(x,y) is all Lagrangian points in the two-dimensional isosurface function, that is, the coordinates of all Lagrangian points are used as the initial value of the function.

[0088] The second step is to solve the partial differential governing equation of the Level-set signed distance function.

[0089] The signed distance function Φ(x,y; t=0) changes with time and space. By taking the time derivative of Φ(x,y; t), the motion control equation of the interface evolution is constructed. In order to improve the solution accuracy, a high-order upwind discretization format (such as ENO and WENO formats) is used, including the use of extended velocity methods in the computational domain.

[0090] The third step is to introduce the reinitialization process.

[0091] Judge the updated Level-set symbol distance function. If it does not meet the (Distance function properties), reinitialization is performed in the area around the interface.

[0092] The fourth step is to correct the quality of the area within the interface after reinitialization.

[0093] During the reinitialization process, the interface area is divided according to the signed distance function Φ(x, y; t), which causes a slight adjustment to the signed distance function around the interface, making it difficult to maintain consistent quality before and after the adjustment. The quality of the area within the interface is corrected to ensure quality conservation.

[0094] The fifth step is to output the interface position at the final moment and implement explicit processing of the implicit interface.

[0095] In summary, steps 1-4 are solved in a loop according to the signed distance function on the background Euler grid nodes to realize the evolution of the interface over time. The relevant numerical discretization process is as follows:

[0096] (1) The interface is characterized by the signed distance function, Γ(x, y; t) is {(x, y)|Φ(x, y; t) = 0, (x, y)∈Γ(x, y)}, and the entire region Ω is divided by the interface Γ(x, y; t). The area outside the interface is represented as {(x, y)|Φ(x, y; t) > 0}, and the area inside the interface is represented as {(x, y)|Φ(x, y; t) < 0};

[0097] (2) The unit normal vector at the interface Γ(x, y; t) is represented by the signed distance function It is expressed as:

[0098]

[0099] (3) The curvature κ at the interface Γ(x, y; t) is expressed as:

[0100]

[0101] (4) The interface reinitialization process refers to forcing the Level-set function Φ(x, y; t) to re-satisfy the distance

[0102] Function d(x, y; t), that is, the distance function is positive outside the interface and negative inside the interface. The interface reinitialization process can effectively smooth the interface and suppress the degree of gradient change of the function around the interface, so that d(x, y; t) satisfies at the interface Γ(x, y; t) The virtual time τ is introduced into the calculation (the time iteration step is selected to satisfy τ = O(ε)). In the virtual time step, only in the area with width ε around the interface (ε is generally 1.5 times the grid width), the interface reinitialization equation is solved, that is,

[0103]

[0104] Introducing the Delta Function Smooth the interface only around the interface. Extend the smoothing of interface quantities to neighboring mesh elements, such as extending the normal velocity v n Extending from the interface to the adjacent isosurface, by introducing the virtual time τ, the following equation is solved:

[0105]

[0106] v n (x,y;τ=0)=v n (x,y;t)

[0107] By using the above-mentioned level-set method and according to the three-dimensional isosurface function, multiple dynamic cross-section icing models are aggregated to obtain a three-dimensional icing model, such as Figure 7 As shown in the figure, the ice state of the blade surface at three consecutive moments is described, and the interface is characterized by the signed distance function. i ,y i ;t i ) is {(x i ,y i )Φ(x i ,y i ;t i )=0,(x i ,y i )∈Γ(x i ,y i )}, and the interface is represented as {(x i ,y i )Φ(x i ,y i ;t i )>0}, which is represented in the interface as {(x i ,y i )Φ(x i ,y i ;t i )<0}, an empirical parameter of ice accumulation is introduced. By correcting the mass of the area within the interface, the conservation of the ice mass on the airfoil surface is ensured. In this way, the ice layer shape characteristics and related data of the three-dimensional blade in the dynamic icing process can be captured and described in real time.

[0108] Finally, it should be noted that the above embodiments are only used to illustrate the technical solution of the present invention rather than to limit it. Although the present invention has been described in detail with reference to the preferred embodiments, those skilled in the art should understand that the technical solution of the present invention can be modified or replaced by equivalents without departing from the purpose and scope of the technical solution of the present invention, which should be included in the scope of the claims of the present invention.

Claims

1. A method for constructing a dynamic icing model of a wind turbine blade under three-dimensional rotational motion, characterized in that: include: A three-dimensional blade model is constructed, and the Lagrange points on the three-dimensional blade model are fitted using RBF to obtain a three-dimensional isosurface function. Acquire multiple two-dimensional cross-sectional models of the three-dimensional blade model, and for each two-dimensional cross-sectional model, use RBF to fit the Lagrange points at the edge of the two-dimensional cross-sectional model to obtain a corresponding two-dimensional isosurface function; For each two-dimensional isosurface function, Level-set is used to construct the corresponding dynamic cross-section icing model. Then, based on the three-dimensional isosurface function, all dynamic cross-section icing models are aggregated to obtain a three-dimensional icing model. The characteristics of the ice layer during the blade icing process are described according to the dynamic icing model.

2. The method for constructing a dynamic ice coating model of a wind turbine blade under three-dimensional rotational motion according to claim 1, characterized in that: Methods for obtaining three-dimensional isosurface functions by fitting the Lagrange points on a three-dimensional blade model using RBF include: The three-dimensional blade model is converted into a Lagrange point representation for representing the interface position information, and the field function is constructed using RBF based on all Lagrange points; Based on the distance between two adjacent Lagrange points, a three-dimensional isosurface function is constructed.

3. The method for constructing a dynamic ice coating model of a wind turbine blade under three-dimensional rotational motion according to claim 1, characterized in that: The field function is as follows: The three-dimensional isosurface function is as follows: ψ(x)=-[d(x)-λ] Where M is the total number of Lagrange points, (x i ,y,z) are the coordinates of any Lagrange point, (x j ,y j ,z j ) is the coordinate of the interpolation point of the RBF isosurface function, β j is the weight coefficient of RBF, R(||xx j ||) is the kernel function.

4. The method for constructing a dynamic ice coating model of a wind turbine blade under three-dimensional rotational motion according to claim 1, characterized in that: The method of using RBF to fit the Lagrange points at the edge of a two-dimensional cross-section model to obtain the corresponding two-dimensional isosurface function is as follows: The finite element method is used to construct a finite element mesh of a two-dimensional cross-section model, and the finite element equation is solved to obtain the node information of the finite element nodes on the finite element mesh; Based on the rigid interface assumption of the beam model and the node information on the finite element nodes, the flow field variables of the Lagrangian points at the edge of the two-dimensional cross-sectional model are calculated; Based on the flow field variables of the Lagrangian points at the edge of the two-dimensional cross-sectional model, RBF is used for fitting to obtain a two-dimensional isosurface function.

5. The method for constructing a dynamic ice coating model of a wind turbine blade under three-dimensional rotational motion according to claim 1, characterized in that: After obtaining the two-dimensional isosurface function, the two-dimensional cross-sectional model is divided into multiple grid areas. The grid areas are divided into gas grids, solid grids and virtual grids according to the two-dimensional isosurface function. Then, the flow field variables of the virtual grids are reconstructed using the RBF difference format, and the flow field variables of the virtual grids of all two-dimensional cross-sectional models are used as the boundary conditions of the two-dimensional isosurface function.

6. The method for constructing a dynamic ice coating model of a wind turbine blade under three-dimensional rotational motion according to claim 1, characterized in that: The steps of constructing multiple dynamic cross-section icing models using the Level-set method include: Construct the Level-set signed distance function Φ(x,y;t), and use the coordinates of all Lagrange points in the two-dimensional isosurface function as the initial values ​​of the signed distance function; Solve the partial differential governing equation ▽Φ of the signed distance function Φ(x,y;t) to obtain the ice cross-section position corresponding to each Lagrangian point; All ice-covered section positions are fitted to obtain a dynamic section ice-covered model.

7. The method for constructing a dynamic ice coating model of a wind turbine blade under three-dimensional rotational motion according to claim 1, characterized in that: The interface area is divided according to the signed distance function Φ(x, y; t). For the partial differential control equation ▽Φ corresponding to each Lagrangian point, it is also necessary to determine whether |▽Φ|=1 is satisfied. If not, the signed distance function Φ(x, y; t) is reinitialized and a mass conservation correction is performed on the interface region.

8. The method for constructing a dynamic ice coating model of a wind turbine blade under three-dimensional rotational motion according to claim 1, characterized in that: Reinitialize the equation as follows: Among them, τ is the virtual time step.

Images