Wind turbine foundation ring fatigue reliability analysis under wind-load using SCADA system

Overview of the wind turbine engineering

ABAQUS is utilized to conduct finite element modeling and analysis of a 5 MW wind turbine generator. The main parameters of the wind turbine are as follows. The height of the generator rotor is 80 m, the diameter of the rotor is 115 m, the rated wind speed is 11.9 m/s, and the design life is 20 years. The nacelle of the turbine weighs 131.142 t. The rotor consists of three blades and a hub, with each blade’ length being 55 m. The total weight of the rotor is 68.283 t, and the rated rotational speed is 14.3 rpm. The wind turbine tower is composed of three sections of variable cross-section steel tower, with geometric parameters shown in Table 1. The wind turbine foundation adopts a circular pedestal-type expanded foundation, with a ground diameter of 18.8 m and a burial depth of 3.4 m. The tower and foundation ring are made of Q345 steel, and the foundation concrete is made of C40 concrete.

Calculation of wind speed-time series

In order to obtain real-time operational condition factors of wind turbines, most wind farms are equipped with SCADA systems, which can monitor and collect data in real time, with a data collection interval of seconds. There are limitations in obtaining data through SCADA systems, such as incomplete and inaccurate data, etc. In order to alleviate these problems, the data used in this paper is the average wind speed ν₈₀ = 16.74 m/s at the rotor of the wind turbine collected by the SCADA system during normal operation within a day.

The wind speed time series at a certain height z is equal to the average wind speed plus the fluctuating wind speed time series. The average wind speed at each height is obtained using wind speed conversion formulas. The random fluctuating wind field is calculated and simulated using the orthogonal expansion method^[22] and number theory point selection method^[35], and the wind speed time series at a certain height z is obtained accordingly. The wind speed time series is calculated by

where $$\bar{V}(z)$$ is the average wind speed at height z; $$\tilde{V}$$(θ; z, t) is the pulsating wind speed-time series at height z.

In order to simplify the calculation, the tower of a 75 m high wind turbine tower is divided into 0~15, 15~45 and 45~75 m, and the wind speed at the center of the three sections is taken as the average wind speed of the tower in the section, i.e., the average wind speed at 7.5, 30 and 60 m. The conversion for the average wind speed at each height is expressed as follows.

where λ is the ground roughness coefficient; k is the conversion coefficient of different landforms. According to the “Standard for design of high-rising structure”^[36], it can be seen that the wind farm is located in the landform of class B; the value of λ is 0.16, and the value of k is 1.0.

In the case where the power spectral density (PSD) function of the pulsating wind speed-time series $$\tilde{V}$$(θ; z, t) does not vary with height, the pulsating wind speed-time series $$\tilde{V}$$(θ; z, t) can generally be represented as^[22]

where U(θ;z) is a random field that responds to the spatial correlation of pulsating wind speeds; ν(t) is a stochastic process that responds to the pulsating wind properties.

According to the Literature^[22], it can be obtained that the stochastic process ν(t) that responds to the pulsating wind characteristics can be represented by an orthogonal expansion as

where N is the number of unfolding terms taking the value of 600; r is the truncation coefficient taking the value of 10 (The first 10 eigenvalues account for about 90% of the total energy, and the random variables of their corresponding eigenvalues can reflect the main probabilistic properties of the virtual pulsating wind displacement process^[22]); χ_n+1(n = 1, 2…, N) is a set of harmonic energy adjustment coefficients, which is mainly introduced to consider the error caused by truncation; λ_j is the j-th eigenvalue of the correlation matrix R of the stochastic process of pulsating wind speed based on the Fourier mean spectrum; φ_j,n+1 is the (n+1)-th row element of the j-th standard eigenvector Φ_j of the correlation matrix R; $$\dot{\phi}_{n}(t)$$ is the derivative of the Hartley orthogonal basis function; {ξ_j, j = 1, 2, 3…r} is a set of mutually independent standard Gaussian random vectors obtained using the number theory point selection method, in this paper r = 10.

