Aerodynamic Force and Comprehensive Mechanical Performance of a Large Wind Turbine during a Typhoon Based on WRF/CFD Nesting

Featured Application

The main research conclusions of this paper could provide references for load prediction of similar wind turbine systems during the occurrence of typhoon, and deepen the understanding on meso/microscale wind field nesting mechanism.

Abstract

Compared with normal wind, typhoons may change the flow field surrounding wind turbines, thus influencing their wind-induced responses and stability. The existing typhoon theoretical model in the civil engineering field is too simplified. To address this problem, the WRF (Weather Research Forecasting) model was introduced for high-resolution simulation of the Typhoon “Nuri” firstly. Secondly, the typhoon field was analyzed, and the wind speed profile of the boundary layer was fitted. Meanwhile, the normal wind speed profile with the same wind speed of the typhoon speed profile at the gradient height of class B landform in the code was set. These two wind speed profiles were integrated into the UDF (User Defined Function). On this basis, a five-MW wind turbine in Shenzhen was chosen as the research object. The action mechanism of speed was streamlined and turbulence energy surrounding the wind turbine was disclosed by microscale CFD (Computational Fluid Dynamics) simulation. The influencing laws of a typhoon and normal wind on wind pressure distribution were compared. Finally, key attention was paid to analyzing the structural response, buckling stability, and ultimate bearing capacity of the wind turbine system. The research results demonstrated that typhoons increased the aerodynamic force and structural responses, and decreased the overall buckling stability and ultimate bearing capacity of the wind turbine.

1. Introduction

A large wind turbine structure is more flexible [1] and has a more prominent wind-induced dynamic effect [2]. It is a typical wind-sensitive structure. Wind load is the control load of its structural design [3,4]. In particular, the wind-induced failure of large wind turbine structures occur frequently during strong typhoons [5,6]. For example, the Typhoon “Cuckoo” in 2003 affected the normal operation of wind power plants in south China, which halted 13 wind turbines. The Typhoon “Saomai” in 2006 landed on Zhejiang, and caused different degrees of failure of 20 units in the Cangnan Wind Power Plant. In 2014, the Typhoon “Rammasun” landed at Hainan Province, which caused the collapse of—and serious damage to—the 1.5-MW wind turbine with a 77-m blade diameter in the wind power plant. Compared with normal wind, typhoons form a more complicated near-ground wind field. The wind field characteristics of high turbulence, changeable direction, and great changes of shearing wind speed will intensify the fluctuation of the wind turbine structure. Under this circumstance, airflow movement on the near-wall surface is disorganized, and the surface pressure is changed significantly, which further caused significant changes in the wind-induced responses and the stability of the wind turbine system [7,8]. For these reasons, a contrastive analysis of aerodynamic performances and comprehensive performances of the large wind turbine system during typhoons and normal wind was of important theoretical significance and engineering values.

For studies concerning wind turbine performance during typhoon, Takeshi et al. [9] studied influences of blade, wind direction, and the relative position of a doorway on the collapse of wind turbines according to wind turbine collapse data and measured typhoon speed in Japan. Tang et al. [10] and Chou and Tu [11] analyzed stresses on the wind turbine structure under typhoon load by the simplified static analysis method. Li et al. [12] analyzed typhoon characteristics and the damages of different parts of wind turbines. They pointed out that the design of different parts of a wind turbine has to strengthen the structural strength and pay attention to the yaw of the wind turbine under typhoon load. Wang et al. [6] and Luo et al. [13] carried out a systematic analysis on wind-induced static/dynamic responses of the tower in the wind turbine system during a typhoon by combining a theoretical deduction and fluctuating wind spectrum. Some suggestions regarding the design of the typhoon resistance of the tower body in the wind turbine system were given. Lian et al. [14] and Zhou et al. [15] carried out a systematic analysis on the aerodynamic load on blades of the wind turbine system under typhoon loads by combining a theoretical analysis and measured data. They established a typhoon-resistance aerodynamic blade design optimization model for wind turbines. Cao [16] and An [17] analyzed the flow field and pressure field and structural stresses of the wind turbine under typhoon loads by theoretical studies and CFD numerical simulation. A new typhoon-resistance wind turbine design philosophy and structural design was proposed. Utsunomiya et al. [18] and Ma et al. [19] carried out a statistical analysis on the dynamic response and mooring linear tension response of a submarine platform of an offshore wind turbine under typhoon loads. They attached high attentions to the aerodynamic shape design and typhoon resistance characteristics of blades.

