Seismic soil–structure interaction analysis of wind turbine support structures using augmented complex mode superposition response spectrum method

. In this study, the seismic soil–structure interaction (SSI) of wind turbine support structures is investigated based on the complex mode superposition approach. For accurate and efﬁcient estimation of seismic loadings on wind turbine support structures, an augmented complex mode superposition response spectrum method (RSM) is developed, where the maximum shear force and bending moment of the non-classically damped system are analytically derived. An empirical formula of the modal damping ratios with a threshold value for the allowable damping ratio is also proposed to improve the prediction accuracy of the shear force acting on the footing. Furthermore, additional loadings to consider the contribution of the mass moment of inertia of rotor and nacelle assembly and P − (cid:49) effect to the bending moment on the tower are analytically derived. The proposed formulae are ﬁrst demonstrated upon a 2 MW wind turbine supported by two different types of foundations. A parametric study is then carried out by changing the tower geometries and soil conditions to propose the threshold value for the allowable damping ratio.


Introduction
In recent years, the expansion in wind energy has increased the construction of wind turbines in seismically active regions, including Japan, and damages on wind turbine support structures caused by huge earthquakes have been reported. A piled foundation was damaged at Kashima wind farm during Great East Japan Earthquake (March 11, 2011) (Ashford et al., 2011) 20 and a wind turbine tower was buckled at Kugino wind farm during Kumamoto Earthquake (April 16, 2016) (Harukigaoka Wind Power Inc., 2016). In order to assure the structural integrity of wind turbine support structures against such huge earthquakes, an accurate and efficient method to estimate seismic loadings on wind turbine support structures is thus desired to be developed.
The response spectra of the structures with low damping ratios show large fluctuations, whereas damping correction factors 30 in the above codes cannot accurately capture the uncertainty in the response spectra due to such fluctuations. To cope with this issue, Ishihara et al. (2011) have proposed a new damping correction factor for wind turbine support structures, in which the uncertainty in the response spectra is considered by employing a quantile value. Moreover, Kitahara and Ishihara (2020) have recently extended it to be applicable to MW class large sized wind turbine support structures as well.
Moreover, the mass ratio between the super-and sub-structures of wind turbines is very different, and the footing mass 35 can reach about six times total masses of the tower, rotor, and nacelle, particularly in areas of intense seismic activities (Ishihara ed., 2010). Therefore, seismic responses of wind turbine support structures are severely affected by soil-structure interaction (SSI) (Wolf, 1989;Zhao et al., 2019). An efficient method to account for the effect of SSI is to represent the soilstructure system as a seismic SSI model with a set of springs and dashpots at the soil-foundation interface, substituting the soil-foundation system. This method has been widely employed for the seismic analysis of wind turbine support structures 40 (Bazeos et al., 2002;Butt and Ishihara, 2012;Stamatopoulos, 2013;Kitahara and Ishihara, 2020). However, by introducing the dashpots, the seismic SSI model is a non-classically damped system and thus classical modal damping models, such as Rayleigh damping model, cannot accurately represent its modal damping properties. Although the modal damping ratios of non-classically damped systems can be theoretically obtained by the complex eigenvalue analysis, they cannot be directly used for RSM, since complex eigenmodes do not generally coincide with real eigenmodes, which is used in conventional 45 RSM. To overcome this obstacle, Kitahara and Ishihara (2020) have proposed a modal decomposition method to identify equivalent modal damping ratios for the real eigenmodes from the modal damping ratios obtained by the complex eigenvalue analysis, and estimated seismic loadings on wind turbine support structures by conventional RSM using the identified modal damping ratios.
On the other hand, Gao et al. (2019) have recently proposed complex mode superposition RSM for the non-classically 50 damped systems to estimate maximum relative displacements based on the complex eigenmodes. This method employs the modal damping ratios obtained by the complex eigenvalue analysis, and thus it does not require the additional procedure to estimate the equivalent modal damping ratios for the real eigenmodes in Kitahara and Ishihara (2020). In this study, seismic loadings, i.e., maximum shear forces and bending moments, will be newly derived by this method to estimate those on wind turbine support structures. Gao et al. (2019) have demonstrated this method upon 3-and 6-story shear structures. However, 55 in these structures, highly damped modes are not dominant in their seismic responses, and thus its accuracy in case of highly damped modes have not been investigated. Meanwhile, considering wind turbine support structures, a mode corresponding to the sway motion of the footing is dominant in the shear force acting on footings, and this mode shows significantly large damping ratio due to the soil radiational damping in case of soft soil profiles.
The aim of this study is consequently to further investigate the applicability of complex mode superposition RSM to 60 wind turbine support structures, where highly damped modes are dominant in their seismic responses, and to propose an accurate and efficient method to estimate seismic loadings acting on towers and footings. To achieve this objective, a seismic SSI model of wind turbine support structures is constructed, in which the effect of SSI is considered by a pair of springs and https://doi.org/10.5194/wes-2021-81 Preprint. Discussion started: 3 December 2021 c Author(s) 2021. CC BY 4.0 License. dashpots in the sway and rocking directions, respectively. Seismic loadings on this model are first derived by Complex mode superposition RSM. To improve the prediction accuracy of the shear force acting on footings, complex mode superposition 65 RSM is then augmented by introducing the upper limit of modal damping ratios. In addition, the bending moment at the hub height due to the mass moment of inertia of rotor and nacelle assembly (RNA) is accounted for as an additional loading by the angular acceleration at the hub height. The proposed method is validated by time history analysis (THA) accounting for different types of foundations and different tower geometries.

