Seismic Analysis of 10-MW Off Shore Wind Turbine with Large-Diameter Monopile in Consideration of Seabed Liquefaction

12 With the increasing construction of large-scale wind turbines in seismically active coastal areas, the 13 survivability of these high-rated power offshore wind turbines (OWTs) in marine and geological 14 conditions becomes extremely important. Although research on the dynamic behaviors of OWTs 15 under earthquakes has been conducted in consideration of soil-structure interaction, attention paid to 16 the impact of earthquake-induced seabed liquefaction on OWTs supported by large-diameter 17 monopiles is limited. In view of this research gap, this study carries out dynamic analyses of a 18 10-MW OWT under the combined wind, wave, and earthquake loadings. This study uses a 19 pressure-dependent multi-surface elastoplastic constitutive model to simulate the soil liquefaction 20 phenomenon. Results indicate that the motion of the large-diameter monopile leads to more extensive 21 soil liquefaction surrounding the monopile, specifically in the zone near the pile toe. Moreover, 22 compared with earthquake loading alone, liquefaction becomes more severe under the coupled wind 23 and earthquake loadings. Accordingly, the dynamic responses of the OWT are apparently amplified, 24 demonstrating the importance of considering the coupling loadings. Compared with wind loading, the 25 effect of wave loading on the dynamic response and liquefaction potential is relatively insignificant.


31
The contradiction between the plan to increase energy consumption and reduce the carbon footprint 32 has motivated many countries to explore more sustainable and renewable energy sources, including 33 wind, solar, hydropower, geothermal, and biomass. In the past several decades, wind energy has 34 rapidly developed worldwide since the oil crisis in the early 1970s. Moreover, many onshore and 35 offshore wind farms have been constructed or are being planned, specifically in Europe, the United 36 States, and Asia. Fig. 1 presents the development history of the total installation of wind turbines 37 (WTs) in the world. The accumulated wind power capacity reached approximately 651 GW in 2019, 38 with 28.9 GW for offshore WTs, accounting for 4.5% of the total capacity (GWEC 2020).
where monopile m is the physical mass of monopile and add m is the added mass, a C is the added 167 mass coefficient, which is assumed to be 1.0 as suggested in (DNVGL 2017 PressureDependMultiYield03, and LadeDuncanMultiYield). These models were developed in the 184 UCSD based on multi-surface models suggested by Prevost (1977;1985). They presented the model 185 applications in some typical pile-supported multi-span bridges and wharves. Then, Zhang et al. (2018) 186 employed this elastic-perfectly plastic model to study the effects of peak ground velocity on the 187 pile-soil system's response on a liquefiable site. Cheng and Jeremić (2009) adopted the boundary 188 surface model presented by Dafalias and Manzari (2004) and Manzari and Dafalias (1997) and u-p-U 189 constitutive formulation to validate its effectiveness for liquefaction analysis. Rahmani and Pak (2012) 190 further demonstrated the efficiency of this model by comparing the numerical results and 191 experimental outcomes by Wilson (1998   As shown in Fig. 3 and Table 1, the length and width of the soil domain are assumed to be 100 m, 211 which is ten times the monopile diameter, whereas the depth is 55 m, approximately 1.5 times the pile 212 length. This depth is adopted per the suggestion of Barari et al. (2017) to achieve a good balance 213 between computing efficiency and accuracy. The soil domain is modeled by eight-node BrickUP 214 elements with four DOFs for each node, wherein DOFs 1-3 are for node displacements and DOF 4 is 215 for fluid pressure. The mesh size in the vertical direction is set to 5.0 m, whereas different mesh sizes 216 (2.5, 5.0, and 10.0 m) are adopted in the horizontal direction, as suggested by Chiaramonte et al. 217 (2013). Fig. 3 shows the detailed mesh information, where fine and coarse meshes are used in the 218 domains close to and far from the monopiles, respectively. A simplified approach, the rigid 219 beam-column links (radially rigid "spikes") connecting the pile and soil nodes, is employed to model 220 the interaction between the monopile and surrounding soil. This method was initially adopted by Law 221 and Lam (2001)  Eigenvalue analyses are conducted to examine the dynamic properties (i.e., modal frequencies and 238 mode shapes) of the OWT system. The obtained modes 1, 2, and 6 with the frequencies of 0.287, 239 1.256, and 3.091 Hz are corresponding to the first, second, and third tower fore-aft frequencies, 240 respectively. The first frequency is slightly higher than the value of 0.25 Hz reported by Velarde 241 (Velarde 2016). The discrepancy may be attributed to different SSI simulation approaches. The 242 simulated fundamental frequency in this study still lies in between the rotational frequency range 1P 243 (0.104-0.167 Hz) and blade passing frequency range, 3P (0.312-0.501Hz). 244