Most of the above studies concerning the typhoon resistance of a wind turbine are based on theoretical analysis and measured data. The theoretical systems in these studies were too simplified to reflect the specific wind profile and landing attenuation effect of a mesoscale typhoon field. Nowadays, mesoscale models have been well applied in case analyses of typhoons, early warnings, and disaster forecasts. Common models include Weather Research Forecasting (WRF), MM5 (Mesoscale Model Version 5), and ETA. Among these models, MM5 enjoyed widespread applications in earlier studies, because it was studied and used earlier. However, it has been never adopted by the NCEP (National Centers for Environmental Prediction), on the grounds that its dynamical framework is obsolete, and its procedures are not fairly standard. Utilized by the NCEP for operational forecast, ETA can hardly promptly assimilate the outstanding research findings of scientific research institutes and universities, although its simulation results are reliable. As a result, its promotion is restricted. Comparatively, the WRF model based on fluid dynamics and thermodynamics can not only effectively simulate wind speed, pressure, and temperature in a typhoon field, it can also comprehensively consider the evolution process, strong specificity, and attenuation characteristics of typhoons. In spite of its parametric schemes and the more complicated setting of some of its parameters, it takes more complete and accurate environmental parameters into account. Therefore, WRF is used most widely among all of the existing models. However, due to the large sphere of scope (hundreds of kilometers) of typhoons, grid resolution of the typhoon field is generally at the magnitude of kilometers. However, the overall dimension of a large wind turbine system is only at the magnitude of hectometers. To predict the aerodynamic load at blade edges accurately, it has to deepen into the near-ground boundary layer. The minimum grid size of the near-wall surface is generally lower than 10⁻² m. The mesoscale model in WRF fails completely at this magnitude. Hence, it will be helpful for observing the airflow motion in a typhoon field by a more refined simulation of microscale typhoon fields based on the mesoscale WRF. WRF and Computational Fluid Dynamics (CFD), which are new-generation mesoscale and microscale forecast models respectively, are employed most widely as technologies for high-resolution and refined wind field forecasts. In the literature [20,21], refined numerical simulations were conducted on the bridge sites of mountainous areas by multi-scale coupling. The results suggested that the results of WRF could be more adaptable to inlet boundaries during numerical simulations of these wind fields by CFD. Thus, the existing problem about the mean wind speed at the inlet during the numerical simulations was successfully solved. In the literature [22,23,24], the complicated wind fields of lakes on the underlying surface were simulated meticulously by coupling WRF and CFD. According to the simulation results, WRF–CFD coupling better eliminates the mean deviation and the root mean squared error than WRF. Besides, it reflects the atmospheric characteristics of actual microscale topographical changes more accurately. Therefore, the meso/microscale nesting WRF/CFD has to be used for the high-precision simulation of a multi-scale flow field of large wind turbines under typhoon loads in order to solve various problems, such as the high-precision transition of parameters and flowing structure, multilayer and multi-scale grid nesting, multi-time scale control, and multi-scale mutation [25].

To solve these problems, the mesoscale WRF model was introduced for high temporal spatial resolution simulation of Typhoon “Nuri”. After the typhoon field information was analyzed, the wind speed profile in the simulation region was acquired based on the nonlinear least square fitting. On this basis, the normal wind speed profile with the same wind speed of the typhoon speed profile at the gradient height of the class B landform was set. Later, the aerodynamic distribution on the surface of a five-MW horizontal axis wind turbine in a wind power plant in Shenzhen (China) under typhoon loads and normal wind was studied by microscale CFD numerical simulation. Finally, the mechanical performance and stability performance of a large wind turbine system under typhoon loads and normal wind were compared by combining them with the finite element method.

2. Mesoscale Typhoon Field Simulation and Analysis

2.1. WRF Model

WRF is the mesoscale weather research and forecasting system that was developed by the America Environmental Forecasting Center, National Center of Atmospheric Research of the United States and Rainstorm Analysis and Forecasting Center of University of Central Oklahoma [25]. The Euler equation model based on a non-static force equilibrium was used as the WRF–ARW (Advanced Research WRF) model of the dynamic framework. The horizontal computational domain adopted the Arakawa-C meshing scheme [26] (Figure 1a). The Arakawa-C grid can express scalar and vector parameters simultaneously, which is conducive to increase the accuracy of a high-resolution simulation. Besides, the Arakawa-C grid has good dispersion properties and conservativeness. The vertical computational domain applied the Euler mass coordinate η along the terrain (Figure 1b):

$$\eta =\left(p_h-p_{ht}\right)/\mu ,\mu =p_{hs}-p_{ht}$$

(1)

where p_h is the air pressure at the target point p_ht is the air pressure of the top layer, which is a constant, and p_hs is the near-ground air pressure.

The WRF model is transferable when it is run on the Linux system. With considerations to the influences of physical processes such as water vapor, long-wave radiations, short-wave radiations, cumulus cloud, and the underlying surface, the WRF model can simulate airflow, air pressure, and wind field characteristics in a large region. The results can be used as the input boundary conditions for CFD numerical simulation.

2.2. Selection of Physical Scheme and Parameter Settings

The #12 Typhoon “Nuri” in 2008 was chosen as the simulation object. Based on the information provided by the Typhoon Website of China Meteorological Administration, it landed on coastal regions of Saigon, Hong Kong on 16:55 on 22 August and the maximum wind near the center was 12-level (33 m/s). On the same day, it landed again at Nanlang Town, Zhongshan City, Guangdong Province at 22:10, and the maximum wind near the center was eight-level (25 m/s). For effective numerical simulation of typhoon path as well the air pressure field and wind speed field, the simulation covers the whole process from the first landing at Saigon, Hong Kong to the second landing at Zhongshan City, Guangdong Province. The initial atmospheric conditions and time-dependent boundary conditions, which were applied in the WRF model, were all based on the FNL (Final Global Data Assimilation System) [27,28,29,30]. It covered 27 atmospheric layers. The sampling interval was set at 6 h, and the spatial resolution was 1° × 1°. The WRF simulation region had 37 vertical layers, and the air pressure at the top layer was set 5000 Pa.

The simulation center chose a district in north Shenzhen. Typhoon “Nuri” ran through this region at 20:00 on 22 August. With comprehensive considerations to data demands and computational conditions, the three-dimensional (3D) nesting program was used to calculate typhoon field in the WRF model. The horizontal mesh resolution was 13.5 km, 4.5 km, and 1.5 km, respectively. The longitude and latitude of center of the simulation region was (114.1° E, 22.5° N). The coarse grid space D01 in the external surface was 13.5 km, and the number of grids was 211 × 211. The grid space D02 in the second layer was 4.5 km, and the number of grids was 217 × 217. The grid space D03 in the inner layer was 1.5 km, and the number of grids was 241 × 241. The map projection applied the Lambert program. A total of 37 vertical layers were set, forming upper-sparse and lower-dense stratification patterns. There were 19 layers below 1000 m. The simulation domain is shown in Figure 2. The altitudes of different vertical layers are shown in

Table 1. Altitudes of different vertical layers.

