Yaw bearing dynamics simulation method and system

Abstract

The invention discloses a yaw bearing dynamics simulation method and system, and belongs to the technical field of simulation, and the yaw bearing dynamics simulation system comprises a calculation module, a modeling module, a first solving module, a determination module, a second solving module and a construction module. By combining the three-dimensional power spectral density of turbulent wind, multi-system effects such as pneumatic and structural effects of a fan are coupled, the comprehensive exciting force of a yaw bearing is obtained, finally, the vibration response of the four-point contact ball bearing under the random wind load is calculated by adopting a dynamic modeling technology considering a contact judgment criterion, and the yaw bearing vibration response is calculated. The turbulence characteristic of the wind speed, the dynamic transmission process of the load and the coupling effect among the modules of the fan are comprehensively considered, so that the simulation precision is improved.

Description

Technical Field

[0001] The invention belongs to the technical field of simulation, and in particular relates to a yaw bearing dynamics simulation method and system. Background Art

[0002] In wind power generation systems, wind turbine yaw bearings, as key components, are subject to complex random wind loads. These loads are transmitted to the yaw bearings layer by layer through the wind rotor, hub, nacelle and tower. However, due to the complex and changeable turbulent characteristics of wind speed and the dynamic transfer process of loads, the existing dynamic analysis of wind turbine yaw bearings is mostly based on simplified models or deterministic wind loads, which cannot accurately simulate the dynamic response under random wind loads. At present, some studies have attempted to improve the simulation accuracy by introducing random wind load models, but most of these models fail to fully consider the turbulent characteristics of wind speed and its correlation in space and time. At the same time, these studies often ignore the coupling effects between the various modules of the wind turbine, resulting in a large deviation between the simulation results and the actual situation. Summary of the invention

[0003] The object of the present invention is to provide a yaw bearing dynamics simulation method to solve the problem of low simulation accuracy of the existing yaw bearing mentioned in the above background technology.

[0004] To achieve the above object, the present invention provides the following technical solution: a yaw bearing dynamics simulation method, comprising:

[0005] Obtain the power spectrum density function of any position point in the wind field in different directions of the coordinate system, and calculate the wind speed time series of any position point in the wind field in different directions based on the power spectrum density function;

[0006] Based on the wind speed time series, the blade aerodynamic load is calculated through the iterative calculation of blade element momentum theory, and the tower aerodynamic load is obtained based on the blade aerodynamic load;

[0007] Model the tower and solve the elastic force of the tower based on the tower model;

[0008] Construct the dynamic equation of the tower based on inertial force, elastic force and aerodynamic load, and solve the inertial force of the tower based on the dynamic equation of the tower;

[0009] The yaw bearing boundary conditions are obtained based on the tower inertia force;

[0010] Construct a four-point contact model of the yaw bearing and solve the contact load based on the four-point contact model of the yaw bearing;

[0011] The dynamic equation of the yaw bearing is constructed based on the boundary conditions and contact loads of the yaw bearing.

[0012] Furthermore, the obtaining of the power spectrum density function of any position point in the wind field in different directions in the coordinate system includes:

[0013] Construct spatial correlation functions in different directions of the coordinate system, and calculate the integral scales in different directions of the coordinate system based on the spatial functions in different directions;

[0014] The power spectral density function in different directions of the coordinate system is calculated based on the integral scale in different directions of the coordinate system.

[0015] Furthermore, the spatial correlation function in any direction on the coordinate system is as follows:

[0016]

[0017] Where i is the wind speed point to be calculated, j is k is the set of points that will produce wind speed correlation with the calculation point, r is the distance between two points in the selected direction, unit m, a k is the coherence intensity control coefficient between two points, is the set value, b k is the distance coherence intensity control coefficient between two points, is the set value, z m is the average height between the two points, f is the frequency in Hz, and U is the average wind speed between the two points in m / s.

[0018] Furthermore, the integral scale calculation formula in any direction in the coordinate system is as follows:

[0019]

[0020] r represents the distance between two points in any direction of the coordinate system, in meters.

[0021] Furthermore, the calculation formula of the power spectral density function in any direction of the coordinate system is as follows:

[0022]

[0023] Where S(f) represents the power spectral density in the selected direction, σ represents the standard deviation of wind speed in the selected direction, L represents the integral in the selected direction, U is the average wind speed in m / s, and f is the frequency in Hz.

[0024] Furthermore, the blade aerodynamic load is obtained by iterative calculation based on the wind speed time series through blade element momentum theory, including:

[0025] Given the initial value of the axial induction factor, the tangential induction factor is calculated;

[0026] The blade lift and drag calculation formula 1 based on the lift coefficient and drag coefficient is constructed through the blade element theory;

[0027] The second blade lift and drag calculation formula based on the axial induction factor and the tangential induction factor is constructed through momentum theory;

[0028] Based on the equal relationship between blade lift and drag calculation formula 1 and calculation formula 2, the lift coefficient and drag coefficient are calculated by substituting the axial induction factor and the tangential induction factor;

[0029] Construct a calculation formula for the axial induction factor based on the lift coefficient and the drag coefficient;

[0030] Substitute the calculation formula of the axial induction factor based on the lift coefficient and the drag coefficient into the lift coefficient and the drag coefficient to obtain the iterative value of the axial induction factor;

[0031] The convergence of the axial induction factor is judged. When the axial induction factor converges, the aerodynamic load of the blade is calculated using the current iterative value of the axial induction factor.

[0032] Furthermore, the calculation formula of the axial induction factor based on the lift coefficient and the drag coefficient is as follows:

[0033]

[0034] In the formula represents the inflow angle, C l is the lift coefficient, C d is the drag coefficient, a is the axial induction factor, σ′ is the solidity of the wind wheel, and the calculation formula is as follows:

[0035] σ′=Bc / C

[0036] Where C is the circumference of the wind wheel, c is the local chord length of the blade, and B is the number of blades.

[0037] Furthermore, the yaw bearing includes a rolling element and a ring having an inner ring and an outer ring, and solving the contact load based on the yaw bearing four-point contact model includes:

[0038] The contact geometry approach between the rolling element and the inner ring is calculated, and the contact load is calculated based on the contact geometry approach.

[0039] Furthermore, the formula for calculating the contact load by contact geometry approach is as follows:

[0040] Q=Kδ 1.5

[0041] Where: K is the Hertz contact coefficient, Q is the contact load, and δ is the contact geometry approach.

[0042] On the other hand, the present application discloses a yaw bearing dynamics simulation system, comprising:

[0043] An acquisition module is configured to acquire a power spectrum density function of any position point in the wind farm in different directions of the coordinate system, and calculate a wind speed time series of any position point in the wind farm in different directions based on the power spectrum density function;

[0044] A calculation module is configured to obtain a blade aerodynamic load by iterative calculation based on a wind speed time series through blade element momentum theory, and obtain a tower aerodynamic load based on the blade aerodynamic load;

[0045] A modeling module configured to model the tower and solve the tower elastic force based on the tower model;

[0046] A first solving module is configured to construct a dynamic equation of the tower based on inertial force, elastic force and aerodynamic load, and solve the inertial force of the tower based on the dynamic equation of the tower;

[0047] a determination module configured to obtain a yaw bearing boundary condition based on a tower inertia force;

[0048] A second solving module is configured to construct a four-point contact model of the yaw bearing and solve the contact load based on the four-point contact model of the yaw bearing;

[0049] A building module is configured to build the dynamic equations of the yaw bearing based on the yaw bearing boundary conditions and contact loads.

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

[0051] This application combines the three-dimensional power spectral density of turbulent wind, couples the aerodynamic, structural and other multi-system effects of the wind turbine, and obtains the comprehensive excitation force of the yaw bearing. Finally, the dynamic modeling technology considering the contact judgment criterion is used to calculate the vibration response of the four-point contact ball bearing under random wind load. The simulation accuracy is improved by comprehensively considering the turbulent characteristics of the wind speed and the dynamic transfer process of the load, as well as the coupling effect between the various modules of the wind turbine. BRIEF DESCRIPTION OF THE DRAWINGS