System damping 245
The damping mechanism of an OWT is quite complicated. The damping usually comprises 246 structural, aerodynamic, hydrodynamic, and soil damping, which accounts for the contributions of the 247 superstructure, wind, wave, and current, and soil, respectively. Aerodynamic damping results from the 248 interaction between the wind and rotating blades and varies with operational conditions. In the parked 249 state, the blades are pitched to the maximum pitch angle, and the aerodynamic damping is minimal 250 and negligible. By contrast, the aerodynamic damping is high in normal operating conditions. The 251 aerodynamic damping ratio in the fore-aft direction for an in-operation WT is generally in the range of Rayleigh damping is the most common model in the RHA, and the following equations could be 267 adopted to estimate mass and stiffness coefficients: 268 where α and β are the mass and stiffness coefficients, and a  and b  are two selected frequencies 269 with the target damping ratio ζ. The fundamental frequency 1  of structures and the predominant 270 frequency ep  of ground motions, as shown in Table 3, are suggested by QUAD4M to be used to 271 consider excitation characteristics, where QUAD4M is a computer program to evaluate seismic 272 responses of soil structures (Hudson et al. 1994). 273

276
In IEC 61400-3-1 (IEC 2019), five different types of loads should be included for the design 277 calculation: (1) gravitational and inertial loads, including static and dynamic loads resulting from 278 gravity, vibration, rotation, and earthquake; (2) aerodynamic loads, involving static and dynamic loads 279 caused by the airflow and its interaction with the stationary and moving parts of WTs; (3) actuation 280 load, resulting from the control system of WTs; (4) hydrodynamic loads caused by water flow and its 281 interaction with the support structure; and (5) ice loads acting on the foundation. As this study focuses 282 on the influence of earthquakes on the liquefaction risk, the dynamic complexity arising from the ice 283 loads and control system (e.g., pitch and yaw system) is not included. 284

Inertial loading 285
As listed in Table 1, the structural weight is composed of the masses of different components, 286 including the masses of RNA, tower, monopile, and transition section. In this study, the tower and 287 monopile are assumed to interface at the mean sea level. The water surrounding the monopile is 288 modeled by the added mass method with a coefficient of 1.0, as described in subsection 2.2. The soil 289 domain is simulated by the continuum elements, in which the gravitational forces are added 290 automatically in OpenSees. 291

Wind loading 292
The wind load along the tower is applied to the tower nodes in the x-direction. The thrust force 293 imposed on the hub is extracted from the Fatigue, Aerodynamic, Structures, and Fatigue (FAST) code 294 (Jonkman and Jr. 2005), which was developed by NREL in the US. A regular wind profile at a rated 295 wind speed is considered, which is more likely to occur during the entire lifetime of WTs and is more 296 realistic than storms and hurricanes when an earthquake occurs. 297 The wind action along the tower is proportional to the wind velocity profile, as described by the 298 following formula in IEC 61400-3-1 (IEC 2019). 299 where hub h is the height of the hub center measured from the sea surface, hub v is the velocity at the 300 hub height and is assumed to be the rated wind speed (i.e., 11.4 m/s), and z and v(z) are the height and 301 wind velocity of the concerned tower node, respectively. Then, the velocity could be employed to 302 calculate the horizontal forces through the following equation: 303 where a  is the air density, typically equal to 1.225 kg/m 3 , and ( ) A z is the tributary area.