Layers	Altitude/m	Layers	Altitude/m
0	0	19	997.85
1	5.41	20	1441.03
2	21.64	21	1864.44
3	43.28	22	2383.35
4	64.96	23	2929.90
5	86.67	24	3508.54
6	108.42	25	4124.07
7	130.20	26	4780.99
8	152.01	27	5484.65
9	179.34	28	6243.02
10	217.70	29	7067.05
11	272.71	30	8489.99
12	350.13	31	10,070.22
13	439.15	32	11,332.60
14	528.81	33	12,805.37
15	619.11	34	14,592.36
16	710.08	35	17,014.78
17	824.83	36	19,617.33
18	987.29	37	——

.

The hub height of the wind turbine is generally lower than 150 m, which is below the thickness of the boundary layer. The hub height is significantly influenced by airflow mass, humidity, and heat transmission in the boundary layer [31,32,33,34,35]. Therefore, the settings of the physical parameterization scheme can affect the simulation accuracy of boundary layer of the typhoon field directly. In this paper, the physical parameters of WRF model are listed in

Table 2. Physical parameter settings of the WRF model.

WRF Parameters	Main Zone (D01)	Nesting Zone (D02)	Nesting Zone (D03)
Horizontal resolution	13.5 km	4.5 km	1.5 km
Integral time step	180 s
Microphysical process scheme	Lin
Long-wave radiation	RRTM
Short-wave radiation	Dudhia
Near-ground layer scheme	Monin-Obukhov
Land surface process scheme	Noah
Planet boundary layer scheme	MYJ
Cumulus convective parameterization scheme	Kain-Fritsch

.

2.3. Simulation Results and Discussions

To verify the validity of the proposed simulation model of Typhoon “Nuri”, the typhoon path and minimum sea level pressure throughout the simulation are shown in Figure 3. The simulated path was drawn according to the minimum pressure in the cyclonic center of the WRF model. The measured path was the best path provided by the Typhoon Website of the China Meteorological Administration, and it could also refer to some literature cites [36,37]. The simulated path and the measured path were stable, and both paths moved toward the northwest. In early simulations, the simulated path was close to the measured path. However, the typhoon path deflected gradually as time went on, which might be caused by the intensifying initial error with the increase of simulation time. The evolution of typhoon intensity was mainly represented by the minimum sea level pressure. The simulated value and measured value of minimum sea level pressure in the simulation period differed slightly, but the distribution trend was basically consistent. Both the simulated value and the measured value increased slowly with time, and increased sharply in the late simulation period. At this moment, the typhoon became weak and disappeared gradually.

Simulation results and measured results in the simulation regions at typical height for the landing of the typhoon at 20:00 on 22 August are shown in Figure 4 and Figure 5. These nephograms are conducive to analyzing the wind field information when the Typhoon “Nuri” ran through the center of the simulation region. It can be concluded from the analysis that:

(1)

The typhoon had a low-pressure center. The land underlying the surface was rougher than the ocean after landing, thus cutting the heat source. Hence, the upflow was weakened, while the surrounding airflows continued radiation toward the typhoon center, which gradually further increased the air pressure. After comparison, it was found that the simulated pressure distribution agreed well with the measured results.

(2)

Wind speed increased due to the invasion of the peripheral nephsystem of the typhoon. The wind speed increased continuously as the typhoon approached. The wind speed closer to the center was the higher.

(3)

Influenced by the downward transmission of atmospheric heats, the MYJ boundary layer scheme contributed a high computational accuracy of turbulent fluctuation fluxes, such as heats and momentum in the boundary layer of the region. This implied that the low atmosphere heats were mainly carried into the boundary layer by the upper atmosphere through high temperature, which caused the hybrid heating of the low atmosphere.

(4)

Cloud precipitation close to the typhoon belongs to the long-term convective precipitation. The large rainfall volume in this region is related to the typhoon strength and the central distance. The Noah pavement scheme can simulate the surface flux and airflow convergence field well.

(5)

There is a minor difference between the simulation results and the measured results of the typhoon, but the distribution trend is basically consistent. In the central area of the typhoon simulation, the ranges of wind pressure, wind speed, temperature, and rainfall were 860–920 Pa, 22–26 m/s, 294–297 K, and 40–48 mm/h, respectively. These data are close to the measured data on the website of the National Meteorological Center. The maximum error is lower than 15%.

The wind speed streamlines before, at, and after the landing of the typhoon are shown in Figure 6. No strong convection was observed before the landing of the typhoon. However, it turned to a southeaster gradually at the landing of the typhoon. After the landing, surrounding airflows continued to radiate toward the typhoon center. Wind speed streamlines are disorganized, without a regular development trend.

The wind speed profiles close to the typhoon center at different moments are shown in Figure 7a. According to analysis, wind speed distributes uniformly in the near-ground region at different moments, but wind speed distributes disorderly at high altitude. The typhoon speed at 20:00 on 22 August was significantly lower than those at other moments. This might be because the typhoon had landed at Shenzhen at this moment, but it was on the sea at other moments. The near-ground wind speed in the center of the simulation region and the corresponding nonlinear least square fitting curves are shown in Figure 7b. The fitting value of the wind speed profile index was 0.118, which was smaller than the surface roughness index (0.15) of class B landform under normal wind, indicating the good fitting effect in the near-ground typhoon field (simulation accuracy reached 95.97%). The typhoon speed at 10 m of height was high, and increased slowly with the increase of height. Moreover, it defined that the typhoon field at 350 m altitude of class B landform in the code [38] was equal to the wind speed in the normal wind field for the purpose of qualitative and quantitative analyses of the difference between the normal wind field and the typhoon field. The near-ground typhoon profile and the normal wind profile were integrated into the UDF and used as the initial boundary conditions in the follow-up microscale CFD numerical simulation.

3. Microscale CFD Numerical Simulation

3.1. Brief İntroduction to the Project

The main structural design parameters and model of a five-MW windward horizontal axis three-blade wind turbine in a wind power plant in Shenzhen, China are shown in

