Patent Publication Number: US-2022228568-A1

Title: Method for analyzing wind turbine blade coating fatigue due to rain erosion

Description:
CROSS-REFERENCE TO RELATED APPLICATIONS 
     The present application is a continuation of International Application No. PCT/CN2021/072812, filed on Jan. 20, 2021, which claims priority to Chinese Application No. 202110055626.2, filed on Jan. 15, 2021, the contents of both of which are incorporated herein by reference in their entireties. 
    
    
     TECHNICAL FIELD 
     The present application relates to the field of wind turbine blade design, in particular to a method for analyzing wind turbine blade coating fatigue due to rain erosion. 
     BACKGROUND 
     Wind turbine blades are frequently exposed to impacts from high-relative-speed objects such as rain, atmospheric particles and hail, especially at the tips of the blades. These impacts may lead to damage and peeling of the leading edge of the blade, thereby reducing the aerodynamic performance and power output of the wind turbine. With the continuous growth of the tip speed and rotor diameter of the wind turbine, it becomes more and more important to predict the fatigue life of the wind turbine blade coating due to rain erosion in the design stage. 
     At present, there is no effective solution in this aspect. The present application combines the stochastic rain field model, smoothed particle hydrodynamics and fatigue crack propagation theory to predict and calculate the fatigue life of the wind turbine blade coating. The existing impact approach and energy approach for calculating the impact stress of raindrops on the blade coating have some defects. It is difficult for the impact approach to consider the fluid-solid interaction in the process of raindrops impact, while it is difficult for the energy method to quantify the total transfer energy from the stochastic rain field to the blade coating of the wind turbine. At present, the calculation of the fatigue life of the wind turbine blade coating under rain erosion is usually implemented by using stress-life curve and Miner&#39;s hypothesis of linear accumulation, but the life calculated by this method is only limited to the fatigue crack initiation period. Usually, the fatigue failure of materials will go through three stages: crack initiation, crack stable propagation and crack unstable propagation. Traditional fatigue analysis and calculation methods cannot completely calculate the fatigue life of the wind turbine blade coating. 
     SUMMARY 
     The present application aims to overcome the shortcomings of the prior art, and provides a method for analyzing wind turbine blade coating fatigue due to rain erosion. This method combines a stochastic rain field model, smoothed particle hydrodynamics and the fatigue crack propagation theory to predict and calculate the fatigue life of a wind turbine blade coating. Through effective modeling of natural rainfall, accurate analysis of raindrop impact stress on the wind turbine blade and comprehensive calculation of the fatigue life of the wind turbine blade coating, the fatigue life of a wind turbine blade coating due to rain erosion is accurately and effectively calculated. 
     The purpose of the present application is realized by the following technical solution. 
     The present application relates to a method for analyzing wind turbine blade coating fatigue due to rain erosion, a stochastic rain field model is used to effectively model the natural rainfall condition, smoothed particle hydrodynamics and stress interpolation are used to accurately analyze the stress of raindrop impacting the blade, and the fatigue crack propagation theory is used to comprehensively calculate the fatigue life of the blade coating; the method specifically comprises the following steps: 
     S1: establishing a plurality of stochastic rain field models according to different rain intensities I and rainfall durations t s ; 
     S2: calculating stresses caused by different raindrops impacting the blades by using finite element analysis; 
     S3: calculating the impact stress of a coating in a stochastic rain field; 
     S4: calculating blade coating fatigue lives t I  under different rain intensities I; 
     S5: counting an annual rainfall duration t A  and a probability P I  of occurrence of each rain intensity; 
     S6: repeating steps S3 and S4 to obtain the fatigue life of the blade coating under different rain intensities I and calculating a wind turbine blade coating fatigue life t f  by using the following formula according to the calculation results of S4 and S5. 
     
       
         
           
             
               D 
               
                 1 
                 ⁢ 
                 
                     
                 
                 ⁢ 
                 year 
               
             
             = 
             
               
                 ∑ 
                 I 
               
               ⁢ 
               
                 
                   
                     P 
                     I 
                   
                   × 
                   
                     t 
                     A 
                   
                 
                 
                   t 
                   I 
                 
               
             
           
         
       
       
         
           
             
               t 
               f 
             
             = 
             
               
                 1 
                 
                   D 
                   
                     1 
                     ⁢ 
                     
                         
                     
                     ⁢ 
                     year 
                   
                 
               
               . 
             
           
         
       
     
     Furthermore, S1 specifically comprises the following steps: firstly, determining a number k of raindrops in a stochastic rain field, then determining parameters of each raindrop, including a diameter, a shape, an impact angle θ and an impact position of each raindrop, and constructing the stochastic rain field model according to relevant attributes of the k raindrops; 
     (1) the number k of the raindrops is calculated by the following formula: 
     
       
         
           
             
               P 
               ⁡ 
               
                 ( 
                 
                   
                     N 
                     ⁡ 
                     
                       ( 
                       V 
                       ) 
                     
                   
                   = 
                   k 
                 
                 ) 
               
             
             = 
             
               
                 
                   
                     ( 
                     
                       λ 
                       ⁢ 
                       
                           
                       
                       ⁢ 
                       V 
                     
                     ) 
                   
                   k 
                 
                 ⁢ 
                 
                   e 
                   
                     
                       - 
                       λ 
                     
                     ⁢ 
                     
                         
                     
                     ⁢ 
                     V 
                   
                 
               
               
                 k 
                 ! 
               
             
           
         
       
       
         
           
             λ 
             = 
             
               48.88 
               ⁢ 
               
                   
               
               ⁢ 
               
                 I 
                 0.15 
               
             
           
         
       
     
     where λ is an estimated number of raindrops per unit volume, P(N(V)=k) is a probability that there are k raindrops in a volume V, and I is a rain intensity in mm h −1 ; raindrops are uniformly distributed in a space of the volume V; 
     the rainfall space volume Vis calculated by the following formula: 
     
       
      
       V=S×ν×t 
       s  
      
     
     where S is a rainfall projection area, i.e., a blade coating area; ν is a relative velocity of raindrop impact, i.e., addition of a blade linear velocity with a raindrop velocity; t s  is a rainfall duration; 
     (2) the diameter of each raindrop is calculated by the following formula: 
     
       
         
           
             F 
             = 
             
               1 
               - 
               
                 exp 
                 ⁡ 
                 
                   [ 
                   
                     - 
                     
                       
                         ( 
                         
                           d 
                           
                             1.3 
                             ⁢ 
                             
                               I 
                               0.232 
                             
                           
                         
                         ) 
                       
                       2.25 
                     
                   
                   ] 
                 
               
             
           
         
       
     
     where F is a cumulative distribution function of a raindrop size d, d is the raindrop size in mm, and I is the rain intensity in mm h −1 ; 
     (3) determination of the shape of the raindrop is to determine a type of raindrop according to an occurrence probability of a type of raindrop, and to carry out geometric modeling according to specific types; 
     the raindrop has the shape of one of flat ellipsoid, spindle ellipsoid or perfect sphere, the occurrence probabilities of which are 27%, 55% and 18% respectively; for perfect-sphere raindrops, modeling is directly implemented according to a raindrop radius; for flat-ellipsoid and spindle-ellipsoid raindrops, geometric modeling of raindrops is completed by an axial ratio formula: 
         a= 1.030−0.124r 0  
 
     where α is the axial ratio of a minor axis to a major axis, r 0  is the equivalent perfect-sphere raindrop radius, i.e., r 0 =d/2 ; 
     (4) the impact angle 0 of the raindrop follows a uniform distribution of [0, 90°]; 
     (5) the impact position of the raindrop is any position in the blade coating area and is evenly distributed. 
     Furthermore, S2 specifically comprises the following substeps: 
     S2.1: constructing a blade model, meshing, setting properties of related composite materials, and setting constraint conditions: 
     S2.2: according to different sizes and shapes of raindrops, construing different single raindrops, meshing, setting an impact speed and impact angle of the raindrops, using finite element analysis software (e.g., Abaqus) combined with a smooth fluid dynamics method to implement simulation analysis, and calculating the impact stress of a single raindrop; 
     S2.3: obtaining Von Mises stresses of various sites on the blade coating in the finite element analysis as the impact stresses; as an embodiment, MATLAB can be used for obtaining the stress; 
     S2.4: repeating steps S2.2-S2.3 to calculate the impact stress of the raindrops under various conditions, including a combination of different raindrop diameters, different raindrop shapes, different impact angles and different impact speeds, for example, 9 types of raindrop diameters (d=1, 2, 3, 4, 5, 6, 7, 8, 9 mm), 3 types of raindrop shapes (flat ellipsoid, spindle ellipsoid, perfect sphere) and 6 types of impact angles (θ=15°, 30°, 45°, 60°,75°, 90°) and one impact speed (v=90 ms −1 ) 
     Furthermore, S3 specifically includes the following substeps: 
     S3.1: according to the rain field model constructed by S1, after determining the size, shape, impact angle and speed of a single random raindrop, a circular domain with the impact point as a center and N times of the raindrop diameter as the radius being considered as the area influenced by raindrop impact, and N is 9-11; 
     S3.2: choosing the same type of raindrop shape according to the stress of the raindrop impact in a series of cases obtained in step S2, and searching for stress results of the impact cases that have the closest raindrop diameter, impact angle, and impact speed to interpolate the stress in the circular area; 
     S3.3: repeat steps S3.1-S3.2 for each raindrop until all the impact stresses caused by k raindrops on the blade are calculated. 
     Furthermore, S4 specifically includes the following substeps: 
     S4.1: selecting the rain intensity I and the rainfall duration t s  of a single simulation, and calculating the impact stress of the coating in the stochastic rain field according to steps S1 to S3; 
     S4.2: selecting a local maximum stress and a neighboring minimum stress, or selecting a local minimum stress and a neighboring maximum stress to constitute a half stress cycle, and splitting an impact stress curve into a plurality of half-cyclic stresses with constant amplitudes; 
     S4.3: for each half-cyclic stress in S4.2, calculating the number of allowable stress cycles N f  by using the following formula: 
     
       
         
           
             
               σ 
               a 
               ′ 
             
             = 
             
               
                 
                   σ 
                   a 
                 
                 ⁢ 
                 UTS 
               
               
                 UTS 
                 - 
                 
                   σ 
                   m 
                 
               
             
           
         
       
       
         
           
             
               N 
               f 
             
             = 
             
               
                 ( 
                 
                   
                     σ 
                     a 
                     ′ 
                   
                   
                     σ 
                     f 
                   
                 
                 ) 
               
               
                 1 
                 / 
                 b 
               
             
           
         
       
     
     where σ′ a  is the corrected stress amplitude, σ a , is the stress amplitude, σ m  is the mean stress, UTS is the ultimate tensile strength, σ f  is the fatigue strength coefficient, b is the fatigue strength exponent, UTS, σ f  and b are all inherent properties of a coating material, which can be obtained through experiments, while σ a , and σ m  can be calculated according to the maximum stress and minimum stress of the half-cyclic stress; 
     S4.4: repeating step S4.3 until the number of the allowable stress cycles N f  of all half-cyclic stresses is calculated; according to Miner&#39;s rule for damage accumulation, a fatigue damage caused by all half-cyclic stresses caused by a raindrop impacting the blade is 
     
       
         
           
             D 
             = 
             
               
                 ∑ 
                 i 
               
               ⁢ 
               
                 0.5 
                 
                   N 
                   f 
                   i 
                 
               
             
           
         
       
     
     S4.5: repeating steps S4.2-S4.4 until the fatigue damage D s  caused by the impact stress of k raindrops on the blade in the rainfall duration t s  is calculated, and calculating the fatigue life t initiation  of a crack initiation period by the following formula: 
     
       
         
           
             
               t 
               initiation 
             
             = 
             
               
                 t 
                 s 
               
               
                 D 
                 s 
               
             
           
         
       
     
     S4.6: for each half-cyclic stress in S4.2, iteratively calculating a crack length by the following formula: 
         a   i+1   =a   i +0.5 ×C [ Y (σ max −σ min ) √{square root over (π a   i )}] m  
 
     where a i+1  is a crack length after the half-cyclic stress, a i  is a crack length before the half-cyclic stress; C and m are inherent properties of the material, which are obtained through material fatigue experiments; a value of Y is determined by a crack shape, σ max  is the maximum stress of the half-cyclic stress and σ min  is the minimum stress of the half-cyclic stress; 
     S4.7: repeating steps S4.2 and S4.6 until the crack length a caused by the impact stress of k raindrops on the blade in the rainfall duration t s  is calculated; 
     S4.8: if the rain intensity I is greater than or equal to 10 mm h −1 , proceeding to step S4.9; if the rain intensity I is less than 10 mm h −1 , proceeding to step S4.1; 
     S4.9: repeating steps S4.1, S4.2, S4.6, S4.7, and the rainfall duration increases continuously, while the crack length increases continuously until the crack length meets the following formula or the crack length is greater than a coating thickness, considering that a crack stable propagation period is completed: 
         Yσ   max √{square root over (π a   now )}&gt; K   C  
 
     where a now  is a current crack length, K C  is a fracture toughness, which is an inherent property of the material and can be measured by experiments; when the crack length meets the above conditions, the rainfall duration is the fatigue life during the crack stable propagation period; 
     S4.10: when the rain intensity I is low, calculating an equivalent stress range Δσ within the rainfall duration t s  by the following formula, and using a constant amplitude cyclic stress of the equivalent stress range Δσ to replace all varied-amplitude cyclic stresses within the rainfall duration I′ s : 
     
       
         
           
             
               Δ 
               ⁢ 
               
                   
               
               ⁢ 
               σ 
             
             = 
             
               { 
               
                 
                   
                     
                       
                         
                           { 
                           
                             
                               2 
                               
                                 
                                   
                                     N 
                                     t 
                                   
                                   ⁡ 
                                   
                                     ( 
                                     
                                       m 
                                       - 
                                       2 
                                     
                                     ) 
                                   
                                 
                                 ⁢ 
                                 
                                   
                                     C 
                                     ⁡ 
                                     
                                       ( 
                                       
                                         Y 
                                         ⁢ 
                                         
                                           π 
                                         
                                       
                                       ) 
                                     
                                   
                                   m 
                                 
                               
                             
                             ⁡ 
                             
                               [ 
                               
                                 
                                   a 
                                   0 
                                   
                                     ( 
                                     
                                       1 
                                       - 
                                       
                                         m 
                                         2 
                                       
                                     
                                     ) 
                                   
                                 
                                 - 
                                 
                                   a 
                                   
                                     ( 
                                     
                                       1 
                                       - 
                                       
                                         m 
                                         2 
                                       
                                     
                                     ) 
                                   
                                 
                               
                               ] 
                             
                           
                           } 
                         
                         
                           1 
                           m 
                         
                       
                       , 
                     
                   
                   
                     
                       m 
                       ≠ 
                       2 
                     
                   
                 
                 
                   
                     
                       
                         
                           [ 
                           
                             
                               1 
                               
                                 
                                   
                                     CN 
                                     t 
                                   
                                   ⁡ 
                                   
                                     ( 
                                     
                                       Y 
                                       ⁢ 
                                       
                                         π 
                                       
                                     
                                     ) 
                                   
                                 
                                 m 
                               
                             
                             ⁢ 
                             
                               ln 
                               ⁡ 
                               
                                 ( 
                                 
                                   a 
                                   
                                     a 
                                     0 
                                   
                                 
                                 ) 
                               
                             
                           
                           ] 
                         
                         
                           1 
                           m 
                         
                       
                       , 
                     
                   
                   
                     
                       m 
                       = 
                       2 
                     
                   
                 
               
             
           
         
       
     
     where a 0  is an initial crack length, a is a crack length after the rainfall duration t s  and N t  is the number of all stress cycles in the rainfall duration t s ; 
     the number of allowable stress cycles N c  during the crack stable propagation period is calculated by the following formula: 
     
       
         
           
             
               N 
               c 
             
             = 
             
               { 
               
                 
                   
                     
                       
                         
                           
                             
                               2 
                               
                                 
                                   ( 
                                   
                                     m 
                                     - 
                                     2 
                                   
                                   ) 
                                 
                                 ⁢ 
                                 
                                   
                                     C 
                                     ⁡ 
                                     
                                       ( 
                                       
                                         Y 
                                         ⁢ 
                                         
                                             
                                         
                                         ⁢ 
                                         Δ 
                                         ⁢ 
                                         
                                             
                                         
                                         ⁢ 
                                         σ 
                                         ⁢ 
                                         
                                           π 
                                         
                                       
                                       ) 
                                     
                                   
                                   m 
                                 
                               
                             
                             ⁡ 
                             
                               [ 
                               
                                 
                                   a 
                                   0 
                                   
                                     ( 
                                     
                                       1 
                                       - 
                                       
                                         m 
                                         2 
                                       
                                     
                                     ) 
                                   
                                 
                                 - 
                                 
                                   a 
                                   c 
                                   
                                     ( 
                                     
                                       1 
                                       - 
                                       
                                         m 
                                         2 
                                       
                                     
                                     ) 
                                   
                                 
                               
                               ] 
                             
                           
                           , 
                         
                       
                       
                         
                           m 
                           ≠ 
                           2 
                         
                       
                     
                     
                       
                         
                           
                             
                               1 
                               
                                 
                                   C 
                                   ⁡ 
                                   
                                     ( 
                                     
                                       Y 
                                       ⁢ 
                                       
                                           
                                       
                                       ⁢ 
                                       Δ 
                                       ⁢ 
                                       
                                           
                                       
                                       ⁢ 
                                       σ 
                                       ⁢ 
                                       
                                         π 
                                       
                                     
                                     ) 
                                   
                                 
                                 m 
                               
                             
                             ⁢ 
                             
                               ln 
                               ⁡ 
                               
                                 ( 
                                 
                                   
                                     a 
                                     c 
                                   
                                   
                                     a 
                                     0 
                                   
                                 
                                 ) 
                               
                             
                           
                           , 
                         
                       
                       
                         
                           m 
                           = 
                           2 
                         
                       
                     
                   
                   ⁢ 
                   
                     
 
                   
                   ⁢ 
                   
                     a 
                     C 
                   
                 
                 = 
                 
                   
                     
                       ( 
                       
                         
                           K 
                           C 
                         
                         
                           Y 
                           ⁢ 
                           
                               
                           
                           ⁢ 
                           
                             σ 
                             MAX 
                           
                         
                       
                       ) 
                     
                     2 
                   
                   / 
                   π 
                 
               
             
           
         
       
     
     where σ MAX  is the maximum stress in the rainfall duration t s ; 
     the fatigue life of the crack stable propagation period is calculated by the following formula: 
     
       
         
           
             
               t 
               propagation 
             
             = 
             
               
                 
                   N 
                   c 
                 
                 
                   N 
                   t 
                 
               
               ⁢ 
               
                 t 
                 s 
               
             
           
         
       
     
     S4.11: calculating the fatigue life of the coating at a certain point under the rain intensity I by the following formula: 
     
       
      
       t 
       IP 
       =t 
       initiation 
       +t 
       propagation  
      
     
     S4.12: repeating steps S4.1-S4.12 to calculate the fatigue life of each point of the coating, ranking the fatigue lives of all points from small to large, and taking the fatigue life of the 84th percentile as the fatigue life t I  of the whole coating. 
     Furthermore, S5 specifically includes the following substeps: 
     S5.1: obtaining annual rainfall data of an area where the wind turbine is located according to relevant statistical data; 
     S5.2: statistically processing the rainfall data, and obtaining an annual rainfall duration t A  and a probability of occurrence of each rain intensity P I  in the area (i.e., a probability density function PDF or a probability mass function PMF). 
     The present application has the following beneficial effects: 
     (1) the stochastic rain field model proposed by the present application takes the raindrop shapes (perfect sphere, flat ellipsoid and spindle ellipsoid) and the real raindrop sizes into consideration, and the stochastic rain field model well reflects a real rain field situation. 
     (2) The present application uses the smooth particle hydrodynamics (SPH) and stress interpolation method to calculate the impact stress of raindrops in the process of random rainfall, and the method can effectively and accurately calculate the impact stress of raindrops on the coating, while ensuring that the calculation time is not too long. 
     (3) According to the fatigue crack propagation theory, the fatigue life of the coating in the crack initiation period and the fatigue life in the crack stable propagation period are completely calculated, so that the calculated fatigue life is more accurate. 
    
    
     
       BRIEF DESCRIPTION OF DRAWINGS 
         FIG. 1  is a flow chart of the method of the present application; 
         FIG. 2  is a schematic diagram of the method of the present application; 
         FIG. 3  is a schematic diagram of the raindrop shape and impact angle; 
         FIG. 4  is a simulation diagram of a stochastic rain field under four rain intensities, (a) 1 mm h −1 , (b) 10 mm h −1 , (c) 20 mm h −1 , and (d) 50 mm h −1 ; 
         FIG. 5  is a model view of a tip panel of a blade; 
         FIG. 6  is a stress cloud of a single raindrop impacting the blade at eight intervals (0 μs, 10 μs, 20 μs, 30 μs, 40 μs, 50 μs); 
         FIG. 7  is an interpolation calculation result diagram of an impact stress of a raindrop with a diameter of 2.5 mm and an impact angle of 80°, in which (a) is a comparison diagram of an stress interpolation calculation result and raindrop impact stresses under four closest impact cases, and (b) is a comparison diagram of a stress interpolation calculation result and a finite element analysis calculation result; 
         FIG. 8  is a probability mass function diagram of rain intensity in Miami, Florida. 
     
    
    
     DESCRIPTION OF EMBODIMENTS 
     The purpose and effect of the present application will become more clear from the following detailed description of the present application according to the drawings and preferred embodiments. It should be understood that the specific embodiments described here are only used to explain, not to limit, the present application. 
     According to the method for analyzing wind turbine blade coating fatigue due to rain erosion, a stochastic rain field model is used to effectively model the natural rainfall condition, smoothed particle hydrodynamics and stress interpolation are used to accurately analyze the stress of raindrops impacting the blades, and the fatigue crack propagation theory is used to comprehensively calculate the fatigue life of the blade coating. The fatigue life of a blade coating of the wind turbine located in Miami, Florida is predicted and calculated. The specific flow chart is shown in  FIG. 1  and the schematic diagram is shown in  FIG. 2 , which specifically includes the following steps: 
     S1, a plurality of stochastic rain field models are established according to different rain intensities I and rainfall durations t s ; 
     S1.1, the number k of the raindrops is calculated by the following formula: 
     
       
         
           
             
               
                 
                   
                     
                       P 
                       ⁡ 
                       
                         ( 
                         
                           
                             N 
                             ⁡ 
                             
                               ( 
                               V 
                               ) 
                             
                           
                           = 
                           k 
                         
                         ) 
                       
                     
                     = 
                     
                       
                         
                           
                             ( 
                             
                               λ 
                               ⁢ 
                               
                                   
                               
                               ⁢ 
                               V 
                             
                             ) 
                           
                           k 
                         
                         ⁢ 
                         
                           e 
                           
                             
                               - 
                               λ 
                             
                             ⁢ 
                             
                                 
                             
                             ⁢ 
                             V 
                           
                         
                       
                       
                         k 
                         ! 
                       
                     
                   
                   ⁢ 
                   
                     
 
                   
                   ⁢ 
                   
                     λ 
                     = 
                     
                       4 
                       ⁢ 
                       
                         8 
                         . 
                         8 
                       
                       ⁢ 
                       8 
                       ⁢ 
                       
                           
                       
                       ⁢ 
                       
                         I 
                         0.15 
                       
                     
                   
                 
               
               
                 
                   ( 
                   V 
                   ) 
                 
               
             
           
         
       
     
     where λ is an estimated number of raindrops per unit volume, P(N(V)=k) is a probability that there are k raindrops in a volume V, and I is a rain intensity (mm h −1 ); raindrops are uniformly distributed in a space of the volume V; 
     the rainfall space volume Vis calculated by the following formula: 
     
       
      
       V=S×ν×t 
       s  
      
     
     where S is a rainfall projection are (i.e., a blade coating area), ν is a relative velocity of raindrop impact (addition of a blade linear velocity with a raindrop velocity), and t s  is a rainfall duration; a random number conforming to the above probability distribution is generated in MAILAB to obtain the raindrop number k; 
     S1.2, the diameter of each raindrop is calculated by the following formula: 
     
       
         
           
             F 
             = 
             
               1 
               - 
               
                 exp 
                 ⁡ 
                 
                   [ 
                   
                     - 
                     
                       
                         ( 
                         
                           d 
                           
                             1.3 
                             ⁢ 
                             
                               I 
                               0.232 
                             
                           
                         
                         ) 
                       
                       2.25 
                     
                   
                   ] 
                 
               
             
           
         
       
     
     where F is a cumulative distribution function of a raindrop size d, d is the raindrop size (mm), and I is the rain intensity (mm h −1 ); raindrops are uniformly distributed in a space of volume V, and the raindrop size d is obtained by generating a random number conforming to the above probability distribution in MATLAB; 
     S1.3, the shape of the raindrop includes perfect sphere, flat ellipsoid and spindle ellipsoid; for ellipsoidal raindrops, there is a minor axis and a major axis, an axial ratio of which is α which can be calculated by the following formula: 
         a= 1.030−0.124 r   0  
 
     where r 0  is an equivalent perfect sphere raindrop radius, i.e., r 0 =d/2 ; 
     Flat-ellipsoid raindrops have the longest axis in the horizontal plane, while spindle-ellipsoid raindrops have the longest axis perpendicular to the horizontal plane. The horizontal cross-sectional area of flat-ellipsoid raindrops and spindle-ellipsoid raindrops is assumed to be a circle and the vertical cross-sectional area is an ellipse, so the geometric modeling of raindrops can be completed by the axis ratio formula. According to relevant data, the probabilities of occurrence of flat-ellipsoid raindrops, spindle-ellipsoid raindrops and perfect-sphere raindrops are 27%, 55% and 18% respectively, as shown in  FIG. 3 ; 
     S 1 . 4 : the impact angle θ of the raindrop follows a uniform distribution of [0, 90°]; the impact positions of raindrops are any positions in the blade coating area, which are evenly distributed, as shown in  FIG. 3 ; 
     S 1 . 5 : for each raindrop, steps S1.2-S1.4 are repeated to determine the related attributes of each raindrop until the related attributes of k raindrops are determined, as shown in  FIG. 4 ; 
     S 2 : stresses caused by different raindrops impacting the blades are calculated and analyzed by using finite element analysis; 
     S2.1: a blade model is constructed, and meshing is carried out; in order to control the calculation amount, only a finite element model is built for one panel at the tip of the blade, as shown in  FIG. 5 , and properties of a composite material are set; as shown in Table 1 below, epoxy resin is used as the coating, and the bottom and side of the panel are set to be completely constrained: 
     
       
         
           
               
             
               
                 TABLE 1 
               
             
            
               
                   
               
               
                 Attribute table of the composite material of the blade 
               
            
           
           
               
               
            
               
                   
                 Material type 
               
            
           
           
               
               
               
               
            
               
                 Material attribute 
                 Coating 
                 QQ1 
                 Foam 
               
               
                   
               
            
           
           
               
               
               
               
            
               
                 Longitudinal Young’s modulus E 1  (GPa) 
                 3.44 
                 33.1 
                 0.256 
               
               
                 Transverse Young’s modulus E 2  (GPa) 
                 3.44 
                 17.1 
                 0.256 
               
               
                 Poisson ratio’s v 12   
                 0.3 
                 0.27 
                 0.33 
               
               
                 Shear modulus G 12  (GPa) 
                 1.38 
                 6.29 
                 0.098 
               
               
                 Density ρ (kg/m 3 ) 
                 1235 
                 1919 
                 200 
               
               
                   
               
            
           
         
       
     
     S2.2: according to different sizes and shapes of raindrops, different single raindrops are constructed, meshing is carried out, an impact speed and impact angle of the raindrops are set, a smoothed particle hydrodynamics (SP1-1) method in Abaqus finite element analysis software is used to calculate the impact stress of a single raindrop, as shown in  FIG. 6 ; 
     S 2 . 3 : Von Mises stresses of various sites on the blade coating in the finite element analysis is obtained using MATLAB as the impact stresses; 
     S2.4: steps S2.2-S2.3 are repeated, and the raindrop impact stresses under 162 cases are simulated and calculated, namely, nine raindrop diameters (d=1, 2, 3, 4, 5, 6, 7, 8, and 9 mm), three raindrop shapes (flat ellipsoid, spindle ellipsoid, and perfect sphere), six impact angles (θ=15°, 30°, 45°, 60°,75°, 90°)and one impact speed (v =90 ms  1 ); 
     S3: an impact stress of a coating in a stochastic rain field is calculated; 
     S3.1: according to the rain field model constructed by S1, after determining the size, shape, impact angle and speed of a single random raindrop, a circular domain with an impact point as a center and N times of the raindrop diameter as the radius is considered as the area influenced by raindrop impact; 
     S3.2: the same type of raindrop shape is chosen according to the stress of the raindrop impact in a series of cases obtained in step S2, and stress results of the impact cases calculated in S2 that have the closest raindrop diameter, impact angle, and impact speed are searched to interpolate the stress in the circular area; 
     S3.3: steps S3.1-S3.2 are repeated for each raindrop until all the impact stresses caused by k raindrops on the blade are calculated; 
     S4: blade coating fatigue lives t 1  under different rain intensities I are calculated; 
     S4.1: the rain intensity I and the rainfall duration t s  (eg., 10 min) of a single simulation are selected, and the impact stress of the coating in the stochastic rain field is calculated according to steps S1 to S3; 
     S4.2: The actual impact stress has varied stress amplitudes due to the randomness of raindrop impact ; for cycle-by-cycle fatigue analysis, a simple range counting method is used to count all half-cyclic stresses, i.e., selecting a local maximum (minimum) stress and a neighboring minimum (maximum) stress to constitute a half stress cycle, by which, a complex impact stress curve is split into a plurality of half-cyclic stresses with constant amplitudes; 
     S4.3: for each half-cyclic stress in S4.2, the number of allowable stress cycles N f  is calculated by using the following formula: 
     
       
         
           
             
               σ 
               a 
               ′ 
             
             = 
             
               
                 
                   σ 
                   a 
                 
                 ⁢ 
                 UTS 
               
               
                 UTS 
                 - 
                 
                   σ 
                   m 
                 
               
             
           
         
       
       
         
           
             
               N 
               f 
             
             = 
             
               
                 ( 
                 
                   
                     σ 
                     a 
                     ′ 
                   
                   
                     σ 
                     f 
                   
                 
                 ) 
               
               
                 1 
                 / 
                 b 
               
             
           
         
       
     
     where σ′ a  is the corrected stress amplitude, σ a  is the stress amplitude, σ m  is the mean stress, UTS is the ultimate tensile strength, σ f  is the fatigue strength coefficient, b is the fatigue strength exponent, UTS=73.3 MPa, σ f =83.3 MPa, b=−0.117. σ a  and σ m  can be calculated according to the maximum stress and minimum stress of the half-cyclic stress; 
     S4.4: step S4.3 are repeated until the number of the allowable stress cycles N f  of all half-cyclic stresses is calculated; according to Miner&#39;s rule for damage accumulation, a fatigue damage caused by all half-cyclic stresses caused by a raindrop impacting the blade is 
     
       
         
           
             D 
             = 
             
               
                 ∑ 
                 i 
               
               ⁢ 
               
                 
                   0 
                   . 
                   5 
                 
                 
                   N 
                   f 
                   i 
                 
               
             
           
         
       
     
     S4.5: steps S4.2-S4.4 are repeated until the fatigue damage D s  caused by the impact stress of k raindrops on the blade in the rainfall duration t s  is calculated, and the fatigue life t dinitiation  of a crack initiation period is calculated by the following formula: 
     
       
         
           
             
               t 
               initiation 
             
             = 
             
               
                 t 
                 s 
               
               
                 D 
                 s 
               
             
           
         
       
     
     S4.6: for each half-cyclic stress in S4.2, a crack length is iteratively calculated by the following formula: 
         a   i+1   =a   i +0.5 ×C [ Y (σ max −σ min )√{square root over (π a   i )}] m  
 
     where a i+1  is a crack length after the half-cyclic stress, a i  is a crack length before the half-cyclic stress; C=9.7, m=0.08; the value of Y is determined by a crack shape, and Y=1 in this embodiment; σ a  is the maximum stress of the half-cyclic stress and σ min  is the minimum stress of the half-cyclic stress; 
     S4.7: steps S4.2 and S4.6 are repeated until the crack length a caused by the impact stress of k raindrops on the blade in the rainfall duration t s  is calculated; 
     S4.8: if the rain intensity I is greater than or equal to 10 mm h −1 , proceeding to step S4.9; if the rain intensity I is less than 10 mm h −1 , proceeding to step S4.1; 
     S4.9: steps S4.1, S4.2, S4.6, S4.7 are repeated, and the rainfall duration increases continuously, while the crack length increases continuously until the crack length meets the following formula or the crack length is greater than a coating thickness, it is considered that a crack stable propagation period is completed: 
         Y σ max √{square root over (π now )}&gt; K   C  
 
     where, a now  is a current crack length, K C  is a fracture toughness, which is an inherent property of the material, K c =0.59 MPa h 1/2  in this embodiment; when the crack length meets the above conditions, the rainfall duration is the fatigue life during the crack stable propagation period; 
     S4.10: when the rain intensity I is low, the method of S4.9 needs a lot of iterative calculation, which takes a long time, so the method of S4.10 is proposed; an equivalent stress range Δσ within the rainfall duration t s  is calculated by the following formula, and a constant amplitude cyclic stress of the equivalent stress range Δσ is used to replace all varied-amplitude cyclic stresses within the rainfall duration t s : 
     
       
         
           
             
               Δ 
               ⁢ 
               σ 
             
             = 
             
               { 
               
                 
                   
                     
                       
                         
                           { 
                           
                             
                               2 
                               
                                 
                                   
                                     N 
                                     t 
                                   
                                   ⁡ 
                                   
                                     ( 
                                     
                                       m 
                                       - 
                                       2 
                                     
                                     ) 
                                   
                                 
                                 ⁢ 
                                 
                                   
                                     C 
                                     ⁡ 
                                     
                                       ( 
                                       
                                         Y 
                                         ⁢ 
                                         
                                           π 
                                         
                                       
                                       ) 
                                     
                                   
                                   m 
                                 
                               
                             
                             ⁡ 
                             
                               [ 
                               
                                 
                                   a 
                                   0 
                                   
                                     ( 
                                     
                                       1 
                                       - 
                                       
                                         m 
                                         2 
                                       
                                     
                                     ) 
                                   
                                 
                                 - 
                                 
                                   a 
                                   
                                     ( 
                                     
                                       1 
                                       - 
                                       
                                         m 
                                         2 
                                       
                                     
                                     ) 
                                   
                                 
                               
                               ] 
                             
                           
                           } 
                         
                         
                           1 
                           m 
                         
                       
                       , 
                     
                   
                   
                     
                       m 
                       ≠ 
                       2 
                     
                   
                 
                 
                   
                     
                       
                         
                           [ 
                           
                             
                               1 
                               
                                 C 
                                 ⁢ 
                                 
                                   
                                     
                                       N 
                                       t 
                                     
                                     ⁡ 
                                     
                                       ( 
                                       
                                         Y 
                                         ⁢ 
                                         
                                           π 
                                         
                                       
                                       ) 
                                     
                                   
                                   m 
                                 
                               
                             
                             ⁢ 
                             
                               ln 
                               ⁡ 
                               
                                 ( 
                                 
                                   a 
                                   
                                     a 
                                     0 
                                   
                                 
                                 ) 
                               
                             
                           
                           ] 
                         
                         
                           1 
                           m 
                         
                       
                       , 
                     
                   
                   
                     
                       m 
                       = 
                       2 
                     
                   
                 
               
             
           
         
       
     
     where a 0  is an initial crack length, a 0 =12 μm, a is a crack length after the rainfall duration t s  and N t  is the number of all stress cycles in the rainfall duration t s ; 
     the number of allowable stress cycles N c  during the crack stable propagation period is calculated by the following formula: 
     
       
         
           
             
               N 
               c 
             
             = 
             
               { 
               
                 
                   
                     
                       
                         
                           
                             
                               2 
                               
                                 
                                   ( 
                                   
                                     m 
                                     - 
                                     2 
                                   
                                   ) 
                                 
                                 ⁢ 
                                 
                                   
                                     C 
                                     ⁡ 
                                     
                                       ( 
                                       
                                         Y 
                                         ⁢ 
                                         
                                             
                                         
                                         ⁢ 
                                         Δ 
                                         ⁢ 
                                         
                                             
                                         
                                         ⁢ 
                                         σ 
                                         ⁢ 
                                         
                                           π 
                                         
                                       
                                       ) 
                                     
                                   
                                   m 
                                 
                               
                             
                             ⁡ 
                             
                               [ 
                               
                                 
                                   a 
                                   0 
                                   
                                     ( 
                                     
                                       1 
                                       - 
                                       
                                         m 
                                         2 
                                       
                                     
                                     ) 
                                   
                                 
                                 - 
                                 
                                   a 
                                   c 
                                   
                                     ( 
                                     
                                       1 
                                       - 
                                       
                                         m 
                                         2 
                                       
                                     
                                     ) 
                                   
                                 
                               
                               ] 
                             
                           
                           , 
                         
                       
                       
                         
                           m 
                           ≠ 
                           2 
                         
                       
                     
                     
                       
                         
                           
                             
                               1 
                               
                                 
                                   C 
                                   ⁡ 
                                   
                                     ( 
                                     
                                       Y 
                                       ⁢ 
                                       
                                           
                                       
                                       ⁢ 
                                       Δ 
                                       ⁢ 
                                       
                                           
                                       
                                       ⁢ 
                                       σ 
                                       ⁢ 
                                       
                                         π 
                                       
                                     
                                     ) 
                                   
                                 
                                 m 
                               
                             
                             ⁢ 
                             
                               ln 
                               ⁡ 
                               
                                 ( 
                                 
                                   
                                     a 
                                     c 
                                   
                                   
                                     a 
                                     0 
                                   
                                 
                                 ) 
                               
                             
                           
                           , 
                         
                       
                       
                         
                           m 
                           = 
                           2 
                         
                       
                     
                   
                   ⁢ 
                   
                     
 
                   
                   ⁢ 
                   
                     a 
                     C 
                   
                 
                 = 
                 
                   
                     
                       ( 
                       
                         
                           K 
                           C 
                         
                         
                           Y 
                           ⁢ 
                           
                               
                           
                           ⁢ 
                           
                             σ 
                             MAX 
                           
                         
                       
                       ) 
                     
                     2 
                   
                   / 
                   π 
                 
               
             
           
         
       
     
     where σ MAX  is the maximum stress in the rainfall duration t s ; 
     the fatigue life of the crack stable propagation period is calculated by the following formula: 
     
       
         
           
             
               t 
               propagation 
             
             = 
             
               
                 
                   N 
                   c 
                 
                 
                   N 
                   t 
                 
               
               ⁢ 
               
                 t 
                 s 
               
             
           
         
       
     
     S4.11: the fatigue life of the coating at a certain point under the rain intensity I is calculated by the following formula: 
     
       
      
       t 
       IP 
       =t 
       initiation 
       +t 
       propagation  
      
     
     S4.12: steps S4.1-S4.12 are repeated to calculate the fatigue life of each point of the coating, the fatigue lives of all points are ranked from small to large, and the fatigue life of the 84th percentile is taken as the fatigue life t 1  of the whole coating; 
     S5: an annual rainfall duration t A  and a probability P I  of occurrence of each rain intensity are counted; 
     S5.1: annual rainfall data of an area where the wind turbine is located are obtained according to relevant statistical data; 
     S5.2: the rainfall data are statistically processed to obtain an annual rainfall duration t A  and a probability of occurrence of each rain intensity P I  in the area (i.e., the probability density function PDF or probability mass function PMF, as shown in  FIG. 8 .) 
     S6: steps S3 and S4 are repeated to obtain the fatigue life of the blade coating under different rain intensities; 
     
       
         
           
               
             
               
                 TABLE 2 
               
             
            
               
                   
               
               
                 Fatigue life of a wind turbine blade coating under various rain intensities 
               
            
           
           
               
               
               
               
            
               
                 Rain intensity 
                 Fatigue life 
                 Rain intensity 
                 Fatigue life  
               
               
                 (mm h −1 )) 
                 (h) 
                 (mm h −1 ) 
                 (h) 
               
               
                   
               
            
           
           
               
               
               
               
            
               
                 20 
                 4.2 
                 10 
                 192.7 
               
               
                 19 
                 6.9 
                 9 
                 470.4 
               
               
                 18 
                 8.3 
                 8 
                 1254.5 
               
               
                 17 
                 14 
                 7 
                 1989.2 
               
               
                 16 
                 15.5 
                 6 
                 4155.7 
               
               
                 15 
                 31.3 
                 5 
                 14463 
               
               
                 14 
                 45.4 
                 4 
                 53673.3 
               
               
                 13 
                 46.4 
                 3 
                 200250 
               
               
                 12 
                 79 
                 2 
                 1590481.9 
               
               
                 11 
                 142.5 
                 1 
                 44960142.3 
               
               
                   
               
            
           
         
       
     
     According to the statistical results of S5, combined with the fatigue life of the wind turbine blade coating under various rain intensities in Table 2, the following formula is used to calculate the fatigue life t f  of the wind turbine blade coating. 
     
       
         
           
             
               
                 D 
                 
                   1 
                   ⁢ 
                   
                       
                   
                   ⁢ 
                   year 
                 
               
               = 
               
                 
                   ∑ 
                   I 
                 
                 ⁢ 
                 
                   
                     
                       P 
                       I 
                     
                     × 
                     
                       t 
                       A 
                     
                   
                   
                     t 
                     I 
                   
                 
               
             
             ⁢ 
             
               
 
             
             ⁢ 
             
               
                 t 
                 f 
               
               = 
               
                 1 
                 
                   D 
                   
                     1 
                     ⁢ 
                     
                         
                     
                     ⁢ 
                     year 
                   
                 
               
             
           
         
       
     
     The fatigue life of the wind turbine in Miami, Florida is calculated to be 1.3 years. 
     In order to verify the accuracy of the proposed analysis method, according to the above calculation flow, the total fatigue life of the blade coating is recalculated according to the rainfall data in the relevant experimental research made by Bech et al. And the total fatigue life of the blade coating is compared with the fatigue life calculation results in the relevant experimental research of the foreign scholars Bech et al., as shown in Table 4, the annual wind turbine life loss percentage is the annual rainfall time under each rain intensity divided by the fatigue life. Under the condition of using the same rainfall data, the expected fatigue life calculated by the method of the present application is 2.1 years, slightly longer than that obtained by Bech. This is mainly because the calculation flow proposed by the present application involves more complicated and realistic calculation methods, for example, various impact angles and raindrop shapes are considered in the simulation of a stochastic rain field. 
     
       
         
           
               
             
               
                 TABLE 3 
               
             
            
               
                   
               
               
                 Comparison of the total fatigue life in this study and from  
               
               
                 Bech&#39;s result under different rain intensities 
               
            
           
           
               
               
               
               
               
               
               
            
               
                   
                   
                   
                   
                 Annual 
                 Fatigue lives 
                   
               
               
                   
                   
                   
                 Fatigue life 
                 wind 
                 under various 
                 Annual wind 
               
               
                   
                   
                   
                 under various 
                 turbine life 
                 rain 
                 turbine life 
               
               
                   
                   
                   
                 rain 
                 loss 
                 intensities 
                 loss 
               
               
                 Rain  
                 Annual 
                 Tip 
                 intensities 
                 percentage 
                 (the result of 
                 percentage 
               
               
                 intensity 
                 rain time  
                 speed  
                 (Bech&#39;s 
                 (Bech&#39;s 
                 this method) 
                 (result of this 
               
               
                 (mm h −1 ) 
                 (h yr −1 ) 
                 (ms −1 ) 
                 result) (h) 
                 result) (%) 
                 (h) 
                 method) (%) 
               
               
                   
               
            
           
           
               
               
               
               
               
               
               
            
               
                 20 
                 1.8 
                 90 
                 3.5 
                 51 
                 4.2 
                 42.9 
               
               
                 10 
                 8.8 
                 90 
                 79 
                 11 
                 192.7 
                 4.6 
               
               
                 5 
                 88 
                 90 
                 3600 
                 2.4 
                 14463 
                 0.6 
               
               
                 2 
                 263 
                 90 
                 7.5 × 10 5   
                 3.5 × 10 −2   
                 1.6 × 10 6   
                 1.6 × 10 −2   
               
               
                 1 
                 438 
                 90 
                 2.8 × 10 9   
                 1.6 × 10 −5   
                 4.5 × 10 7   
                 9.7 × 10 −4   
               
            
           
           
               
               
               
               
            
               
                 Total percentage of annual wind turbine life loss (%): 
                 64.4 
                   
                 48.1 
               
               
                 Expected fatigue life of wind turbine (year): 
                 1.6 
                   
                 2.1 
               
               
                   
               
            
           
         
       
     
     The example effectively shows that the fatigue life of the wind turbine blade coating in a certain area can be effectively predicted by the calculation method of the present application in combination with the historical rainfall data of the area. 
     It can be understood by those skilled in the art that the above description is only the preferred examples of the present application, and is not used to limit the present application. Although the present application has been described in detail with reference to the foregoing examples, those skilled in the art can still modify the technical solutions described in the foregoing examples or replace some of their technical features equivalently. Within the spirit and principle of the present application, the modification, equivalent replacement and the like should be included within the scope of protection of the present application.