[0052] Figure 1 This is a flow chart of the application method;

[0053] Figure 2 It is the time series diagram of wind speed in X direction of the wind turbine;

[0054] Figure 3 Schematic diagram of the wind wheel and blade model Figure 1 ;

[0055] Figure 4 Schematic diagram of the wind wheel and blade model Figure 2 ;

[0056] Figure 5 Schematic diagram of iterative calculation of blade aerodynamic loads;

[0057] Figure 6 The aerodynamic load diagram of the yaw bearing after raindrop counting;

[0058] Figure 7 Contact analysis diagram of a four-point contact ball bearing. DETAILED DESCRIPTION

[0059] The following will be combined with the drawings in the embodiments of the present invention to clearly and completely describe the technical solutions in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0060] A yaw bearing dynamics simulation method, referring to Figure 1 ,include:

[0061] S100: obtaining a power spectrum density function of any position point in the wind field in different directions of the coordinate system, and calculating a wind speed time series of any position point in the wind field in different directions based on the power spectrum density function;

[0062] S200: Based on the wind speed time series, the blade aerodynamic load is calculated through the blade element momentum theory iteration, and the tower aerodynamic load is obtained based on the blade aerodynamic load;

[0063] S300: Modeling the tower and solving the elastic force of the tower based on the tower model;

[0064] S400: construct the dynamic equation of the tower based on inertial force, elastic force and aerodynamic load, and solve the inertial force of the tower based on the dynamic equation of the tower;

[0065] S500: obtaining a yaw bearing boundary condition based on the tower inertia force;

[0066] S600: Construct a four-point contact model for the bearing and calculate the contact load based on the four-point contact model for the bearing;

[0067] S700: Construct the dynamic equations of the yaw bearing based on the yaw bearing boundary conditions and contact loads.

[0068] In S100, the specific steps of obtaining the power spectrum density function of any position point in the wind farm in different directions in the coordinate system are as follows:

[0069] S101: constructing spatial correlation functions in different directions (e.g., X, Y, and Z axis directions in the inertial coordinate system) of a coordinate system (e.g., an inertial coordinate system), and calculating integral scales in different directions of the coordinate system based on the spatial functions in different directions;

[0070] S102: Calculate the power spectral density function in different directions of the coordinate system based on the integral scale in different directions of the coordinate system.

[0071] Specifically, the spatial correlation function in any direction on the coordinate system is as follows:

[0072]

[0073] In formula (1), i is the wind speed point to be calculated, j is k is the set of points that will produce wind speed correlation with the calculation point, r is the distance between the two points in the selected direction, unit is m; generally speaking, the smaller r is, that is, the closer the distance is, the stronger the correlation between the wind speeds tends to be. k is the coherence strength control coefficient between two points, and its size directly affects the strength of wind speed correlation. It is the set value, b k is the distance coherence intensity control coefficient between two points, which controls the influence of distance on coherence intensity and is the set value, z m is the average height between the two points, reflecting the effect of the difference in spatial position in the vertical direction on the wind speed correlation. f is the frequency, which reveals the influence of the characteristics of wind speed change on the correlation from another dimension. The unit is Hz. U is the average wind speed between the two points. The unit is m / s, which provides a benchmark reference for evaluating the wind speed correlation.

[0074] After obtaining the spatial correlation function in any direction of the coordinate system through the above formula (1), the integral scale of the corresponding direction in the coordinate system is calculated through the following formula (2):

[0075]

[0076] In formula (2), r represents the distance between two points in any direction of the coordinate system, with the unit of m.

[0077] After obtaining the integral scale and spatial correlation function, the power spectral density function in different directions in the coordinate system is calculated by the following formula:

[0078]

[0079] Among them, S u (f) S v (f) S w (f) represents the power spectrum density along the wind direction (u), across the wind direction (v) and perpendicular to the wind direction (w), σ u , σ v , σ w Represents the standard deviation of wind speed in each direction, L u , L v , L wThey respectively indicate that the integral scales in each direction are related to the spatial coherence, U is the average wind speed in m / s, and f is the frequency in Hz.