Table 3. Main structural design parameters and model of the five-MW wind turbine.

Parameters	Numerical Value	3D Blade Model	3D Wind Turbine Model
Tower height	124 m
Radius at tower top	3.0 m
Radius at tower bottom	3.5 m
Thickness at tower top	0.04 m
Thickness at tower bottom	0.09 m
Blade length	60 m
Nacelle size	18 m × 6 m × 6 m
Cut-out speed	25 m/s
Yaw rotation speed	0.5°/s

. This five-MW wind turbine was developed in the Jiangsu Key Laboratory of Hi-Tech Research for Wind Turbine Design of Nanjing University of Aeronautics and Astronautics [39]. The tower body was a long thickness-variable structure. The upper wall thickness was 40 mm, and the bottom wall thickness was 90 mm. The nacelle size was 18 m (Length) × 6 m (Width) × 6 m (Height). The dip angle of the blade was 5°, and the cut-out speed was 25 m/s. The blades were distributed uniformly along the circumferential direction at an angle of 120°, and the blade length was 60 m. The parameters of the different blade sections along the span are listed in

Table 4. Blade parameters of wind turbine.

Position/%	Blade Span/m	Chord Length/m	Installation Angle/°	Blade Pitch Angle/°	Position/%	Blade Span/m	Chord Length/m	Installation Angle/°	Blade Pitch Angle/°
5	3	2.9	0.823	37.14	55	33	1.95	0.169	−0.293
10	6	3.66	0.64	26.672	60	36	1.75	0.156	−1.072
15	9	4.41	0.507	19.069	65	39	1.58	0.144	−1.736
20	12	4.56	0.414	13.692	70	42	1.42	0.134	−2.31
25	15	4.25	0.346	9.83	75	45	1.27	0.125	−2.81
30	18	3.91	0.296	6.976	80	48	1.12	0.118	−3.25
35	21	3.59	0.258	4.802	85	51	0.98	0.111	−3.64
40	24	3.05	0.229	3.103	90	54	0.83	0.105	−3.987
45	27	2.63	0.205	1.742	95	57	0.69	0.099	−4.299
50	30	2.29	0.186	0.63	100	60	0.54	0.095	−4.58

. The airfoil types of the blade were NH02_40, NH02_35, NH02_30, NH02_25, NH02_21, NH02_18, and NH02_15. Given the coordinates, chord lengths, and installation angles of the discrete points on the different blade profile sections in the wind turbine blades, the three-dimensional (3D) coordinates of different blade profile sections were gained through coordinate transformation. Finally, the 3D model of the blade was constructed.

Specific steps of coordinate transformation included [40]:

(1)

The coordinates of discrete points on different blade profile sections were calculated from relevant software. The normalized coordinates (x₀, y₀) of different blade elements were transformed to the coordinates (x₁, y₁) of discrete point under the coordinate system 1:

$$\left(x_1,y_1\right)=c\times \left(x_0,y_0\right)$$

(2)

(2)

The coordinates were further transformed to the coordinates (x₂, y₂) under the coordinate system 2:

$$(x_2,y_2)=(x_1,y_1)-(X,Y)$$

(3)

where (X, Y) are the coordinates of the pressure center.

(3)

The coordinates were further rotated to the actual spatial coordinates (x, y, z):

$$x=\sqrt{x_2^2+y_2^2}\times \cos \left(\arctan \frac{y_2}{x_2}+\theta \right)$$

(4)

$$y=\sqrt{x_2^2+y_2^2}\times \sin \left(\arctan \frac{y_2}{x_2}+\theta \right)$$

(5)

$$z=r$$

(6)

where θ is the installation angle, and r is the turning radius of the blade elements.

The inflow angle is the included angle between the plane of rotation and inflow speed in the traditional theory of blade momentum. However, the calculation formula of the inflow angle during the blade design stage of the wind turbine according to the Baez limit is:

$$\cot \phi =\frac{3}{2}\lambda$$

(7)

where λ is the tip velocity ratio.

According to the relationship between the tip velocity ratio and the power coefficient, the optimal tip velocity ratio can be determined, and the inflow angle of each blade element is:

$$\cot \phi =\frac{3}{2}{\lambda }_r$$

(8)

where λ_r is the optimal tip velocity ratio.

The tower body, blade, and nacelle models were constructed successively. Finally, the 3D physical model of a large wind turbine was formed by Boolean operation.

3.2. Computational Domain and Meshing

To assure the full development of the wake flow of the wind turbine [40,41], the computational domain was determined as 12D × 5D × 5D (flow direction X × extension direction Y × vertical direction Z, and D is the diameter of the wind wheel). The wind turbine was placed 3D to the entrance of the computational domain. To achieve satisfying computational efficiency and accuracy [42,43], due to the complicated blade appearance, the hybrid grid discrete form was applied to divide the whole computational domain into the internal and external parts. The local encrypted region covered the wind turbine model, which used the non-structured meshing. Non-structured meshing produces a tetra/mixed tetrahedral mesh in ICEM CFD automatically by Robust (Octree). The peripheral region had a regular shape that used the high-quality structured meshing. The grid number, grid quality, and windward pressure coefficient of the tower body under different meshing schemes are shown in

Table 5. Grid number, grid quality, and windward pressure coefficient on the tower body under different meshing schemes.

Meshing Schemes	1	2	3	4	5
Total number of grids	1.1 million	4.5 million	7.4 million	9.3 million	28.4 million
Minimum orthogonal quality of grids	0.13	0.36	0.53	0.60	0.64
Grid skewness	0.95	0.87	0.82	0.74	0.71
Windward pressure coefficient	0.92	0.88	0.85	0.80	0.79