Wind turbine support structures under earthquake 70
Wind turbine support structures subjected to a horizontal earthquake motion are investigated in this study. As the foundation type, the gravity and piled foundations are considered. In section 2.1, a sway-rocking (SR) model is constructed as the seismic SSI model of wind turbine support structures. In this model, the effect of SSI is considered by a pair of springs and dashpots in the sway and rocking directions, respectively. In addition, an input acceleration response spectrum is defined accounting for the effects of soil amplification and damping correction in Section 2.2. In Section 2.3, seismic loadings on 75 this model is derived by complex mode superposition RSM, and this method is then augmented to improve the prediction accuracy of the shear force acting on footings. The relationship between the structural damping ratios of steel towers and their natural frequencies, as well as additional loadings due to RNA and P-Δ effect are finally summarized in Section 2.4.

Seismic SSI model for wind turbine support structures
The seismic SSI model of wind turbine support structures is constructed as the SR model shown in Fig. 1. Here, x and z 80 express the horizontal and vertical coordinates, respectively. The (= 1, ⋯ , − 1) th node represents each degree of freedom (DOF) of the steel tower and footing. In this model, RNA is simplified as a lumped mass at the hub height ( = ) and is connected to the tower by using a rigid beam. This simplification has been verified not to influence the prediction accuracy of seismic loadings acting on the tower and footing, excluding that the bending moment at the hub height is underestimated because of no consideration of the mass moment of inertia of RNA (Kitahara and Ishihara, 2020). This 85 bending moment though enables to be obtained as an additional loading by the angular acceleration at the hub height, which is further investigated in Section 2.4. In addition, the tower and footing are modeled by lumped masses and Euler-Bernoulli beam elements. The number of beam elements is set as = 27, which fulfills the requirement in the guidelines for design of wind turbine support structure and foundations by Japan Society of Civil Engineers (JSCE) (Ishihara ed, 2010).
As the foundation type, the gravity and piled foundations are considered. The gravity foundation is typically employed 90 in case of stiff soil profiles, whilst the piled foundation is necessary to be installed in case of soft soil profiles. Regardless of the foundation type, the soil-foundation system is substituted with a pair of springs and dashpots in the sway and rocking directions, respectively, and is connected to the footing bottom. It should be noted that, for simplification, the frequency dependence of the springs and dashpots, cross-coupling between the sway and rocking springs, and mass moment of inertia of the foundation are neglected in this model. The slightly embedded footing with a small embedment ratio, which is defined as the ratio of the footing depth to width, is employed in this study, and it ensures the representation of SSI by the uncoupled sway and rocking springs and dashpots is a reasonable simplification. The dynamic finite element equation of the seismic SSI model with respect to relative motions can be written as: 100 at each DOF of the tower relative to the footing; is the unit column vector; ̈0 is the input acceleration time histories at the footing bottom.
For the gravity foundation, the stiffness constants and damping coefficients of the soil-foundation system in Eq. (1), , , , and , can be obtained using the cone model as detailed in Architectural Institute of Japan (AIJ) (2006). For the 110 piled foundation, on the other hand, the stiffness constants of the springs in the sway and rocking directions can be calculated by Francis and Randolph models, respectively (Francis, 1964;Randolph, 1981), whereas the damping coefficients of the dashpots can be obtained by Gazetas model (Gazetas and Dobry, 1984). The detailed derivation of these stiffness constants and damping coefficients is given by Ishihara and Wang (2019). Except for the damping coefficients of these dashpots, the damping matrices in Eq. (1) can be obtained by the Rayleigh damping model. 115