For the correlation matrix R = [ρ_ij]_N×N, the autocorrelation function ρ_ij (only related to the time interval τ = t₂ - t₁) is given as follows^[22]

where the deterministic function ϕ_k(t) is the Hartley orthogonal basis function, where cas(t) = cos(t) + sin(t); T is the length of time which takes the value of 600 s.

Construct the corresponding autocorrelation function of the virtual pulsating wind displacement random process u(t) based on the equivalent pulsating wind characteristic random process ν(t), as given in^[22]

where $$A=6.0 k_{0} v_{10}^{2}$$; α = 3.6239$$\frac{v_{10}}{c}$$; β = γi, where γ = 2.6853$$\frac{v_{10}}{c}$$; k₀ is the dimensionless surface roughness coefficient, which takes the value of 0.0033; ν₁₀ is the average wind speed at a height of 10 m, which can be obtained from Eq. (2); the coefficient C = 600 m; $$\mathit{i} = \sqrt{-1}$$.

Because the wind speed at the center point of each tower section and at the impeller is mainly considered, for Eq. (3), there is^[22]

where {ξ_j(θ), j = 1, 2, 3…r_z} is consistent with Eq. (4), and it is also a set of mutually independent standard Gaussian random vectors; for the building within 150 m in this paper, r_z = 5^[22] can meet the required accuracy; $$\sqrt{\lambda_{z j}}$$ and f_zj(z) are the eigenvalues and eigenfunctions of the Fredholm integral equation, satisfying^[22]

The reaction pulsating wind characteristic stochastic process ν(t) and the reaction pulsating wind speed spatial correlation random field U(θ;z) where the dimensionless set of random variables {ξ_j, j = 1, 2, 3…r} and {ξ_j(θ), j = 1, 2, 3…r_z} satisfy^[22]

where δ_ij is the Kronecker-delta symbol.

If the stochastic process of pulsating wind speed and its virtual wind displacement can be assumed to be the Gaussian process, then the sets of random variables {ξ_j, j = 1, 2, 3…r} and {ξ_j(θ), j = 1, 2, 3…r_z} are both mutually independent sets of standard Gaussian random variables. Utilizing the number theory point selection method^[35], a data table that provides vectors according to Literature^[37] can be used to select points directly. In this paper, the MATLAB program is used to implement the number theory point selection method, and the number of selections in the hypersphere is N_sel = 180; i.e., 180 standard Gaussian random variable sets are obtained from the screening. Based on 180 standard Gaussian random variable groups, MATLAB software is used to write programs. Since the SCADA system collects data at one-second intervals, and wind speeds are converted to reflect the 10 m height standard in our country, we used the 10 min average annual maximum wind speed as the basic wind speed ν₀; so, the pulsating wind speed simulation is performed with a time step of 1 s, over a total simulation period of 600 s. Using Eq. (3) through Eq. (13), the appropriate pulsating wind speed-time series $$\tilde{V}$$(θ; z, t) can be obtained. Four groups of excellent pulsating wind speed-time series were selected among them, and combining with Eq. (1) and Eq. (2), further calculations were performed to obtain the wind speed-time series at 7.5, 30, 60 m, and at the impeller, i.e., the wind speed-time courses acting on 0~15, 15~45, 45~75 m and the blades of the wind turbine, and the results of the calculations are shown in Figures 1-4.

Figure 1. 0~15 m wind speed-time series.

Figure 2. 15~45 m wind speed-time series.

Figure 3. 45~75 m wind speed-time series.

Figure 4. Wind speed-time series at the impeller under normal operating conditions of the wind turbine.

Calculation of wind load-time series

The standard value of wind loads perpendicular to the surface of the wind turbine tower according to “Standard for design of high-rising structure”^[36] is calculated according to