. Grid quality includes the minimum orthogonal quality of the grids and grid skewness. The minimum orthogonal quality of the grids shall be larger than 0.1, and it is best to be higher than 0.2. No negative volume is allowed. The grid skewness shall be lower than 0.95, and it is suggested to be lower than 0.9 [44]. It can be seen from Table 5 that with the increase of the total number of grids, the grid quality and windward pressure coefficient converged gradually. No significant differences were observed in grid quality and computational results under a 9.3 million meshing program and a 28.4 million meshing program. Considering computational accuracy and efficiency, this paper chose the 9.3 million meshing program. When processing flow in the near-wall region, the standard wall function was applied to connect the physical variables on the wall and the unknown variables in the turbulence core region directly. In this paper, the cooperation of the standard wall function and standard k-ε model (which will be introduced in Section 3.3) can increase the simulation accuracy more significantly. Essentially, the flowing conditions in the turbulence core region were solved by the k-ε model. The flow in the wall surface was not solved. Instead, the physical parameters (e.g., speed) of this region were solved by a semi-empirical formula directly. Since the near-wall flow was processed by the standard wall function in this paper, no encryption of the wall surface was needed during meshing, and it only has to place the first internal node in the logarithmic law region. The calculated results demonstrated that the y+ of the model wall could assure bottom grids in the logarithmic law, which met common engineering requirements. The meshing schemes of the overall computational domain and model are shown in Figure 8. Both the number of grids and grid quality met the requirements of calculation.

3.3. Selection of Turbulence Model

The 3D steady implicit algorithm was used in this paper. The air wind field chose the incompressible flow field and the turbulence model chose the standard k-ε model [45]. The speed vector and pressure coupling problem in the momentum equation was solved by the SIMPLEC algorithm. The turbulent kinetic energy, turbulent dissipation term, and momentum equation in the calculation all adopted the second-order upwind discretion, and the calculation residual error of the governing equation (Navier-Stokes(N-S) equation) was set 1 × 10−6. During the simulation, the enhanced wall surface function model was started to assure that the logarithmic distribution of the bottom grids was true.

The continuity equation and momentum equation of the viscous incompressible N-S equation based on Reynolds average are described by Equations (9)–(13).

$$\frac{\partial {\mu }_i}{\partial x_i}=0$$

(9)

$$\frac{\partial (\rho \overline{{\mu }_i})}{\partial t}+\frac{\partial (\rho \overline{{\mu }_i}\overline{{\mu }_j})}{\partial x_j}=-\frac{\partial (\overline{p})}{\partial x_i}+\frac{\partial }{\partial x_j}(\mu \frac{\partial \overline{{\mu }_i}}{\partial x_j}-\rho \overline{{\mu }_i^{\prime }}\overline{{\mu }_j^{\prime }})$$

(10)

where i, j = 1, 2; ρ is the air density, which is generally 1.225 kg/m³L and μ is the dynamic coefficient of viscosity, which is generally 1.7894 × 10⁻⁵ kg/(m·s)

$$\frac{\partial ({\rho }_{m}k)}{\partial t}+\frac{\partial ({\rho }_{m}k{\mu }_j)}{\partial x_j}=P-{\rho }_m\epsilon +\frac{\partial }{\partial x_j}[(\mu +\frac{{\mu }_t}{{\sigma }_k})\frac{\partial k}{\partial x_j}]$$

(11)

$$\frac{\partial ({\rho }_m\epsilon )}{\partial t}+\frac{\partial ({\rho }_m\xi {\mu }_j)}{\partial x_j}=C{\xi }_1\frac{\xi }{k}{P}_t-C{\xi }_2{\rho }_m\frac{{\xi }^2}{k}+\frac{\partial }{\partial x_j}[(\mu +\frac{{\mu }_t}{{\sigma }_k})\frac{\partial \xi }{\partial x_j}]$$

(12)

$${\mu }_t=\frac{C_{\mu }{\rho }_m{k}^2}{\xi }$$

(13)

where k and ε are the turbulence energy and the turbulence dissipation rate, $P_t$ is the generating term of turbulence energy, and ${\mu }_t$ is the viscosity coefficient of turbulence. The model constants are $C_{\zeta 1}$ = 1.44, $C_{\zeta 2}$ = 1.92, ${\sigma }_{\zeta }$ = 1.3, ${\sigma }_k$ = 1.0, and $C_{\mu }$ = 0.09.

Under the premise of meeting the above mass and momentum conservation equations, the basic flow characteristics of turbulence are described by Equations (14)–(16).

$$\frac{\partial H}{\partial t}+\frac{\partial (\mu H)}{\partial x}+\frac{\partial (vH)}{\partial y}+\frac{\partial (wH)}{\partial z}=0$$

(14)

$$\frac{\partial (\mu H)}{\partial t}+\frac{\partial }{\partial x}\left|u^{2}H+\frac{1}{2}g{H}^2\right|+\frac{\partial (\mu vH)}{\partial y}=Hfv-gH\frac{\partial Z_b}{\partial x}-\frac{H}{\rho }\frac{\partial p}{\partial x}+\frac{{\tau }_{ax}-{\tau }_{bx}}{\rho }$$

(15)

$$\frac{\partial (vH)}{\partial t}+\frac{\partial (\mu vH)}{\partial x}+\frac{\partial }{\partial y}\left|v^{2}H+\frac{1}{2}g{H}^2\right|=-Hfu-gH\frac{\partial Z_b}{\partial y}-\frac{H}{\rho }\frac{\partial p}{\partial y}+\frac{{\tau }_{ay}-{\tau }_{by}}{\rho }$$

(16)

where u, v, and w are components of the speed vector u(z) on the x, y, and z directions; H is the reference height; g is the gravitational acceleration; w is the angular velocity of ground rotation; $\phi$ is the geological latitude of the calculation point; $Z_b$ is the elevation of ground; ${\tau }_{ax}$ and ${\tau }_{ay}$ are the wind stress components along the x and y directions; ${\tau }_{bx}$ and ${\tau }_{by}$ are the ground frictional force components along the x and y directions; p is the pressure difference; and $\rho$ is the air density. f can be expressed as:

$$f=2w\sin \phi$$

(17)

The k-ε model is a turbulence model that is applicable to the calculation of high Reynolds numbers. It can predict the fully developed turbulence flow well and it is suitable to the simulation of the surrounding flow fields of similar megawatt wind turbines [46].

3.4. Settings of Boundary Conditions

The inlet of the computational domain was used as the speed inlet boundary, and the outlet was the pressure outlet boundary. Two side walls and the top surface used the symmetric boundaries. The wind turbine and ground set the wall boundary. The overlapping surface of the local and peripheral computational domain was set as the interface. The wind field computational domain and its boundary conditions are shown in Figure 9. Wind velocity profiles were set as the indexes of the atmospheric boundary layer in the computational domain, which are based on the profiles of typhoon and normal wind that are shown in Figure 7b of the revised manuscript, respectively. The velocity profiles of the typhoon and normal wind were determined by formulas (18) and (19) as follows:

$$V_{\mathrm{T}}=16.49\times {\left(h/10\right)}^{0.118}$$

(18)

$$V_{\mathrm{N}}=14.70\times {\left(h/10\right)}^{0.15}$$

(19)

where, V_T and V_N were velocity of typhoon and normal wind at the height of h. They equaled to 16.49 and 14.70 respectively at the height of 10 m. 0.118 was the index on the velocity profiles of typhoon fit by nonlinear least squares, while 0.15 was the roughness index of surface on Type B landforms as specified by Chinese Code [38].

In this study, turbulence of wind fields near the ground was determined by Formulas (20) and (21) under typhoon and normal wind respectively.

$$I_{\mathrm{T}}=0.20\times {(h/10)}^{-0.118}$$

(20)

$$I_{\mathrm{N}}=0.14\times {(h/10)}^{-0.15}$$

(21)

where 0.20 was the nominal turbulence of the typhoon “Nuri” at the height of 10 m, which was determined based on the synchronous monitoring results of the typhoon [37], and 0.14 was the nominal turbulence on Type B landforms at the height of 10 m, as specified by Chinese code [38].

In typhoon field simulation, fluid parameters (e.g., mean wind profile, scale of turbulence, turbulence energy, turbulence integral scale, and specific dissipation rate) in the fluctuating wind field were integrated into Fluent by UDF based on the above k-ε model and fitting typhoon profile (Figure 10). This was to accomplish the simulation of the typhoon field in a meso/microscale nesting large wind turbine.

3.5. Analysis of Aerodynamic Results

The surface pressure coefficient on the wind turbine tower is a main index for evaluating the aerodynamic performances of a wind turbine system. Figure 11 shows the typical cross-section circumferential curve of the pressure coefficient for a wind turbine tower under normal wind conditions, and compares the coefficient with the code value [38]. The analysis suggested that the numerical simulation results of the wind turbine tower under normal wind conditions were generally in accord with the circumferential distribution laws and numerical values specified by Chinese codes. They were only lower than the value specified by the codes out of wind, but the maximum error was below 10%. This demonstrated that the numerical simulation results of CFD were effective.

The wind speed streamlines and turbulence energy distribution on typical sections of the tower body under normal wind and typhoon loads are shown in Figure 12 and Figure 13. H was the height of the tower body. The wind speed streamlines presented relatively consistent distribution under normal wind and typhoon loads. The wind speed surrounding the tower body and blade under the typhoon loads was higher than that under normal wind. The inflow split at 0° of the windward side of the tower body, thus generates reflux and different sizes of eddies on the leeside of the tower body. Additionally, the inflow moved along the blade surface and ran through the front and rear edges of the blades, which formed small-sized eddies on the blade leeside. At the same time, the turbulence energy surrounding the tower body and blades under typhoon loads was significantly stronger than that under normal wind. Wakes on the blade and tower body developed the large-scaled eddy increment region. Besides, the interaction between the blade and the tower body was significant on the middle and upper sections of the tower body.

The distribution curves of the wind pressure coefficient at the typical section of the tower body under normal wind and typhoon loads are shown in Figure 14. Under normal wind and typhoon loads, the circumferential wind pressure coefficient on the tower body was symmetric along the windward side. With the increase of circumferential angle, it decreased firstly and then increased until it stabilized close to the leeside. However, extremum positive pressure on the windward side and extremum negative pressure on the lateral surface under the typhoon were significantly higher compared with those under normal wind. The maximum increase of the wind pressure coefficient on the windward side was 61.1%, and the maximum increase of the extreme negative pressure on the lateral surface was 44.9%. Both were achieved at bottom of the tower body.

The straight-up blade was defined Blade A. The remaining two clockwise rotating blades were Blade B and Blade C. Wind pressure coefficient nephograms on the windward side and leeside of blade as well as mean along the blade span under normal wind and typhoon loads are shown in Figure 15. The distribution curves of the overall wind pressure coefficient of the blades under normal wind and typhoon loads are shown in Figure 16. It can be seen that:

(1)

The wind pressure on the windward side of all of the blades was positive under normal wind and typhoon loads. Negative wind pressure only occurred at the front blade edges and rear edges of the blade root. However, the wind pressure on the leeside of blades was basically negative.

(2)

The wind pressure on the windward side of Blade A generally increased along the blade span. However, the positive pressure on the windward sides of Blade B and Blade C changed slightly along the span, without irregular changes. Negative pressure on the leeside of the blades decreased firstly, and then increased along the span direction. If these two effects were superposed, the overall wind pressure coefficient of the different blades decreased firstly, and then increased along the span direction.

(3)

Positive wind pressure on the windward side of the blades under typhoon loads was higher than that under normal winds. The negative wind pressure on the leeside was slightly different from that under normal winds. As a result, the overall wind pressure coefficient on different blades under typhoon loads was higher than that under normal winds. The maximum increase was 23.6%, which occurred at the tip of Blade C.