[0080] In some embodiments, the formula (3) for calculating the wind speed time series of any point in the wind field in different directions based on the power spectrum density function is as follows:

[0081]

[0082] In formula (3), U is the average wind speed, in m / s, S(f n ) is the wind speed power spectrum density, is the period, φ n is the phase angle, Figure 2 An example of a wind speed time series diagram with an average wind speed of 8 m / s in the axial direction is given.

[0083] In S200, refer to Figure 5 , and combined with Figure 3 and Figure 4 The schematic diagram of the wind wheel and blade model is shown in Figure 1. Based on the wind speed time series, the specific steps for iteratively calculating the blade aerodynamic load through the blade element momentum theory are as follows:

[0084] First, the initial value of the axial induction factor is given, and the tangential induction factor is calculated;

[0085] The axial induction factor represents the attenuation rate of the axial wind speed, and the initial value can be given as 0, for example, and the definition formula is as follows:

[0086]

[0087] Where a represents the axial induction factor, V ∞ is the axial wind speed at infinity, in m / s, and V is the average axial wind speed, in m / s.

[0088] The tangential induction factor represents the ratio of the tangential wind speed at the rotor plane to the blade speed. The calculation formula is as follows:

[0089]

[0090] Where b is the tangential induction factor, Ω is the impeller angular velocity, in rad / s, is the axial wind speed at the rotor plane, in m / s, and r is the blade length, in m.

[0091] The blade lift and drag calculation formula based on lift coefficient and drag coefficient is constructed through blade element theory as follows:

[0092]

[0093] Formula (4) is the blade lift calculation formula 1, and formula (5) is the blade resistance calculation formula 2. In formulas (4) and (5), W is the relative wind speed at the rotor plane, in m / s, V ∞ is the axial wind speed at infinity, in m / s, is the axial wind speed at the rotor plane, in m / s, B is the number of blades, and ρ is the air density, in kg / m 3 , C l is the lift coefficient, C d is the drag coefficient, c is the local chord length of the blade, in m, dr is the length of the blade's tiny section, represents the inflow angle.

[0094] The blade lift and drag calculation formula 2 based on the axial induction factor and the tangential induction factor is constructed through momentum theory as follows:

[0095] dT=4πrρV ∞ 2 (1-a)adr (6)

[0096] dQ=4πr 3 ρV ∞ Ω(1-a)bdr (7)

[0097] Formula (6) is the blade lift calculation formula 2, and formula (7) is the blade resistance calculation formula 2. In formulas (6) and (7), r is the blade length, unit m, and ρ is the air density, unit Kg / m 3 , V ∞ is the axial wind speed at infinity, in m / s, a is the axial induction factor, and b is the tangential induction factor.

[0098] Based on the equality between blade lift and drag calculation formula 1 and calculation formula 2, that is, the equality between formula (4) and (6), formula (5) and (7), the axial induction factor (starting as the initial value) and the tangential induction factor (that is, the above calculated value) are substituted to calculate the lift coefficient and the drag coefficient;

[0099] The calculation formula (8) for the axial induction factor based on the lift coefficient and the drag coefficient is as follows:

[0100]

[0101] In formula (8), represents the inflow angle, C l is the lift coefficient, C d is the drag coefficient, a is the axial induction factor, σ′ is the solidity of the wind wheel, and the calculation formula is as follows:

[0102] σ′=Bc / C

[0103] Where C is the circumference of the wind wheel, c is the local chord length of the blade, and B is the number of blades.

[0104] Substituting the lift coefficient and drag coefficient calculated above into the calculation formula (8) of the axial induction factor based on the lift coefficient and the drag coefficient, the iterative value of the axial induction factor can be obtained.

[0105] The convergence of the iterative value of the axial induction factor is judged. When the axial induction factor converges, the aerodynamic load of the blade is calculated using the current axial induction factor and the iterative values ​​of the lift coefficient and the drag coefficient.

[0106] The convergence of the axial induction factor is determined based on the difference between the iterative value obtained by calculating the axial induction factor and the axial induction factor value input at the beginning (for example, a given initial value). When the difference between the iterative value obtained by calculating the axial induction factor and the axial induction factor value input at the beginning (for example, a given initial value) is less than a set value, it can be determined that the axial induction factor converges.