where F_ki is the wind load standard value of each section of the wind turbine tower; ω_k is the standard value of wind load acting on the unit projection area at the height z of the wind turbine tower; ω₀ is the basic wind pressure, which takes the value of 0.55 kN/m² according to the site conditions where the wind turbine is located; h_i is the height of each section of the wind turbine tower; d_i is the diameter of the center point of each section of the wind turbine tower; β_z is the wind vibration coefficient at height z; μ_z is the wind pressure height change coefficient at height z of the wind turbine tower; μ_s is the wind load body type coefficient, taking the value of 0.7. β_z and μ_z are obtained by taking values according to the specification^[36].

Given the wind speed time series acting on the tower body, the corresponding wind load acting on the wind turbine tower body is calculated according to

where F_i is the wind load corresponding to each section of the wind turbine tower; ω_p is the wind pressure at the corresponding height z.

According to the Bernoulli's theory the relationship between wind speed and wind pressure can be obtained as

where ρ is the air density; V_t is the wind speed-time course at the impeller under normal operating conditions of the wind turbine.

For the wind load on the tower blades, it can be divided into two specific cases: normal operating conditions and extreme load conditions. Due to the obtained data in this paper being wind speeds obtained during normal operation of the wind turbine, the effects of extreme loads are not considered for now. Under normal operating conditions, the thrust force of the blade due to the wind load can be calculated according to the following equation, based on the momentum-blade element theory.

where A is the wind turbine blade swept area; C_T is the thrust coefficient. The thrust coefficients of the wind turbine are different at different wind speeds and the thrust coefficients of this 5 MW wind turbine are given by the manufacturer as shown in Figure 5.

Figure 5. Thrust coefficients of the wind turbine at different wind speeds.

The wind speed-time series of the wind turbine tower at 0~15, 15~45, 45~75 m is brought into Eq. (15) respectively, and the total wind load-time series acting on each section of the tower can be calculated out, as given in Figures 6-8.

Figure 6. Total wind load acting on 0~15 m tower body-time series.

Figure 7. Total wind load acting on 15~45 m tower body-time series.

Figure 8. Total wind load acting on 45~75 m tower body-time series.

The wind speed-time series at the impeller of the wind turbine under normal operating conditions is brought into Eq. (17) to obtain its corresponding total wind load-time series, as shown in Figure 9.

Figure 9. Total wind load acting on the impeller of the fan under normal operating conditions-time series.

Establishment of wind turbine model in ABAQUS

The focus of this study is on the foundation ring, so the model has been simplified for ease of modeling, mainly divided into three major parts: (1) upper structure; (2) tower body; and (3) foundation. The upper structure consists of the hub, blades, and nacelle, while the foundation includes the concrete foundation and the foundation ring. The tower body is modeled based on the actual geometric parameters in Table 1. Considering that the vibration characteristics of the wind turbine structure mainly depend on the tower body, the modeling of the upper structure is simplified. In order to ensure that the overall stiffness and mass remain unchanged, an equivalent modeling approach is adopted. The material properties used in the equivalent modeling have a density of 258.26 kg·m^-3, an elastic modulus of 200 GPa, and a Poisson’s ratio of 0.3. The material properties for the remaining parts are set according to their actual properties, as shown in Tables 2 and 3.

Since the focus of this study is on the fatigue reliability and fatigue life of the foundation ring, the modeling of the steel reinforcement in the concrete foundation and the modeling below the foundation are abandoned. Tie connections are used between the hub and nacelle, nacelle and tower body, and tower body and foundation ring. The concrete foundation and the foundation ring are in face-to-face contact, with Coulomb friction contact in the tangential direction and a friction coefficient of 0.35. The normal contact is treated as hard contact. Regarding the boundary conditions, the concrete foundation ground is assumed to have a fully fixed boundary condition since the modeling of the foundation is not considered. The mesh division for each component is as follows: the blades and hub use solid elements C3D10, the nacelle, concrete foundation, and tower body use solid elements C3D8R, and the foundation ring uses shell elements S4R. The overall mesh division is shown in Figure 10, and the mesh division for the focus on the foundation ring is shown in Figure 11.

Figure 10. Overall meshing.

Figure 11. Foundation ring meshing.