Input acceleration response spectrum
The design acceleration response spectrum is generally defined at the bedrock condition, from which artificial ground motions are generated by considering a given phase property. The input acceleration time histories at the footing bottom ̈0 used in THA can be then obtained from these artificial ground motions by SHAKE analysis (Schnabel et al., 1972) using the one-dimensional site transfer function. Meanwhile, the input acceleration response spectrum defined at the footing bottom is 120 required in RSM. Several formulae of the acceleration response spectrum are presented, such as Eurocode (2006) where and are the characteristic period and damping ratio; 0 is the peak ground acceleration at the bedrock condition; is the soil amplification factor; is the damping correction factor; 0 is the acceleration response magnification ratio for 125 the region where the acceleration response becomes constant; , , , 1 , and 2 indicate coefficients representing the shape of the response spectrum. Parameters used for defining the input response spectrum in this study are listed in Table 1.
The peak ground acceleration 0 is chosen so that its return period is 475 years as recommended in IEC61400-1 (2019). The design response spectrum at the bedrock condition is typically defined with the damping ratio = 0.05, and the 130 soil amplification factor and damping correction factor will both be one. In the input response spectrum, on the other hand, and should be carefully defined according to the investigated soil profile and structure. In this study, the soil amplification factor is obtained using the one-dimensional site transfer function based on the response spectrum-based method by Okano and Sako (2013). In addition, the damping correction factor employed in this study is expressed as: where is the quantile value. It should be noted that this equation will be one with the damping ratio ζ = 0.05. Even though 135 the input acceleration time histories are generated from the same design response spectrum, their response spectra largely vary due to differences in phase properties especially in the case with the low damping ratio. Therefore, it is essential to quantify the uncertainty in the response spectra for estimating seismic loadings by RSM with the same reliability level as obtained by THA of several input acceleration time histories. The damping correction factor in Equation (3) is capable to estimate seismic loadings according to the reliability level by changing the quantile value . More detailed information about 140 the damping correction factor and quantile value is given by Kitahara and Ishihara (2020).