[0107] Correspondingly, when the axial induction factor has not converged, the iterative value of the axial induction factor is re-substituted into equations (4) to (7) to calculate new lift coefficients and drag coefficients, and the axial induction factor is recalculated through equation (8) using the new lift coefficients and drag coefficients until the axial induction factor converges and the aerodynamic load of the blade is calculated based on the converged axial induction factor.

[0108] At the same time, the aerodynamic load of the blade calculated above is the aerodynamic load of the blade element segment. By synthesizing the aerodynamic loads of the blade element segment, the aerodynamic load of the entire blade can be obtained. Next, the initial aerodynamic load of the blade is transmitted to the hub. The external load on the hub and its own gravity are superimposed to obtain the total aerodynamic load on the hub. This part of the aerodynamic load is then transmitted to the nacelle, and a new total aerodynamic load is obtained by superposition until it is transmitted to the tower to obtain the aerodynamic load on the tower, that is, subsequent calculations are performed based on the obtained tower aerodynamic load.

[0109] In S300, the tower is modeled, and the process of solving the elastic force of the tower based on the tower model is as follows: the tower is regarded as a flexible cantilever beam structure with continuous mass and stiffness distribution, and the modal analysis method is used to solve the modal equation of the structure after parameter modeling to obtain the modal vibration shape, and then the elastic force of the tower is obtained by linear fitting.

[0110] In S400, a dynamic equation of the tower based on inertia force, elastic force and aerodynamic load is constructed, and the inertia force of the tower is solved based on the dynamic equation of the tower. The specific steps are as follows:

[0111] The dynamic equation of the tower based on inertial force, elastic force and aerodynamic load is as follows:

[0112] F i T +F i W +F i G =0 (9)

[0113] In formula (9), F i T Represents elastic force, F i W is the external aerodynamic load, F i G is the inertial force.

[0114] In formula (9), the elastic force and the external aerodynamic load have been solved through the above steps. Correspondingly, based on the tower dynamics equation of formula (9), the inertia force of the tower can be obtained.

[0115] After the inertial force of the tower is calculated, since the output result in step S100 is a wind speed time series, correspondingly, the inertial force of the tower obtained is also a tower inertial force time series. Then, the inertial force time series at the top of the tower is extracted as the load time series of the yaw bearing (i.e., the boundary condition of the yaw bearing), and the variable amplitude load is identified in the time history. Peak-valley filtering is used to retain the peak-valley points and remove the intermediate transition points. Then, the raindrop flow counting method is used to decompose the complex history into multiple load cycles. Finally, the obtained series of load cycles are subjected to amplitude classification, cyclic statistics, etc., and finally a load time series is generated as the boundary condition of the yaw bearing (see Figure 6 ).

[0116] In S600, the construction process of the four-point contact ball bearing model is as follows:

[0117] Wind turbine yaw bearing is a four-point contact ball bearing, refer to Figure 7 The core of the modeling process is to calculate the interaction force between the rolling element and the raceway. The main difference between the modeling of four-point contact ball bearings and ordinary bearings is that there are four groove curvature centers on the inner and outer raceways, and the four contact effects need to be judged separately. The judgment of a single point is used as an example to illustrate the specific steps as follows:

[0118] In order to describe the motion of the rolling element and the ring, a set of coordinate systems is first constructed on the bearing, namely: i(x i ,y i ,z i ) inertial coordinate system; a(x a ,y a ,z a ) rolling element orientation coordinate system; r(x r ,y r ,zr ) Ring fixed body coordinate system, as shown in Figure X is the contact analysis diagram of the four-point contact ball bearing.

[0119] The position vectors of the center of mass of the ring and the rolling element relative to the inertial coordinate system are r r and r a , then the position vector of the ring mass center relative to the rolling element mass center can be expressed as:

[0120] r ra =r a -r r

[0121] The position vector between the center of mass of the ring and the center of curvature of the raceway groove is r rc , then the position vector of the rolling body mass center relative to the raceway groove curvature center can be expressed as:

[0122] r ca =r ra -r rc

[0123] The contact angle between the rolling element and the ring can be expressed as:

[0124]

[0125] Among them: subscripts 1 and 3 represent the first and third components of the vector respectively, that is, the components of the vector in different directions of the coordinate system.

[0126] The contact geometry approach can be expressed as:

[0127] δ=r ca13 -(f-0.5)D

[0128] Where: r ca13 It represents the projection of the vector on the plane formed by the first and third directions, f is the curvature coefficient of the raceway groove, and D is the rolling element diameter.

[0129] Using Hertz contact theory, the contact load can be given by:

[0130] Q=Kδ 1.5

[0131] Where: K is the Hertz contact coefficient.

[0132] Furthermore, the spalling damage model of the four-point contact ball bearing is as follows: Taking the spalling of the inner ring as an example, when the fault area appears at the contact point, the contact geometry approach needs to take into account the additional clearance. The geometric approach δ′ between the rolling element and the inner ring is the original geometric approach δ and the additional clearance δ d The difference between the inner ring and the inner ring is as follows: Figure 6 shown.

[0133]

[0134] Where: θ0 is the damage starting angle, θ1 is the damage width angle, θ g is the rolling element revolution angle.

[0135] After the contact load is calculated, the dynamic equation of the yaw bearing can be constructed through step S700, that is, the dynamic equation is constructed based on the boundary conditions and contact load of the yaw bearing, that is, a corresponding dynamic model in which the inner ring rotates and the outer ring is fixed is established, and only the radial translation of the inner ring in the x and y directions and the axial rotation in the z direction, the translation of the roller in the x, y and z directions and the rotation in the z direction, and the rotation of the retaining cage in the z direction are considered.

[0136] The translational motion of the inner ring and rolling elements can be described by Newton's laws:

[0137]

[0138] Where: m is the mass of the bearing element, is the acceleration vector, F is the external force on the component, and is the resultant force of the yaw bearing boundary condition and the contact load.

[0139] The rotation of the inner ring, rolling elements and cage in the z direction can be described by the Euler equation:

[0140]

[0141] Where: I is the moment of inertia of the bearing element; ω is the angular velocity; and M is the moment of force on the component, calculated from the resultant of the yaw bearing boundary conditions and the contact load.

[0142] Although embodiments of the present invention have been shown and described, it will be appreciated by those skilled in the art that various changes, modifications, substitutions and variations may be made to the embodiments without departing from the principles and spirit of the present invention, and that the scope of the present invention is defined by the appended claims and their equivalents.

Claims

1. A yaw bearing dynamics simulation method, characterized in that: include: Obtain the power spectrum density function of any position point in the wind field in different directions of the coordinate system, and calculate the wind speed time series of any position point in the wind field in different directions based on the power spectrum density function; Based on the wind speed time series, the blade aerodynamic load is calculated through the iterative calculation of blade element momentum theory, and the tower aerodynamic load is obtained based on the blade aerodynamic load; Model the tower and solve the elastic force of the tower based on the tower model; Construct the dynamic equation of the tower based on inertial force, elastic force and aerodynamic load, and solve the inertial force of the tower based on the dynamic equation of the tower; The yaw bearing boundary conditions are obtained based on the tower inertia force; Construct a four-point contact model of the yaw bearing and solve the contact load based on the four-point contact model of the yaw bearing; The dynamic equation of the yaw bearing is constructed based on the boundary conditions and contact loads of the yaw bearing.

2. A yaw bearing dynamics simulation method according to claim 1, characterized in that: The method of obtaining the power spectrum density function of any position point in the wind field in different directions in the coordinate system includes: Construct spatial correlation functions in different directions of the coordinate system, and calculate the integral scales in different directions of the coordinate system based on the spatial functions in different directions; The power spectral density function in different directions of the coordinate system is calculated based on the integral scale in different directions of the coordinate system.