Since the wind speed collected in real time by the SCADA system in this study is the wind speed on the windward side, the wind load is only applied to the windward side of the turbine rotor and tower body. The total wind load received at heights of 0-15, 15-45, 45-75 m, and at the rotor is established as “amplitude” based on the corresponding relationship between time and load values. In the ABAQUS finite element software, the 75 m high tower body is divided into three sections, and reference points RP-1 to RP-3 are established at the center of each section. Each reference point is coupled with the corresponding surface of the tower body, and the loads acting on each section of the tower body are applied to the reference points. As for the total wind load-time series acting on the turbine rotor, since the main focus of the study is on the foundation ring, a concentrated load method is used to apply the load to the center of the rotor hub.

The modal analysis of the wind turbine structural model in ABAQUS

Based on the ABAQUS model established in Section "Establishment of wind turbine model in ABAQUS", in order to ensure that the wind turbine does not collapse due to resonance between low-order natural frequencies and rotor rotation frequencies (The difference between the higher order intrinsic frequency of the wind turbine and the rotational frequency of the impeller is large, and due to the role of structural damping, the higher order part of the dynamic response attenuates very quickly), modal analysis is conducted before performing transient modal dynamic analysis in ABAQUS finite element software. Therefore, the wind turbine model is first subjected to modal analysis to obtain its low-order natural frequencies.

Modal analysis is performed in ABAQUS, without considering the influence of loads, only setting boundary conditions. This study only aims to obtain the first four natural frequencies of the wind turbine, so the subspace iteration method is used for faster calculation results. The calculated natural frequencies of the wind turbine are shown in Table 4, and the modal shapes are shown in Figure 12. From the figure, it can be observed that the low-order modal shapes of the wind turbine mainly exhibit bending and swinging forms.

Figure 12. Wind turbine first 4th order mode shapes; (A) First-order mode shapes; (B) Second-order mode shapes; (C) Third-order mode shapes; (D) Fourth-order mode shapes.

The rated speed of the 5 MW wind turbine rotor is known to be 14.3 r/min. From calculations, the rotational frequency of the rotor can be determined as f = 14.3/60 Hz = 0.2383 Hz. The rotor has three blades, and it is known from literature^[38] that the prerequisite for avoiding resonance between the wind turbine tower and the rotor is that the relative difference between the natural frequency of the tower and the rotor rotation frequency or three times the rotor rotation frequency (3f = 0.7149 Hz, i.e., the passing frequency of the blades) should be greater than 10%. The relative differences between the calculated natural frequencies in Table 3 and the rotor rotation frequency are shown in Table 5, and all the relative differences are greater than 10%. Therefore, there will be no resonance between the wind turbine tower and the rotor, thus verifying the correctness of the wind turbine structure model established in this paper.

Transient modal dynamics analysis of the wind turbine model in ABAQUS

In the transient modal dynamic analysis, Rayleigh damping is considered. According to Table 3, the damping coefficients α and β corresponding to the first and second natural frequencies of the wind turbine are calculated by

where ξ is the structural vibration damping ratio, which is taken as 0.05 according to the “Standard for design of high-rising structure”^[36]; ω₁ and ω₂ are the circular frequencies of the 1st and 2nd vibration modes, respectively, and ω_i = 2πf_i, f_i is the i-th order intrinsic frequencies of the wind turbine.

By substituting the first and second natural frequencies of the wind turbine from Table 3 into Eqs. (18) and (19), the calculated Rayleigh damping coefficients α = 0.11850 and β = 0.02110 can be obtained.

The wind load calculated in Section "Calculation of wind load-time series" is applied to the wind turbine model. Then, the Rayleigh damping coefficients α and β are input into the modal dynamics analysis step of the ABAQUS finite element software to perform transient modal dynamic response analysis, thereby obtaining the stress-time series. The stress contour map of the foundation ring, which is the main focus of the study, is shown in Figure 13. By observing the figure, it can be seen that the red region in the figure, which corresponds to the flange edge of the foundation ring, is the area with the highest stress. This is consistent with the occurrence of the highest stress area in the structure in actual engineering. The stress-time series of this maximum stress region is extracted and shown in Figure 14.