Augmented complex superposition RSM
The equation of motion of the seismic SSI model in Eq. (1) can be converted to a first order matrix equation as (Foss, 1958): The complex eigenvalue problem of this first order matrix equation is written as: On the other hand, the first order matrix equation in Eq. (4) can be also decoupled into the following single DOF equations as: where is the displacement response of the single DOF system. By the superposition of the solutions of Eq. (7), the relative displacement in Eq. (4) can be written as: Maximum relative displacements can be obtained based on the complex complete quadratic combination rule as (Gao et al., 2019): where is the highest mode considered in the calculation; and are the relative displacement and relative velocity response spectra of the th mode, respectively; , , and are the displacement-displacement, velocity-displacement, 160 and velocity-velocity correlation coefficients between the th and th modes, respectively. The displacement and velocity response spectra of the th mode can be approximately obtained from the given acceleration response spectrum by using the so-called pseudo spectrum transformation as: where ( , ) is the acceleration response spectrum in Eq.
(2) with the th natural period and modal damping ratio.
Moreover, the correlation coefficients , , and are expressed as: 165 where = ⁄ is the natural frequency ratio of the th to th modes.
Maximum seismic loadings at the th node of the seismic SSI model can be finally obtained based on Eq. (9) as: where 0.1 is the maximum value of the modal damping ratio based on an engineering judgement. This formula will be validated by case studies considering different types of foundations and different tower geometries in Section 3.2. 185

Structural damping of steel towers and additional loadings
Except for the damping coefficients of the dashpots, and , the damping matric in Eq. (4) is obtained by the Rayleigh damping model based on the first and second natural frequencies and modal damping ratios. In thus study, the first and second modal damping ratios are assumed to coincide with the structural damping ratio of steel towers, since these two modes correspond to the sway motion of the tower. The structural damping ratio of steel towers depends on their sizes and 190 an empirical formula was proposed by Oh and Ishihara (2018) as: where is the structural damping ratio, and 0.2 % is its minimum value for unlined welded steel stacks as shown in ISO 4354 (2009); is the characteristic period of the fixed foundation model of wind turbine support structures. In this study, the structural damping ratio is obtained by Eq. (16) and is employed as the first and second modal damping ratios.
Moreover, by the simplification of RNA in the seismic SSI model, this model cannot accurately evaluate the bending 195 moment at the hub height due to neglect of the mass moment of inertia of RNA. Hence, an additional loading by the angular acceleration at the hub height is proposed and this loading at the th node of the seismic SSI model is expressed as: where is the mass moment of inertia of RNA; ̈ is the angular acceleration at the hub height; ̈1 is the maximum acceleration of the first mode at the th node; is the correction factor. The maximum acceleration ̈1 can be obtained as: This additional loading will be validated by comparison with the bending moment at the hub height obtained by THA of the 200 full finite element (FE) model of wind turbine support structures, including the detailed configuration of the rotor and nacelle as shown in Section 3.2. In addition, an additional loading by the P-Δ effect is also proposed in this study as: , for = 1, 2, ⋯ , n − 1 , where is the gravitational acceleration; is the maximum displacement at the th node obtained by Eq. (9).