The distribution curves of the lift coefficient and drag coefficient on the different height sections of the tower body under normal wind and typhoon loads are shown in Figure 17. The lift–drag coefficient ratio is shown in

Table 6. Lift–drag coefficient ratio of the tower body under normal wind and typhoon loads.

CL/CD	Height of the Tower Body
	12.4 m	24.8 m	37.2 m	49.6 m	62 m	74.4 m	86.8 m	99.2 m	111.6 m
Normal wind	−0.040	−0.071	−0.072	−0.055	−0.011	−0.074	−0.107	−0.064	−0.179
Typhoon loads	0.024	−0.019	0.001	0.001	0.042	−0.010	−0.028	−0.002	−0.078

. The calculation formulas of the across-wind lift coefficient (C_L) and downwind drag coefficient (C_D) are shown as follows [4]:

$$C_L=\frac{{\displaystyle \sum _{i=1}^n{C}_{Pi}{A}_i\sin {\theta }_i}}{A_T}$$

(22)

$$C_D=\frac{{\displaystyle \sum _{i=1}^n{C}_{Pi}{A}_i\cos {\theta }_i}}{A_T}$$

(23)

where $Cpi$ is the mean wind pressure coefficient at the measuring point i on the tower body, A_i is the pressure coverage area at the measuring point i, θ_i is the included angle between the pressure at measuring point i and the wind axis, and A_T is the projection area of the overall structure along the wind axis.

The following can be seen from Figure 17 and Table 6. (1) The lift coefficient was negative under normal wind, but it was alternatively positive and negative under typhoon loads. However, the distribution law of the lift coefficient was consistent under normal wind and typhoon loads, and the numerical value was small; (2) The drag coefficient was mainly in the range of 0.35~0.4 under normal winds, and it was in the range of 0.45~0.65 under typhoon loads. The drag coefficient at different heights under the typhoon loads increased dramatically than that under the normal wind. The maximum was about 0.18, which was at the tower top. This was mainly caused by the large disturbance of the tower top to blades.

4. Response Analysis of the Wind Turbine Coupling Structure

4.1. Finite Element Modeling and Dynamic Characteristics

The tower–blade integrated finite element coupling model was constructed by the large universal finite element analysis software ANSYS. The tower body and blade were simulated by Shell63 unit, whereas the nacelle and internal structure were simulated as a whole by using Beam 189 units. The round raft basis unit was Solid65. The interaction between the foundation and the basis was simulated by a Combin14 unit. After finite element models of different wind turbine parts were constructed, a coupling connection may cause the relative slippage of models under loads, which was attributed to the different unit types of different parts. Therefore, the multipoint restriction unit coupling command was used to connect different parts and construct the tower–blade integrated simulation model. According to the principle of efficiency and the accuracy balance, the model was divided into 4122 units.

The natural frequency of vibration and the mode of vibration of the wind turbine were solved by a block Lanczos algorithm after modeling [47,48]. The finite element model of the wind turbine, the first 100 orders of the frequency distribution, and the mode of vibration at typical orders are shown in Figure 18 and Figure 19. The tower–blade integrated model had a low natural frequency of vibration. The fundamental frequency was 0.137 Hz, and the first 10 orders of frequency were lower than 1 Hz. At the same time, the natural frequency of the vibration of the structure was low and distributed densely. The structural frequency achieved approximately linear growth with the number of vibration orders, but the growth rate decreased after 90 orders. In addition, three blades in the tower–blade integrated model made front-to-back movement at the first, fifth, and 10th order. The blades made complicated front-to-back and left-to-right swings at the 30th order and 50th order, which were accompanied by the bending deformation of the tower body. According to the multi-order mode of vibration, the low-order mode of vibration in the tower–blade integrated model was mainly manifested by the front-to-back and left-to-right swing of the blades. Under the high order of vibration, the tower body and blades developed deformations and buckling failure.

4.2. Stress Analysis on the Tower Body

The 3D distributions of the radial displacement of the tower body and distributions of internal stress responses at the tower bottom under normal wind and typhoon loads are shown in Figure 20 and Figure 21. It can be seen that: (1) the radial distribution of the tower body presented relatively similar distribution patterns under normal wind and typhoon loads. The radial displacement increased gradually going up of the tower body, forming two extreme regions on the windward side and leeside. The maximum positive and negative displacements were at 180° and 0°; (2) Different wind field mainly affected the radial displacement at the middle and upper parts of the tower body. Under the typhoon loads, the maximum radial displacement of the lower tower body was 0.117 m, which was increased by 28.6% compared with that under normal wind; (3) Internal stresses at the tower bottom were symmetric along the windward side under normal wind and typhoon loads. Despite the slight difference in the meridian axial force at the tower bottom, other internal stresses presented basically consistent circumferential distribution laws under normal wind and typhoon loads. All of the internal stresses achieved the maximum value close to the circumferential 330° of the tower body; (4) Extreme internal stresses at the tower bottom under typhoon loads were stronger than those under normal wind loads. The maximum meridian axial force, maximum shearing force, maximum circumferential bending moment, and maximum meridian bending moment were increased by 36.1%, 34.8%, 34.2%, and 46.5%, respectively. All occurred close to the circumferential 330°.

4.3. Stress Analysis of Blades

The internal stress responses at blade roots and downwind displacement distribution along the blade span under normal wind and typhoon loads are shown in Figure 22 and Figure 23. Downwind displacements at the blade tips are shown in

Table 7. Downwind displacements at blade tips under normal wind and typhoon loads.

Blade A		Blade B		Blade C
Normal Wind	Typhoon	Normal Wind	Typhoon	Normal Wind	Typhoon
−0.701 m	−1.047 m	−0.599 m	−0.885 m	−0.596 m	−0.884 m