Figure 13. Stress cloud of wind turbine foundation ring.

Figure 14. Stress-time series in the stress-concentrated region of the wind turbine foundation ring.

Discussion for applicability of FE model

The model used above is a 5 MW wind turbine. Because of its high power generation rate, it is more commonly used in the field of wind power generation, especially in some areas rich in wind resources (including the seaside). The dimensions and material properties of the wind turbine are obtained from the manufacturer. In performing the ABAQUS finite element modeling and analysis, the structure is modeled in a simplified manner, such as the superstructure is modeled in an equivalent manner for ease of computation and generalizability, which can be used for the study of this type of 5 MW wind turbine. This paper focuses on the fatigue reliability at the maximum stress at the edge of the flange on the base ring of a wind turbine. The fatigue reliability of the flange edge on the base ring is analyzed by the fatigue reliability calculation method considering the crossing-threshold duration, and a comparative analysis is made by the traditional method. The method is generalized and can be used for other types and sizes of wind turbine structures, and is also applicable to the calculation of stress-fatigue reliability of structures in different environments.

Stress-time histories characterized by Wiener processes

According to the “Standard for design of steel structures”^[39], the stress amplitude of lower than the fatigue cut-off limit [Δσ_L], that is, Δσ = σ_max - σ_min < [Δσ_L], will not lead to fatigue damage. Therefore, it can be concluded that fatigue damage is mainly related to the stress amplitude, and not directly related to the maximum and minimum stress values. Therefore, the stress-time series of the maximum stress region of the foundation ring is obtained by subtracting its mean value, as shown in Figure 15. At the same time, the effect of tensile-averaged stress on fatigue damage cannot be neglected because the stress average σ_m = 60.85 MPa > 0 in the region of maximum stress in the base ring in the original Figure 14 is not zero. The stress amplitude S is to be corrected equivalently in the subsequent rainfall counting method. This stress-time series has the same stress amplitude as the original stress-time series. Therefore, this new stress-time series is used for subsequent fatigue reliability analysis.

Figure 15. Stress-time series in the stress-concentration region of the modified wind turbine foundation ring.

The stress-time series of a structure under load can be considered as a stochastic process. And some stress-time histories of structures can be approximated as Wiener processes, which have the characteristics of a Wiener process. Figure 16 shows a sample curve of a Wiener process (with a duration of 1,200 and 200 s). By comparing it with the stress-time series of the new maximum stress region of the foundation ring in Figure 15, it is found that the long-term Wiener process (e.g., 1,200 s) does not match the stress-time series, but the short-term Wiener process (e.g., 200 s) is more similar to the stress-time series. At the same time, the stress variation process of the maximum stress region of the foundation ring in normal operation of the wind turbine shows similarity; i.e., the entire random stress spectrum can be regarded as a repetition of typical stress-time load spectra. Therefore, it is assumed that the short-term Wiener process can be regarded as a typical stress-time spectrum block, and based on the similarity between the short-term Wiener process and the stress-time series, the stress-time series of the maximum stress region of the foundation ring can be considered as a repetition of stress-time histories based on short-term Wiener processes.

Figure 16. The sample curves for the Wiener process.

Structural fatigue reliability analysis

Based on the progress of the analysis of multiple threshold crossings of stochastic processes in recent years^[33], the structural fatigue reliability analysis method considering the crossing-threshold duration of stochastic processes is used to analyze the fatigue reliability of structures. The calculation formula is as follows^[40].

where x₁ refers to the number of threshold crossings for the cyclic process of stresses within a periodic short time period t_short; in this paper, t_short =100 s; i is the number of periodic short time periods t_short within the normal running time length T, i.e., i = T/t_short; N_e is that the statistically obtained number of stress cycles per day is about 5.00 × 10⁶ times; E(S^m) is the m-order moment of origin of the stress amplitude S; m and C are the correlation coefficients of components and connections, according to the “Standard for design of steel structures”^[39], the components and connection types are selected as Z8, take m = 3, C = 0.72 × 10¹²; q is the length of time taken for the stress to cycle once; b is the damage threshold, according to the “Standard for design of steel structures”^[39], when the stress amplitude is lower than the fatigue cut-off limit will not lead to fatigue damage, so take the stress amplitude at times N = 1 × 10⁸ as the fatigue cut-off limit, i.e., b = $$\left[\Delta \sigma_{\mathrm{L}}\right]_{1 \times 10^{8}}$$ = 29 MPa.