Numerical verification and discussion
The proposed method is demonstrated upon numerical examples accounting for different types of foundations and different 205 tower geometries. A typical 2-MW wind turbine is first investigated with different types of foundations, i.e., the gravity and piled foundations. Next, the rated power is varied from 1-MW to 3-MW for the piled foundation supported wind turbine. The outlines of these wind turbine support structures are summarized in Section 3.1. The description of two types of soil profiles used in this study and the input acceleration response spectra for these soil profiles are also presented in this section. Seismic loadings on the towers and footings are then estimated by the augmented complex mode superposition RSM and compared 210 with the results by THA and the original complex mode superposition RSM in Section 3.2. Table 2 summarizes the outline of the 2-MW wind turbine and its support structures. The structural damping ratio is estimated by Eq. (16) and is used to obtain the Rayleigh damping model. The footing mass is about six times total masses of the tower, rotor, and nacelle. The embedded ratio of the footing is assumed to be 0.2. For the piled foundation case, totally 215 eight piles are embedded and extended to the engineering bedrock. The water depth is assumed to be under the pile bottom.  Table 3 shows the outline of the piled foundation supported wind turbines with the rated power of 1, 1.5, 2, 2.5, and 3-MW constructed based on Xu and Ishihara (2014) and their support structures. Note that, the profiles of the piles are assumed to be same as those shown in Fig.2 for all these wind turbines. As the soil profiles, typical stiff and soft soil profiles, namely Soil type I and II shown in AIJ (2006) are 220 considered in this study. The gravity foundation is used for Soil type I whilst the piled foundation is used for Soil type II.  Table 4 presents the description of one-dimensional layered soil models for these two soil profiles. Table 5 summarizes the stiffness constants and damping coefficients of the soil-foundation systems for these two soil profiles.      Fig. 3 shows the input acceleration response spectra at the footing bottom obtained for both soil profiles, together with the design acceleration response spectrum at the bedrock condition. In these response spectra, the damping ratio is assumed to be = 0.05, the soil amplification factor is obtained by Okano and Sako (2014), and the damping correction factor is 235 assumed to be = 1. Compared with the design response spectrum at the bedrock condition, the input response spectrum at https://doi.org/10.5194/wes-2021-81 Preprint. Discussion started: 3 December 2021 c Author(s) 2021. CC BY 4.0 License. the footing bottom is amplified in the short period range less than 0.5 s for Soil type I, whilst it is amplified in the long period range larger than 0.5 s for Soil type II. Note that, the response spectra shown in Fig. 3 are not directly employed in RSM, but instead, is obtained for each modal damping ratio by Eq. (3) and is multiplied with these response spectra in the RSM procedure. 240 Meteorological Agency), and the other 11 ground motions show random phase properties. Input acceleration time histories at the footing bottom ̈0 are then obtained from these artificial ground motions for both soil types by DYNEQ (Yoshida and Suetomi, 1995), which allows similar analysis as SHAKE (Schnabel et al., 1972), based on the equivalent linearization method. Note that, the shear strain is less than 1 % for both Soil type I and II, and therefore the equivalent linearization method is applicable for both cases. The acceleration response spectra of the input acceleration time histories are also given 250 in Fig. 3 for the case with the phase properties of the four observed earthquake records. It can be seen the response spectra of these input acceleration time histories are good agreement with the input acceleration response spectrum, while they vary due to differences in the phase properties especially for Soil type II. Therefore, seismic loadings on wind turbine support structures obtained by THA using these input acceleration time histories will also vary. No matter what kind of method is employed in the seismic analysis, such uncertainty in seismic loading estimation should be carefully addressed, and in the 255 proposed method, it is achieved by the quantile value in the damping correction factor.