. It was found that: (1) under normal wind and typhoon loads, the internal stresses at the root of Blade A were significantly higher than those of Blade B and Blade C. The numerical values of Blade B and Blade C were similar. Moreover, the internal stresses at the blade roots under typhoon loads were significantly higher than those under normal wind loads. The maximum increase of the shearing force was 41.3%. The maximum increase of the circumferential bending moment was 39.6%. The maximum increase of the meridian bending moment was 39.5%. They were all achieved by Blade C; (2) The downwind displacements of Blade B and Blade C along the span were basically same under normal wind and typhoon loads. However, they were significantly different from those of Blade A. With the extension along the blade span, the downwind displacement of Blade A was higher than those of Blade B and Blade C, and then became smaller. Next, the numerical values that were close to the middle of the blades were the same. However, the numerical values that were close to the middle of Blade A exceeded those of Blade B and Blade C. The numerical values of the three blades were equal at about 35 m of the blade span under normal wind, and at about 30 m of the blade span under typhoon loads; (3) The downwind displacements of the blades were larger under typhoon loads than those under normal wind loads. The maximum increase of downwind displacement was 49.3%, which was achieved at the tip of Blade A.

5. Wind-Induced Stability Analysis of Wind Turbine System

5.1. Buckling Stability Analysis

The buckling mode and eigenvalues of the tower-blade integrated model under normal wind and typhoon loads are shown in

Table 8. Comparison of buckling mode and eigenvalues of wind turbine system under normal wind and typhoon loads.

Buckling Eigenvalue under Normal Wind		Buckling Mode under Normal Wind	Buckling Mode under Typhoon Loads	Buckling Eigenvalue under Typhoon Loads
Buckling coefficient	2.955			Buckling coefficient	2.992
Critical wind speed/(m/s)	24.12			Critical wind speed/(m/s)	31.57
Maximum displacement/m	0.025			Maximum displacement/m	0.032

. Based on comparison, the critical wind speed was higher than the designed wind speed under normal wind and typhoon loads, meeting the design requirements. The buckling position of the wind turbine was at the tower top. The buckling displacement and critical wind speed of the large wind turbine system were significantly influenced by different wind fields. The buckling coefficient, critical wind speed, and maximum displacement under typhoon loads were increased by 1.3%, 30.9%, and 28% compared with those under normal wind loads.

5.2. Ultimate Bearing Capacity

Changes in the buckling displacement of the wind turbine system with wind speed under normal wind and typhoon loads are shown in Figure 24. It was found that with the leveling loading of the wind speed, the buckling displacement of the wind turbine system under normal wind and typhoon loads increased firstly, and then decreased, and finally increased slightly. It reached the peak at about 20 m/s of the reference wind speed, and reached the minimum at about 40 m/s of the reference wind speed. The above trend was due to the “reverse effect” [46,49]. Moreover, the typhoon effect would significantly decrease the ultimate bearing capacity of the wind turbine system. The buckling displacement was increased by 38.5% at 40 m/s of the reference wind speed.

6. Summary and Conclusions

High-resolution simulation of the Typhoon “Nuri” was accomplished based on the mesoscale WRF model. Later, differences regarding the aerodynamic force, wind-induced responses, buckling stability, and ultimate bearing capacity of the wind turbine system under normal wind and typhoon loads were analyzed by microscale CFD numerical simulation and the finite element method. Some major conclusions could be drawn:

• A high-resolution simulation of Typhoon “Nuri” was carried out through the mesoscale WRF model. The typhoon profile index that was based on the least square fitting was 0.118, which was about 20% smaller than the wind profile index under normal wind loads. In this paper, the meso/microscaled nested downscaling method can effectively simulate a three-dimensional typhoon field of similar structures and provide effective input loads for follow-up comprehensive stress analysis.
• The flow field distribution surrounding the wind turbine is basically similar under normal wind and typhoon effects. However, the airflow under typhoon conditions showed higher wind speed and turbulence compared with those under normal wind. Positive pressure on windward surface of the tower body, negative pressure on the crosswind surface of the tower body, the overall pressure coefficient of the blade, the resistance coefficient of the tower body, and the coupling structural response all increased significantly. The buckling stability and ultimate bearing capacity of the coupling system decreased significantly.
• The pressure coefficient on the windward surface of the tower body under typhoon conditions increased by 61.1% to the maximum extent. The maximum increment of the negative pressure extremum on the crosswind surface was 44.9%. Both were observed at the bottom of the tower body. The typhoon influenced the pressure coefficient on the windward surface of the blades significantly, but it only affected the leeside surface slightly. The maximum increment of the overall pressure coefficient of the blade was 23.6%, which occurred at the tip of Blade C. Under typhoon conditions, the resistance coefficient of the tower body increased significantly, and the maximum increment, 67.5%, occurred at the bottom of the tower body.
• Under typhoon conditions, the radial displacement at the tower top of the wind turbine increased by 28.6%. The increments of the maximum meridian axial force, maximum shearing force, maximum circumferential bending moment, and maximum meridian bending moment were 36.1%, 34.8%, 34.2%, and 46.5%, respectively. All were observed close to the circumferential 330°. The maximum increments of the shearing force, circumferential bending moment, and meridian bending moment at the blade root under typhoon conditions were 41.3%, 39.6%, and 39.5%, which were observed on Blade C. The maximum increment of downwind displacement at the blade tip occurred on Blade A, which reached 49.3%.
• Under typhoon conditions, the buckling coefficient, critical wind speed, and maximum buckling displacement of the wind turbine system increased by 1.3%, 30.9%, and 28%, respectively. The wind turbine system achieved the maximum buckling displacement when the basic wind speed was about 20 m/s, and the minimum buckling displacement at about 40 m/s of the basic wind speed. When the basic wind speed was about 40 m/s, the buckling displacement under typhoon conditions increased by 38.5%.