The connection between the foundation ring and the tower is made using a flange connection. The maximum stress generated at the upper flange of the foundation ring can significantly reduce the fatigue strength at the upper flange and bolt locations. The research object is selected as the location of the maximum stress at the upper flange of the foundation ring, and its stress-time series is shown in Figure 17. It is now corrected for equivalent stresses, and the material strength is taken as 345 MPa. Subsequently, the rainfall counting method is applied to this data to obtain its 16-level stress spectrum, as shown in Table 6.

Figure 17. Stress-time series.

According to the data in Table 6 and the member and connection correlation coefficient m = 3, it is possible to obtain:

x₁ is the number of times the stress cycling process crosses the damage threshold b in a periodic short time period t_short, which needs to be generated by simulation with the MATLAB program. Through t_short = 100 s, it can be obtained that in 100 s, in order to accurately require to take dt = 0.01 s, through the MATLAB program simulation to generate 10,000 standard Wiener process curves. The value of the number of crossings x₁ of this sample curve with the damage threshold b was counted, as shown in Table 7.

The next step is to calculate the value on the right side of Eq. (20), i.e., calculating $$U=\frac{2 C \cdot \int_{0}^{q} \frac{b}{\sqrt{2 \pi^{3} t^{3}}} e^{\frac{b^{2}}{2 \pi}} \Gamma\left(0, \frac{b^{2}}{2 t}\right) d t}{i \cdot N_{e} E\left(S^{m}\right) \int_{0}^{q} t \frac{b}{\sqrt{2 \pi^{3} t^{3}}} e^{\frac{b^{2}}{2 t}} \Gamma\left(0, \frac{b^{2}}{2 t}\right) d t}$$^[40] to get the values as shown in Table 8.

The probability of statistically obtaining the number of crossings x₁ and the value of U in Table 8 are uniformly brought into Eq. (20). The fatigue reliability of the stress concentration area of the flange on the base ring over time can be calculated, as shown in Table 9, Figure 18. Similarly, the case of periodic short time duration of 200 s was also calculated, and its U values are shown in Table 8, and the fatigue reliability is shown in Table 10 and Figure 18.

Figure 18. The fatigue reliability curves for flange stress concentration areas on the base ring (t_short = 100 s and 200 s, b = 29).

From Figure 18, it can be observed that the fatigue reliability of the maximum stress region at the upper flange of the foundation ring decreases with time. This means that the probability of fatigue failure at the upper flange of the foundation ring grows with the service time. This also increases the probability of detachment between the foundation ring and the upper connecting structure of the wind turbine tower. Furthermore, this leads to overall instability of the wind turbine structure, resulting in alarm shutdown or collapse. Therefore, during the service life of the wind turbine structure, it is necessary to strengthen monitoring in this area and take timely reinforcement and maintenance measures to avoid accidents and economic losses within the specified service period. In addition, comparing the cases of t_short = 100 s and t_short = 200 s in Figure 18, it can be seen that when the periodic short time duration is shorter, the calculated values are more accurate and closer to reality. This is because the Winer process of the shorter periodic short time duration is closer to the stress-time series. Meanwhile, the fatigue reliability calculated by the conventional linear cumulative damage law is shown in Figure 18. Comparison with the fatigue reliability calculations for t_short = 100 s shows that they are close, but with slight differences. This is due to the fact that the method in this paper considers the effect of the duration of crossing the damage threshold and refines the effect of the stress amplitude above the damage threshold b on the fatigue damage. In contrast, the traditional method simply considers that the stress causes a fixed amount of damage to the structure every cycle, which undoubtedly leads to a large error in the calculation results^[40].