Seismic loading estimation by augmented complex mode superposition RSM
Seismic loadings acting on the towers and footings are estimated by the proposed method. In this study, the quantile value in the damping correction factor is set to be = 0.5 to estimate mean profiles of maximum seismic loadings. The results are https://doi.org/10.5194/wes-2021-81 Preprint. Discussion started: 3 December 2021 c Author(s) 2021. CC BY 4.0 License. compared with mean values of the results obtained by THA using the 15 input acceleration time histories. It should be noted 260 that, while only the mean profiles are considered in this study, the proposed method is also capable to estimate the seismic loadings according to the desired reliability level by changing the quantile value, and more detailed information is given in Kitahara and Ishihara (2020). Table 6 shows the first five natural frequencies and modal damping ratios of the 2-MW wind turbine support structures obtained by the complex eigenvalue analysis. In addition, the corresponding modal participation functions , for = 265 1, ⋯ ,5, are obtained and illustrated in Fig. 4. Considering up to the fifth mode fulfils the criteria of Model Code for Concrete Chimneys (CICIND, 2011). It can be seen that a very large damping ratio larger than 40 % arises at the third mode in case with Soil type II, while all the modal damping ratios are less than 10 % in case with Soil type I. As can be seen in the modal participation function, the third mode in case with Soil type II corresponds to the sway motion of the footing, and the large damping ratio arisen at this mode is caused by soil radiational damping. In the proposed method, this modal damping ratio is 270 substituted with 10 % by Eq. (15) as provided in the parentheses in Table 6 to avoid the underestimation of the shear force acting on the footing.   On the other hand, Fig. 6 shows a comparison of mean values of the seismic loadings acting on the footings obtained by the proposed method as well as CRSM and THA. It can be seen that, except for the shear forces in case with Soil type II, all the methods provide similar results. As shown in Kitahara and Ishihara (2020) with Soil type I and third mode in case with Soil type II, corresponding to the sway motion of the footing, are dominant in the shear force on the footings. In case with Soil type I, the fourth modal damping ratio is less than 10 %, and in such case, it is found complex mode superposition RSM is capable to estimate seismic loadings similar to the THA results. In case with Soil type II, conversely, the third modal damping ratio is larger than 10 %, and CRSM significantly underestimates the shear 295 force. In the proposed method, this modal damping ratio is hence substituted with 10 % to avoid such underestimation, and the estimated shear force is good agreement with the THA result, demonstrating that the proposed upper limit of modal damping ratios of 10 % is a reasonable choice.  Table 7 summarizes prediction errors in the seismic loadings at the tower base and footing by the proposed method and CRSM compared with the results by THA. As can be seen, the prediction accuracy of the proposed method is quite well and prediction errors are less than 6 % for all cases regardless of the foundation type and soil profile, whilst CRSM significantly underestimates the shear force on the footing and the prediction error is larger than 30 %. support structures, including the detail configuration of the rotor and nacelle. It can be seen that those two profiles match 310 well and the proposed additional loadings are capable to accurately evaluate the bending moments at the hub height due to the mass moment of inertia of RNA. At the same time, the additional loadings by the P-Δ effect are also evaluated by Eq.
(18). Table 8 summarizes the estimated additional loadings at the tower base 2 and the ratios of these additional loadings to the bending moments at the tower base 2 | 2 | ⁄ . It can be seen 2 | 2 | ⁄ is less than 3 % for both soil profiles, thus the P-Δ effect can be ignored in the practical use. 315  The accuracy of the proposed method is then further verified considering different tower geometries. Fig. 8 illustrates a comparison of the shear forces and bending moments at the 1/2 height, tower base, and footing estimated by the proposed 320 method and the those mean values obtained by THA using the 15 input acceleration time histories. The results obtained by CRSM are also plotted in this figure. It can be seen that the predicted seismic loadings on the towers by both the proposed method and CRSM show favourable agreement with those by THA. It is because the first and second modal damping ratios are less than 10 % for all cases and the proposed method degrades into CRSM. The same conclusion can be obtained for the bending moments on the footings. The shear forces on the footings, on the other hand, are significantly underestimated by 325 CRSM because very large damping ratios arise at the mode corresponding to the sway motion of the footing, while the proposed method accurately estimates the shear forces, demonstrating the proposed upper limit of modal damping ratios of 10 % is applicable regardless of the tower geometry.

Conclusion
In this study, the seismic SSI of wind turbine support structures is investigated using RSM. The non-classically damped seismic SSI model of wind turbine support structures is constructed, and its seismic loadings are analytically derived from complex mode superposition RSM. Complex mode superposition RSM is herein augmented by introducing the upper limit 335 of modal damping ratios of 10 % to improve the prediction accuracy of the shear force acting on the footing. Moreover, the mass moment of inertia of RNA and P-Δ effect are also considered as additional loadings analytically derived from complex 1. Seismic loadings, i.e., the maximum shear force and bending moment, are analytically derived from complex mode superposition RSM. This method is based on the complex eigenmodes obtained by the complex eigenvalue analysis, thus it is applicable to the non-classically damped systems. 2. A highly damped mode arises at fundamental frequency of wind turbine support structures in case with soft soil profiles, 345 and CRSM significantly underestimates the shear force acting on the footing, which is important for designing piled foundations. Thus, CRSM is augmented by introducing the upper limit of modal damping ratios of 10 % to avoid such underestimation.
3. The proposed method is validated by comparison with time history analysis (THA) accounting for different types of foundations and different tower geometries, and seismic loadings on the towers and footings estimated by the proposed 350 method show favourable agreement with the THA results.