304
Two methods are usually adopted to evaluate the resultant aerodynamic force. One is a simplified 305 method considering the thrust coefficient. In this approach, Eqs. (6) and (7) where T C is the thrust coefficient, a  is the air density, hub v is the velocity at the hub height, 308 swept A is the rotor swept area, and hub V is the rated wind speed. 309 The other is the blade element momentum (BEM) theory considering the appropriate wind speed 310 information. This method incorporated in the FAST is adopted in the following simulation, wherein 311 the Kaimal spectrum is employed to produce the wind field, and the BEM theory is utilized to 312 calculate the aerodynamic loading along the blade. The Kaimal spectrum for three wind components 313 (K = u, v, and w) (Jonkman and Buhl 2006) is shown below. 314 where f is the wind frequency, K L is an integral scale parameter, hub u is the mean wind speed 315 corresponding to the hub height HubHt , U  is the turbulence scale parameter, and K  is the 316 ambient turbulence standard deviation in three directions. In addition, the spatial coherence model is 317 adopted to add the correlations between the same wind components at two spatially separated points. 318 The IEC equations for the three components (Jonkman and Buhl 2006) shown below are considered. 319 0.12 12, 5.67 * min(60 m, ) where f is the wind frequency, d is the distance between two concerned spatial points i and j, and 320 where r is the radial distance of the considered element from the center of the hub, where is the circular frequency in rad/s, is the significant wave height, is the peak wave 339 frequency, and parameter * is calculated by 340 where is the peak enhancement factor (generally taken as 3.3), and is a shape parameter. When 341 ≤ , = 0.07; while = 0.09 for > . In the wave load simulation, the peak wave 342 frequency and the significant wave height are adopted as 0.96 rad/s and 1.9 m, respectively. 343 After the PSD of the sea surface elevation is determined, the sea surface elevation time history can 344 be simulated by using the inverse fast Fourier transform (IFFT) technique. Therefore, the 345 hydrodynamic load on each monopile element can be calculated by using the Morison formula (DNV 346 2014). 347 where w  is the seawater density; representative wave loads at the mean sea level and mudline are shown in Fig. 5b. 352

Earthquake loading 353
The medium sand is considered to reproduce the liquefaction phenomenon, and its corresponding 354 shear wave velocity is 198.7 m/s, which corresponds to the site class D defined in the ASCE 7-10 355 (ASCE 2010). Therefore, five earthquake records for site class D, denoted as E1 to E5, are selected, 356 with their PGAs ranging from 0.18 to 0.48 g. Table 3 Fig. 6 depicts the acceleration response spectra for a damping ratio of 3.0%, wherein the 361 damping ratio corresponds to an operating condition of OWT. The response spectra for the damping 362 ratio (i.e., 1.5%) in a parked state are not shown herein. All the peak frequencies of the five response 363 spectra are considerably lower than the tower's first fore-aft frequency, but they are closer to the 364 second or third fore-art frequencies of the tower. In particular, the Northridge earthquake shows two 365 peaks close to the second and third fore-art frequencies. The simulation time is 5 s beyond the ground motion records to consider a dissipation time. 375

Staged simulation procedure 378
A staged modeling procedure is employed to simulate the pile existence and generate accurate 379 initial conditions for liquefaction assessment (Wang et al. 2016a). In the first stage, the soil domain 380 without the pile is modeled by adopting the BrickUP elements with a significant permeability 381 coefficient and tied boundary condition. An elastic gravity stage is initially carried out. Then, the state 382 of soil is updated using the proposed plastic model to produce the initial effective stress in scenario 1. 383 Next, the superstructure, monopile foundation, and links are introduced, and three translational DOFs 384 of the monopile and link nodes are tied to adjacent soil nodes. Another plastic gravity step is 385 conducted to account for the settlement and consolidation of soil caused by the pile and WT gravity 386 and generate the proper initial stress for liquefaction evaluation in scenarios 2, 3, and 4. After the 387 gravity stages, the earthquake loading or the coupled loadings are applied. The soil permeability 388 coefficient is updated with the appropriate value to simulate the undrained condition of the soil 389 domain. Fig. 7 depicts the detailed procedure for scenarios 2, 3, and 4. 390