3. A yaw bearing dynamics simulation method according to claim 2, characterized in that: The spatial correlation function in any direction on the coordinate system is as follows: Where i is the wind speed point to be calculated, j is k is the set of points that will produce wind speed correlation with the calculation point, r is the distance between two points in the selected direction, unit m, a k is the coherence intensity control coefficient between two points, is the set value, b k is the distance coherence intensity control coefficient between two points, is the set value, z m is the average height between the two points, f is the frequency in Hz, and U is the average wind speed between the two points in m / s.

4. A yaw bearing dynamics simulation method according to claim 3, characterized in that: The integral scale calculation formula in any direction in the coordinate system is as follows: r represents the distance between two points in any direction of the coordinate system, in meters.

5. A yaw bearing dynamics simulation method according to claim 4, characterized in that: The calculation formula of the power spectral density function in any direction of the coordinate system is as follows: Where S(f) represents the power spectral density in the selected direction, σ represents the standard deviation of wind speed in the selected direction, L represents the integral in the selected direction, U is the average wind speed in m / s, and f is the frequency in Hz.

6. A yaw bearing dynamics simulation method according to claim 1, characterized in that: The blade aerodynamic load obtained by iterative calculation of blade element momentum theory based on the wind speed time series includes: Given the initial value of the axial induction factor, the tangential induction factor is calculated; The blade lift and drag calculation formula 1 based on the lift coefficient and drag coefficient is constructed through the blade element theory; The second blade lift and drag calculation formula based on the axial induction factor and the tangential induction factor is constructed through momentum theory; Based on the equal relationship between blade lift and drag calculation formula 1 and calculation formula 2, the lift coefficient and drag coefficient are calculated by substituting the axial induction factor and the tangential induction factor; Construct a calculation formula for the axial induction factor based on the lift coefficient and the drag coefficient; Substitute the calculation formula of the axial induction factor based on the lift coefficient and the drag coefficient into the lift coefficient and the drag coefficient to obtain the iterative value of the axial induction factor; The convergence of the axial induction factor is judged. When the axial induction factor converges, the aerodynamic load of the blade is calculated using the current iterative value of the axial induction factor.

7. A yaw bearing dynamics simulation method according to claim 6, characterized in that: The calculation formula of the axial induction factor based on the lift coefficient and the drag coefficient is as follows: In the formula represents the inflow angle, C l is the lift coefficient, C d is the drag coefficient, a is the axial induction factor, σ′ is the solidity of the wind wheel, and the calculation formula is as follows: σ′=Bc / C Where C is the circumference of the wind wheel, c is the local chord length of the blade, and B is the number of blades.

8. A yaw bearing dynamics simulation method according to claim 1, characterized in that: The yaw bearing comprises a rolling element and a ring having an inner ring and an outer ring, and the contact load is solved based on the four-point contact model of the yaw bearing, comprising: The contact geometry approach between the rolling element and the inner ring is calculated, and the contact load is calculated based on the contact geometry approach.

9. A yaw bearing dynamics simulation method according to claim 8, characterized in that: The formula for calculating the contact load by contact geometry approach is as follows: q=Kδ 1.5 Where: K is the Hertz contact coefficient, Q is the contact load, and δ is the contact geometry approach.

10. A yaw bearing dynamics simulation system, characterized in that: include: An acquisition module is configured to acquire a power spectrum density function of any position point in the wind farm in different directions of the coordinate system, and calculate a wind speed time series of any position point in the wind farm in different directions based on the power spectrum density function; A calculation module is configured to obtain a blade aerodynamic load by iterative calculation based on a wind speed time series through blade element momentum theory, and obtain a tower aerodynamic load based on the blade aerodynamic load; A modeling module configured to model the tower and solve the tower elastic force based on the tower model; A first solving module is configured to construct a dynamic equation of the tower based on inertial force, elastic force and aerodynamic load, and solve the inertial force of the tower based on the dynamic equation of the tower; a determination module configured to obtain a yaw bearing boundary condition based on a tower inertia force; A second solving module is configured to construct a four-point contact model of the yaw bearing and solve the contact load based on the four-point contact model of the yaw bearing; A building module is configured to build the dynamic equations of the yaw bearing based on the yaw bearing boundary conditions and contact loads.

Images