393
The pore water pressure is monitored at seven columns of soil points, including SP1×, SP2×, SP3×, 394 SP4×, and SP5× on the left side (upwind) of the monopile and SP6× and SP7× on the right side 395 (downwind) of the monopile, to assess the liquefaction possibility of the soil domain. The "×" symbol 396 represents the numbers 1-5 for different soil layers in a column. The vertical displacements at points 397 DV1 and DV2 are monitored to evaluate the soil settlement and monopile tilt. The acceleration and 398 displacement at the tower top are also observed, denoted as AT and DT, respectively, to examine the 399 dynamic responses of the OWT system. Fig. 8 shows the detailed layout of the monitoring points. could also be employed to evaluate the foundation bearing capacity. Finally, the EPP, which is defined 409 as the difference between the instantaneous pore water pressure and the hydrostatic pore pressure, at 410 the monitoring points SP11, SP12, and SP13 (on the left side of monopile) is presented to illustrate 411 the EPP accumulation and assess the liquefaction distribution. 412 As the dynamic response and liquefaction development under five selected earthquake records 413 follow a similar trend, the following subsections mainly focus on the Kobe earthquake to explain the 414 dynamic responses and liquefaction potential. This event causes the most severe liquefaction 415 distribution and capacity deterioration. 416  The dynamic response and liquefaction severity are summarized in Table 6. In case 1, the tower top 430 displacement, rotation angle, and maximum bending moment in the Y direction are nearly zero, as 431 both the wind and earthquake loadings are applied in the X direction. In case 2, the dynamic responses 432 in the Y direction intensify with the introduction of the earthquake in this direction. However, the 433 maximum dynamic responses in Table 6 indicate that case 1 with the consistent wind and earthquake 434 direction represents the severest situation. Therefore, the discussions in the following sub-sections are 435 mainly focused on this most critical case. Comparing the results of the half model (shown in Table 5) 436 and the full model -case 1 (shown in Table 6) verifies the accuracy of the half model when the wind 437 and earthquake loading are applied in the same direction. 438 439  The liquefaction areas could be judged by the excess pore water pressure ratio (EPPR) u r , 458 where u(z) is the EPP at a depth of z and ' 0 ( ) of the selected earthquakes (Fig. 6). Compared with other seismic ground motions, the predominant 476 frequency of the Kobe earthquake is closer to the fundamental frequency (0.9 Hz) of the free-field soil, 477 which intensifies the dynamic responses and consequently leads to more severe liquefaction. 478

OWT in the parked condition (scenario 2) 479
The OWT system supported on a large-diameter monopile under earthquakes is simulated in 480 scenario 2. Fig. 10b shows the EPP variation of the soil nodes SP1× at different depths under the 481 Kobe (i.e., E3) earthquake. The development history of EPP exhibits a similar trend as shown in the 482 free field. However, the variation amplitude is more significant in scenario 2 due to the cyclic 483 tension-compression stress induced by the monopile motion. In particular, the EPP around pile toe 484 (i.e., the depth of Z = −35 m) in Fig. 10b accumulates more quickly than that in the free field, as 485 shown in Fig. 10a. Furthermore, the maximum EPP values are approximately 68.9 and 174.7 kPa in 486 scenarios 1 and 2, respectively, which indicates greater stress transfer in this zone in scenario 2. 487 Fig. 11b shows the EPPR distribution under the Kobe earthquake in scenario 2 when the EPPR of 488 SP13 (Z = −15 m) reaches the maximum value. At the selected instant, the soil inside and outside the 489 pile shows different states. Considering the confinement of the pile wall, the EPPR of soil inside the 490 monopile (X = 0 m) is smaller than that outside (X = ±7.5 m), and correspondingly, more severe 491 liquefaction can be found outside of the monopile. Fig. 12a shows the EPP variation in the horizontal 492 direction at a selected depth of Z = −15 m, which corresponds to the maximum liquefaction depth 493 (Table 5). The EPP variation is affected apparently by the monopile within the distance of 12.

OWT in the operating condition (scenarios 3 and 4) 533
Scenarios 3 and 4 consider the OWT operating at a rated wind speed when an earthquake occurs 534 with and without including wave loadings, respectively. As shown in Fig. 14a, the maximum 535 displacement at the tower top in scenario 3 reaches 1.317 m, which is 5.9 times the peak value in 536 scenario 2 but is slightly smaller than the limit required in GB 50135-2019. The maximum value in 537 scenario 4 is slightly larger, around 1.392 m. It is noteworthy that turbulent wind loading has been 538 applied before the earthquake attacks, and thus the initial displacement is not zero in two scenarios. 539 Moreover, Fig. 14b  where the corresponding EPP value is larger due to the introduction of wave loading. Fig. 15a directly 549 compares the EPP development histories between scenarios 2 and 3 at a selected depth of Z = −15 m 550 corresponding to the maximum liquefaction depth in these two scenarios. A noticeable difference can 551 be observed between scenarios 2 and 3 for the soil nodes next to the pile wall (X = −7.5 m). 552 Following the seismic record peak at 8.46 s, the coupling effect of the wind and earthquake loading 553 has a significant influence on EPP accumulation and leads to a larger value. In the remaining time, 554 EPP values in scenario 2 are slightly greater. The effect of wave loading on the EPP accumulation is 555 relatively limited when comparing the EPP values in scenarios 3 and 4, as shown in Fig. 15b. The 556 EPPR distribution in scenarios 3 and 4 when the EPPR of SP13 (Z = −15 m) reaches the maximum is 557 shown in Fig. 11c and Fig. 11d, and similar liquefaction distributions could be found. Parkfield earthquakes (E2, E4, and E5), while the smaller EPP could be found in the other two seismic 578 records. 579 However, the effects of the coupled earthquake and wind loadings (scenario 3) become more 580 evident in the pile rotation angles and bending moments, as shown in Fig. 16b and Fig. 16c,  581 respectively. The increment ratios range from 38.5% to 128.2% under the coupled earthquake and 582 wind loadings for the maximum bending moment, and the increment ratios of the pile rotation angle at 583 the mudline range from 44.8% to 171.4%. The liquefaction severity has a significant impact on the 584 bending moments and rotation angles of the monopile. Given that the liquefaction depth is only 5 m 585 under the Parkfield (E5) earthquake, the bending moment and rotation angle are relatively small. By 586 contrast, with the largest liquefaction depth under the Kobe (E3) earthquake, bending moment and 587 rotation angles are maxima among the five earthquake cases. The significant liquefaction weakens the 588 stiffness and capacity of the soil foundation and thus amplifies the dynamic responses of the monopile 589 and superstructure. When including the wave loading in scenario 4, the rotation angle and bending 590 moment are not increased except for the Kobe earthquake, which corresponds to the severest 591 liquefaction condition. The increment ratios of the rotation angle and maximum bending moment 592 under the Kobe eaerthquake are 7.5% and 4.9%, respectively. 593  3, as shown in Fig. 17a. The corresponding values in scenario 2 are only 69.88 kPa and 1.853%, 602 respectively, which are not shown in this study for brevity. The relatively larger shear stress amplitude 603 and accumulated strain under the coupled loading intensify soil performance deterioration in scenario 604 3. Fig. 17b also indicates that the mean effective stress gradually reduces to zero with the loading 605 cycles during the earthquake, indicating the liquefaction process. In comparison to 10-m depth, 606 although the soil node at 25 m depth in Fig. 17c experiences larger shear stress, the node is not finally 607 liquefied given the positive mean effective stress. 608 A similar stress path could be found in scenario 4 considering the wave loading, but the 609 corresponding shear stress amplitude (104.44 kPa) and maximum shear strain (1.904 %) at shallow 610 depth (Z=-10 m) are smaller than those in scenario 3. For the sake of brevity, the diagrams for 611 scenario 4 are not shown in Fig. 17. (1) Compared with free field soil, the large-diameter monopile considerably increases the EPP of 627 the soil surrounding the pile, particularly at the pile toe. However, the influence of the pile becomes 628 limited for soil with a horizontal plan distance beyond twice the pile diameter. 629 (2) The motions of the OWT under earthquake loadings increase the liquefaction depth and 630 intensify the liquefaction severity. 631 (3) Compared with earthquake loadings, the coupled wind and earthquake loadings have small 632 effects on the EPP values of soil. However, the introduction of wind loading leads to increased shear 633 stress amplitude and accumulated shear strain, which has a critical influence on soil degradation. 634 (4) The liquefaction reduces the foundation stiffness and capacity of the OWT. Consequently, 635 under the coupled wind and earthquake loadings, the pile rotation angles at the mudline and tower top 636 displacements increase distinctly. In particular, the maximum pile rotation angle exceeds the limit 637 specified in a Chinese standard. 638 (5) Considerably larger bending moment in the monopile is observed under the coupled wind and 639 earthquake loadings. The peak bending moment occurs with the liquefaction, with the maximum 640 value appearing close to the interface of the liquefiable and non-liquefiable soil layers. 641 (6) Compared with the wind loading, the influence of wave loading on the dynamic response (tower 642 top displacement, pile rotation angle, and bending moment) and the liquefaction severity is limited. 643