Seismic Analysis of Wind Turbines
Aitor Arrospide Sanz MSc Thesis Danmark Tekniske Universitet Department of Civil Engineering 2016 (DTU)
DTU Civil Engineering June 2016
Supervisors: Postdoc Evangelos Katsanos Prof. Christos T. Georgakis Assoc. Prof. Sebastian Thöns
Fragility analysis of a monopile based offshore wind turbine subjected to multi-hazard environment (wind, waves, currents and earthquakes) is conducted within the framework of this study. The loading scenario comprising multiple hazards is analysied in time domain through the implementation of an aero-servo-elastic code (HAWC2) where the earthquake is fully coupled to the overall system. A wide set of loading scenarios (80 strong ground motions, 16 mean wind speeds and 4 sea-sates) are studied which enables the assessment of the performance of the offshore wind turbine for several exposures. Derived data are used to develop probabilistic seismic demand models (PSDM) for the ultimate limit state based on yielding & local buckling of monopole, and for the serviceability limit state related to accelerations at hub height. These PSDM-s are derived based on the power law regression analysis of the structural response. By means of developed demand models, fragility analysis of the infrastructure is performed and hence, the probabilistic assessment of the reliability of the offshore wind turbine subjected simultaneously to wind, waves, currents and earthquakes is provided. As a result of this MSc thesis, the probability estimate derived from fragility analysis is found rather high (0.5-0.7), even for moderate seismic intensity (𝑆𝑎 ⋍ 0.5𝑔), what highlights the impact of the earthquake on the global reliability of this infrastructure.
II
To my girlfriend, whose support, patience and love never failed. To my parents, whose unconditional support made this possible. To my brother, who always made me grow.
Copenhagen, June 2016
III
Abstract
............................................................................................................................. II
Dedication ............................................................................................................................III Table of contents..................................................................................................................... IV List of figures ........................................................................................................................ VII List of tables ............................................................................................................................ IX Acknowledgements .................................................................................................................. X Summary of the MSc thesis ................................................................................................... XI Introduction ........................................................................................................ 1 1.1
Wind energy: a highly growing industry.................................................................... 1
1.2
Literature review on wind turbines & earthquake hazard .......................................... 5 State-of-the-art ..................................................................................................... 5 State-of-practice ................................................................................................... 9
1.3
Scope, motivation and methodology ........................................................................ 11 5MW OWT: Structural Modelling and Loads Definition ............................ 13
2.1
Description of the 5 MW Offshore wind turbine ..................................................... 14
2.2
Aero-servo-elastic code, HAWC2............................................................................ 15 Structural modelling of the wind turbine based on multi-body-dynamics theory.. ............................................................................................................................ 16 Eigenvalue analysis of the 5MW offshore wind turbine model ......................... 17
2.3
Loads: theory, definition and implementation in HAWC2 ...................................... 19 Environmental loads: wind and waves loads ..................................................... 19 Seismic loads: definition and modelling ............................................................ 27 Validation of HAWC2 for seismic analysis .................................................... 28
3.1
HAWC2 vs SAP2000: validation of seismic analysis .............................................. 28
3.2
Sensitivity analysis, initial conditions and time step for time-domain simulations . 32
IV
Response of the wind turbine subjected to earthquake excitation .............. 35 4.1
Wind turbines operational states .............................................................................. 35
4.2
Seismic analysis of the wind turbine ........................................................................ 36 Fragility analysis .............................................................................................. 41
5.1
Probability of exceeding the limit state .................................................................... 41
5.2
Demand - Capacity considerations ........................................................................... 43 Ultimate Limit State (ULS) ................................................................................ 44 Serviceability Limit State (ULS)........................................................................ 53 Summary of capacity estimate parameters ......................................................... 54
5.3
Definition of the Multi-hazard environment ............................................................ 57 Wind loads.......................................................................................................... 57 Wave loads ......................................................................................................... 58 Earthquake loads ................................................................................................ 58
5.4
Time-Domain analysis strategy ................................................................................ 64 Results ............................................................................................................... 65
6.1
Fragility curves for a reference case ........................................................................ 65 ULS, structural integrity ..................................................................................... 66 SLS, component level integrity .......................................................................... 71 ULS and SLS Probabilistic Seismic Demand Model summary ......................... 72
6.2
Fragility curves and surfaces for a varying multi-hazard environment ................... 74 Ultimate Limit State, fragility estimate .............................................................. 75 Serviceability Limit State, fragility estimate ...................................................... 79 Influence of the sea state for earthquake assessment of an OWT ...................... 81 Concluding remarks ......................................................................................... 83
7.1
Conclusions .............................................................................................................. 83
7.2
Recommendation for future research ....................................................................... 84
References ............................................................................................................................ 86 Suite of 40 pairs of motions ...................................................................... 92 Fragility surfaces and PSDM-s ................................................................ 94 B.1
Significant Wave Height, 𝑯𝒔 = 𝟎. 𝟓𝒎 .................................................................... 95 V
B.2
Significant Wave Height, 𝑯𝒔 = 𝟐. 𝟓𝒎 .................................................................... 98
B.3
Significant Wave Height, 𝑯𝒔 = 𝟒. 𝟓𝒎 .................................................................. 101
B.4
Significant Wave Height, 𝑯𝒔 = 𝟔. 𝟓𝒎 .................................................................. 104
VI
Figure 1-1: Electricity demand increase relative to 2002 (renewable vs any source), data from EIA [2] ................................................................................................................... 1 Figure 1-2: Global wind power installation from 2000 to 2015, data from GWEC [4]............. 2 Figure 1-3: Wind power installed in 2015 & share of total cumulative installed power, data from GWEC [4] .............................................................................................................. 3 Figure 1-4: Capacity of Chinese offshore wind market, data from GWEC [4] ......................... 5 Figure 2-1: Coord. systems for HAWC2 (left) [26] and 5MW OWT model studied (right) ... 13 Figure 3-1: Time history of strong ground motions and response spectra of Kobe and Manjil earthquakes .......................................................................................................... 31 Figure 3-2: Moment demand at base frame, HAWC2 vs. Sap2000, Manjil & Kobe earthquakes. (a) refers to mast and (b) to short column response. ........................................... 32 Figure 3-3: Wind ramp 𝑈𝑚, pitch angle 𝜃 & shaft rotational velocity 𝛺 ................................ 33 Figure 3-4: Time-step influence on HAWC2, relative error vs time cost ................................ 34 Figure 4-1: Time domain response analysis of a wind turbine subjected to an earthquake ..... 37 Figure 4-2: Tower-top deflection in X and Y & moment at the base of MP around X and Y (operating, shut-down and parked states) ............................................................ 38 Figure 4-3: Maximum and minimum response along the height of the wind turbine to a single seismic motion, operating (Op.), shut-down (Sh.d) & parked (Pa.) .................... 40 Figure 5-1: Elastic stress distribution in a tubular cross section .............................................. 45 Figure 5-2: Membrane stresses in shells, from [42] ................................................................. 46 Figure 5-3: Buckling capacity along the monopile (MP) height .............................................. 52 Figure 5-4: Linear Bifurcation Analysis (LBA) result for 5m high can and loading M=1Nmm ............................................................................................................................. 53 Figure 5-5: World map presenting the distribution of the 173 seismic events stored in the PEERNGA database (distribution computed based on the information available online in http://peer.berkeley.edu/peer_ground_motion_database) ............................... 59 Figure 5-6: Earthquake magnitude as function of epicentral distance and peak ground acceleration of the strong ground motion records available in PEER-NGA database (distribution computed based on the information available online in http://peer.berkeley.edu/peer_ground_motion_database) ................................... 60 Figure 5-7: Natural PGA magnitude and direction (around Z axis of the tower) .................... 61 VII
Figure 5-8: Elastic response spectra for operating 𝜁 = 5% and parked 𝜁 = 1% wind turbine [54]& [55], x & y directions ................................................................................ 63 Figure 6-1: Shear failure mode, PSDM .................................................................................... 67 Figure 6-2: Moment failure mode, PSDM ............................................................................... 67 Figure 6-3: vonMises criteria’s failure mode, PSDM .............................................................. 68 Figure 6-4: Buckling failure mode, PSDM (critical height identification) .............................. 69 Figure 6-5: Fragility curves for structural reliability, CDF (ULS) .......................................... 70 Figure 6-6: Acceleration at hub height, PSDM ........................................................................ 71 Figure 6-7: Fragility curves for component level reliability, CDF (SLS) ................................ 72 Figure 6-8: Stress demand in the cross section for wind speeds ranging from 4 to 47 m/s and Hs=6.5m ............................................................................................................... 75 Figure 6-9: Cumulative Distribution for the probability of exceeding the yielding limit for wind speeds ranging from 4 to 47 m/s and Hs=6.5m ................................................... 76 Figure 6-10: Buckling demand in terms of Utilization Ratio at a height of 4 m from the support for wind speeds ranging from 4 to 47 m/s and Hs=6.5m ..................................... 77 Figure 6-11: Cumulative Distribution for the probability of exceeding the buckling limit for wind speeds ranging from 4 to 47 m/s and Hs=6.5m .......................................... 77 Figure 6-12: Comparison of fragility curves for yielding limit state (LS1) and buckling limit state (LS3) according to the cut-in, rated and cut-out wind speeds..................... 78 Figure 6-13: Acceleration demand at the hub for wind speeds ranging from 4 to 47 m/s and Hs=6.5m ............................................................................................................... 79 Figure 6-14: Cumulative Distribution for the probability of exceeding three significant acceleration thresholds at hub height for 4 to 47 m/s and Hs=6.5m .................... 80 Figure 6-15: Comparison of fragility curves for acceleration limits at hub height (SLS) according to the cut-in, rated and cut-out wind speeds (SLS1 – Elec. generators, SLS2 – Inv.&trans. and SLS3 – elec. control.) ................................................... 80 Figure 6-16: Influence of the sea states in the fragility estimate for ULS in yielding ............. 81 Figure 6-17: Influence of sea states in the fragility estimate for ULS in buckling .................. 81 Figure 6-18: Influence of sea states on the fragility estimate for SLS (ele. gen. 𝑆𝑐 = 0.9𝑔) .. 82
VIII
Table 1-1: Global installed wind power capacity in 2015 [GW], continental basis .................. 4 Table 2-1: NREL 5-MW specifications [11]............................................................................ 14 Table 2-2: Structural properties of the monopile and tower [11] ............................................. 15 Table 2-3: Modal natural frequencies and modal damping ratio of the MBD system ............. 18 Table 2-4: Rayleigh damping coefficients for the main bodies of the multi-body-system ...... 19 Table 3-1: Geometrical and material properties of the validation models ............................... 29 Table 3-2: Natural Eigen frequencies of validation models, [Hz] ........................................... 30 Table 3-3: Characteristics of Kobe and Manjil earthquakes applied for validation purposes . 31 Table 4-1: Environmental load specification for reference case scenarios .............................. 36 Table 4-2: Five earthquakes applied for seismic response study of a wind turbine ................. 36 Table 5-1: Nominal and mean yield strength for S235 and S355 ............................................ 55 Table 5-2: Median yield strength and logarithmic standard deviation for S235 and S355 ...... 56 Table 5-3: Capacity estimate for structural and component reliability .................................... 56 Table 5-4: Mean wind speeds simulated .................................................................................. 58 Table 6-1: Correlation comparison 𝐼𝑀𝑃𝐺𝐴 𝑣𝑠. 𝑆𝑎 − 𝑆𝑑 by means of the Coefficient of determination, 𝑅2 ................................................................................................ 66 Table 6-2: Summary of PSDM results, power model definition for operation (OP), parked (PK) and operating & shutdown (SH) .......................................................................... 73
IX
Firstly, I would like to thank my supervisor, Dr. Evangelos Katsanos, for his infinite motivation and unconditional support during this MSc thesis at Denmark Tekniske Universitet. His dedication and guidance have challenged me to progress and evolve in many aspects. For his support, contribution and generosity I would like to show my profound gratitude. In addition, I would also like to express my gratitude to my supervisors Prof. Christos T. Gerogakis and to Assoc. Prof. Sebastian Thöns for their time and experience which have enriched my study. In the same way, I would like thank Eilif Svensson for participating together with my three supervisors on the evaluation committee as an external examiner. Apart from that, I would like to thank sincerely my colleague Albert Meseguer Urbán, for our long discussions and his expertise within the wind energy field what has improved the quality of this study from different perspectives. Finally, I would like to appreciate the input from Ulrik Muurholm Hanse and Fernando Vidal Puche through our technical discussions during my stay in COWI A/S.
X
This MSc thesis presents a numerical study of the response of a 5MW offshore wind turbine subjected to a fully coupled loading scenario composed of environmental loads (i.e., wind, waves, currents) and the earthquake hazard. The analysis is conducted through time-domain simulations performed by means of the aero-servo-elastic code HAWC2. Furthermore, the wind turbine model studied along these lines is the NREL 5MW monopile based offshore wind turbine which was implemented for the Offshore Code Comparison Collaboration (OC3, 2010). The increasing demand for accurate estimation of the reliability of wind energy infrastructure subjected to the seismic hazard is derived from the fact that wind energy industry is spreading in earthquake prone areas as a consequence of its global growth. Hence, in response to this demand, throughout this study probabilistic seismic demand models are developed based on the data obtained from time domain analysis. With that purpose, HAWC2 code is extended to incorporate strong ground motions in order to enable the aero-servo-elastic simulation of wind turbines subjected to seismic motions. Moreover, due to the novelty of this procedure, a validation process is carried out, where the seismic response of two distinct models is compared between HAWC2 and Sap2000, the latest being an already recognized and widely used software for earthquake engineering. Subsequently, the study of critical spots regarding the demand & capacity is conducted based on maximum response studies along the height of the structure. As a consequence, two limit states are selected. The first with regard the Ultimate Limit State for yielding and/or local buckling at the base of the monopile and the second related to the Serviceability Limit State derived from the accelerations at hub height. In addition, a set of 40 earthquakes of different characteristics is selected from PEER-NGA database with the purpose to build a generic suite of earthquakes. Moreover, they are scaled up once, with the aim to produce the most homogeneous possible sample. Similarly, in line with the aim to provide a generic study which could provide relevant data for a broad set of scenarios, an extensive combination of environmental loads (16 mean wind speeds and 4 sea-states) were analysed, which leads to 5120 time domain simulations. By means of this large sample of load case scenarios analysed, fragility analysis is perfomed in order to quantify the relevance of the multi-hazard environment (wind, waves, currents and XI
earthquakes) on the integrity of the wind turbine based on a probabilistic approach. This procedure is applied to determine the fragility of the wind turbine for ULS1 (yielding), ULS3 (local buckling) and SLS (accel. hub height). Consequently, throughout this study, the earthquake is identified as the governing hazard on the response of the wind turbine when subjected to the multi-hazard environment. Moreover, the influence of the mean wind speed and the significant wave height is found to be limited (specially for the case of the wave height). However, the majority of the simulations showed maximum demands for the operating case which corroborates the need to consider the fully coupled loading scenario and the aero-servo-elastic analysis in order to simulate the maximum response of the wind turbine accurately. Additionally, the probability estimates derived from the fragility analysis, even for moderate seismic intensity (𝑆𝑎 ⋍ 0.5𝑔), have been found quite high (0.5-0.7) highlighting the impact of the earthquake on the reliability of this infrastructure.
XII
Modern societies demand incessantly larger quantities of electricity and this motivates the inflation of the electric power consumption worldwide, which is substantially more intensive in the developing countries. However, the environmental concern and ecological awareness of the societies in a worldwide level keep increasing and hence the demand for green energy production systems is becoming of higher importance. The latter has already promoted a considerable increase of global investments in the renewable energy sources like the solar and wind energy. Therefore, due to the expansion of different renewable energy industries, it can be considered that currently green energy sources have become a real alternative.
The relevant role that the renewable energy sector has in modern societies can be identified through Figure 1-1, where the global trend of electricity demand is compared to the growth that electric power generation by means of renewable sources has experienced [1].
Figure 1-1: Electricity demand increase relative to 2002 (renewable vs any source), data from EIA [2]
1
It should be noted that for both cases (i.e., electricity demand produced only by renewable sources and using all resources available), the reference year is considered to be 2002 and hence, the trends are represented in relative perspective. Moreover, in line with the tendency of renewable energy production found in Figure 1-1, in 2015 the maximum renewable power capacity installed in a year was reached exceeding by 30% the installed capacity in 2014 [3]. A significant contribution to this additional power capacity was provided by almost 64GW of wind power and 57GW of solar PV energy. As for the investment rate, mainly due to the upswing observed in China, the African continent, the United States, Latin America and India, a substantial increase was experienced, up to $328.9bn in 2015. The latter represents 3% of growth with respect to the previous record obtained in 2011. In several occasions investments in green energy have been linked to the current market price of fossil fuels. Despite that, the last 18 months till December 2015 have shown the opposite. The price of the crude fell 68%, coal dispatched in the NW Europe plunged 35% and at the same time natural gas in US dropped 48% [3]. And in spite of that, green energy production and its correspondent investment do not cease surging. In the US, the Energy Information Administration [2] predicted that from 2012 to 2040 an increase on renewable electricity production of 69% is going to happen. Comprising a 140% grow in non-hydropower sources, which mainly involves solar and wind energy. Therefore, and based on the trend shown in Figure 1-2 the expansion of wind energy industry seems unavoidable.
Figure 1-2: Global wind power installation from 2000 to 2015, data from GWEC [4]
2
Nowadays, Asia has become the biggest market for wind energy and especially China is leading this expansion. Between 2011 and 2015, China has increased its electricity production by means of wind power in more than 230%, stepping up from 62 GW to 145.36 GWEC [4]. Figure 1-3 reflects the predominant position of China in the wind energy market. The development of the Indian wind energy market is also increasing and this is highlighted from the rate of growth of wind power infrastructures. In addition, South America is also experiencing and intensive growth of the wind energy industry and more specifically, Brazil, being the biggest player in this area, is currently among the first eleven countries in terms of wind energy production. Figure 1-3 depicts the distribution of the installed wind energy in 2015 while the cumulative quantity is also presented. It is clearly visible the leading position of China which currently has the largest capacity of wind power installation and additionally keeps increasing its capacity faster than any other country.
Figure 1-3: Wind power installed in 2015 & share of total cumulative installed power, data from GWEC [4] Hence, based on the current situation and the anticipated future trends regarding the expansion of wind energy industry, it is profound that this highly promising green energy sector is increasing constantly and new installations are spreading in several areas in the world, where numerous hazardous sources may impose intensive threats to the reliability of these energy infrastructures during its lifetime. Particularly, a significant amount of onshore or offshore wind turbines have been already installed in earthquake prone areas including several Chinese areas as well as seismically active regions in South and North America, India, Korea and Southern 3
Europe (see Table 1-1). Thus, the earthquake-resistant design of the wind turbines should be substantially refined. Otherwise, the seismic risk involved in the wind turbines, along with the conventional exposures, like the wind and the waves (in case of offshore installations) may significantly threat the structural integrity and reliability of these critical infrastructures. In addition, due to the increasing number of turbines installed in wind farms and their inherent similarities, the threat derived from seismic loads results in a significant risk which is essential to be included in the reliability assessment of these infrastructures. So far, for the structural analysis of wind turbines and their design, it has been widely assumed that the governing load is produced by the wind inflow. However, some studies (see chapter 1.2.1), have already shown that earthquake-imposed dynamic loading may govern the design. Table 1-1: Global installed wind power capacity in 2015 [GW], continental basis
Africa & Mid. East
Asia
Europe
Lat. Am. &
North
Pacific
Caribbean
America
regions
World
End 2014
2.54
141.97
134.25
8.57
77.93
4.44
369.71
New 2015
0.95
33.86
13.8
3.65
10.82
0.38
63.47
Growth[%]
37.4
23.85
10.28
42.59
13.88
8.56
17.17
On the other hand, offshore wind industry is undergoing a steep increase. In 2015 3.4 GW of new capacity were installed which led to reach the cumulative power capacity of 12 GW. The foremost part was held in Europe where Germany stalled most of it. However, some of the relevant new markets such as China or India have started to explore the offshore area. As an example of it, Figure 1-4 shows the development of the wind power installed offshore in China, from 2007 which was the initiation. Thus, the future trends and forecasting of the wind energy industry tendencies point towards the expansion into the offshore environment and especially in Asia and South America both containing earthquake prone regions. That is the reason why this MSc thesis aims to study the influence of seismic motions in horizontal axis Offshore Wind Turbines (OWT), based on monopile foundations.
4
Figure 1-4: Capacity of Chinese offshore wind market, data from GWEC [4]
This chapter includes the most relevant and recent studies focusing on the wind turbines and their exposure to the earthquake hazard. The vast majority of the existing literature includes numerical and/or analytical studies, while the experimental ones are too rare mainly because of the obstacles to create a realistic set up for a wind turbine that will be subjected to seismic motions. Since this study is a numerical one, the pertinent experimental studies are deliberately excluded by this following review.
As mentioned above, wind turbines have been widely extended to many areas around the globe in the last decades. This spreading has led to broaden the quantity of different hazards that may threaten the infrastructure. Moreover, the size of wind farms as well as their relevance in the energy industry are increasing what represents the importance of these infrastructures. Hence, related to the size and the potential impact of wind energy on the global energy market, the reliability assessment becomes a major player due to the great consequences that might be derived from the new multi-hazard scenario for wind turbines where the earthquake has to be involved.
5
The need for the consideration of the earthquake hazard for risk assessment of the wind turbines is not a new topic. Nevertheless, the implementation of complex aero-elastic codes (see chapter 2.2) facilitated by the increasing computational power through the last decades, has enabled a consistent treatment of the uncertainty related to the environmental loads as well as the control system of the turbines. By means of this brief literature review, some of the relevant studies related to wind turbines subjected to seismic hazards are intended to be presented. Along these lines, two main different approaches identified in the existing literature are shown. Studies based on a deterministic approach for loads definition and studies, where more refined and advanced probabilistic approaches have been adopted.
Within this subchapter, some of the studies that have assessed the impact of earthquakes on wind turbines based on a deterministic approach are referred. Time Domain Analysis (TDA) as well as Frequency Domain Analysis (FDA) have been alternatively applied in order to study the influence of the seismic hazard on wind turbines. Additionally, the accuracy of the model applied for analysis purposes is another characteristic of these studies that is considered along this lines. Full FEM models including blades and nacelle features have been compared to simplified beam models where the hub and blades were modelled as a lumped mass at the tower top. Bazeos et al. [5], performed and analysis where two refinement levels of the 37 m high structural model were applied and the earthquake demand was captured both derived from FDA and TDA analysis methods. Frequency domain analysis was performed applying the complex 3D FEM model while the simple beam model with lumped masses was utilized for the time domain analysis. Both cases experienced very similar results, which concluded that for low seismicity areas the loads derived applying pseudo-static aerodynamic loads were greater than those related to the earthquake. Furthermore, a similar study was conducted by Ritschel et al. [6] where the differences between the FDA and TDA showed the conservativeness of the frequency domain analysis. Despite that, both procedures reflected the governance of the wind with regard to earthquake loading acting over a 60 m high wind turbine. Related to the conservativeness of the FDA, Witcher [7] included the influence of the operational state of the wind turbine in the overall response. From that perspective, the usual 5% damping ratio 6
implemented by the FDA was found too far conservative for the parked turbine case and only acceptable for an operational state. Hence, the need to consider the aerodynamic damping was shown. In reality, wind turbines are interactive devices which are designed intended to achieve the maximum power production despite the turbulence intensity or the mean wind speed. In order to accomplish that, the control system is constantly updating the pitch angle of the blades based on the detected rotational speed of the shaft or the torque generated. By means of acquiring the most constant torque or shaft angular velocity, the maximum rated power production is ensured. Additionally, the control system can even generate a shutdown of the turbine if excessive vibrations are detected at hub height, leading to a totally different scenario in terms of aerodynamic properties of the device. Due to the aforementioned characteristics, it can be assumed that the aerodynamic damping of a wind turbine is not a fixed value, which considering its relevance in the response, it is evident the superiority of the TDA over the FDA when the response of a wind turbine to seismic events is aimed to be investigated. Furthermore, the characteristics of the earthquake might be influential in the analysis as well, and they can only be represented by time domain simulations. In line with the time domain analysis, Wang & Zhang [8] assessed a 70m high wind turbine of 1.65 MW modelled in both ways simplified (lumped masses) and complex (full FEM) forms. These structures were subjected to a single ground motion of 0.51 g of peak ground acceleration (PGA) and even some divergence was identified in the response, both models assured the governance of the seismic motion in the structural demand. Similarly, some other studies i.e. Song et al. [9] or He & Li [10] confirmed the research performed by Wang & Zhang [8] where the relevance of the earthquake for the reliability of a wind turbine was found necessary. Nonetheless, aforementioned studies that applied TDA-s did not accomplish the full coupling of the environmental loads for the time domain analysis of the wind turbines and moreover, most of them applied a narrow sample of seismic motions. Aero-servo-elastic codes (see OC3 study [11]) enabled the implementation of the environmental loads in a precise manner, where the turbine and the controller were active during the analysis. By means of these codes, the interaction between the structural dynamics, the control mechanism and the environmental loads was possible to be captures and analysed in a fully coupled simulation. This method was partially applied by Asareh & Volz [12] who obtained the aerodynamic loads through an aeroelastic code and combined them to the demands derived from 22 different far-fault earthquake motion pairs. They concluded that the rated wind 7
speed produced the highest demands along the height of the support structure. Moreover, in agreement with the last reference, Alati et al. [13] carried out a thorough study of the impact of the earthquake on two wind turbines (jacket and tripod), considering three operational states. In this last study, a fully coupled scenario where the environmental loads (wind and waves) were combined to the earthquake hazard was accomplished through an aeroelastic code. Within that study [13], up to 49 pairs of seismic excitations were applied and evaluated for an offshore environment, which resulted in a confirmation of the fact that the dominating design load was derived from the seismic excitation. In addition to Alati et al. [13] several other studies were carried out by means of the most recent aeroelastic codes where the fully coupled scenario was modelled precisely. That was the case for Prowell [14], who extended the analysis to the comparison of different turbine models. The whole structure was assessed and the relevance of the earthquake was evidenced. 99 motions were selected with the aim to cover different features and that way overcome the inherent uncertainty within the earthquake motions. Apart from modal codes, more complex multi-body models were also applied for the aerolestic simulations as can be found through Jin et al. [15] where the power production was identified to be threatened even for no inflicted structural detriment.
Within the probabilistic framework for the study of the wind turbines excited by an earthquake much less research cases can be detected. Some of those cases are commented in the following. Kiyomiya et al. [16] established a probabilistic method where earthquakes and wind loads were included within the framework of wind turbines reliability assessment. By means of the Joint Proabability Density function derived from Weibull distribution for both (wind and seismic motions), the probability for the joined occurrence of specific intensity measures was determined. Additionally, the demand was obtained through pseudo-static application of wind loads and TDA by means of a single seismic event. This study concluded that the wind was the predominant load for the wind turbine. Nevertheless, the size of the applied turbine model (early-stage) and the reduced amount of motions evaluated makes the relevance of the analysis of limited extent for an accurate prediction of the structural reliability. Further investigation was conducted by Nuta et al. [17] who performed their analysis based on a single-hazard scenario where the earthquake was solely applied in time domain to a simplified 8
FE model which neglected the geometry of the blades and hub. By virtue of an incremental dynamic analysis procedure (IDA) that was applied stepwise by increasing the intensity of the earthquake, they generated a sample of data which was implemented to derive probabilistic estimations for different Damage Stages (DS) for two specific geographical locations. Certainly, the probabilities for an earthquake prone area were found to be of considerable quantity. On the other hand, Mensah et al. [18] run an analysis where wind and seismic waves were coupled through an aeroelastic code. In this study, they simulated a contemporary wind turbine (5MW) exposed to a number of real ground motions and different wind scenarios. Through the research conducted by Mensah et al. [18], it was concluded that the seismic demand governed significantly over the wind loads, hence the need to consider the earthquake hazard for wind turbines was corroborated. Finally, the extended research conducted by Mardfekri & Gardoni [19] is worth of mentioning. Along their research they built a highly developed probabilistic model where the multi-hazard environment of an offshore scenario was included and the reliability of the wind turbine derived. The analysis included variations on the structural model covering a wide range of different structural characteristics as well as fragility analysis of the operational wind turbine. In addition, by joining those fragility curves with the annual occurrence probability of the studied hazards intensities, it was concluded that the actual probability of failure is far from being neglected.
Wind turbines have barely been installed in earthquake prone areas during the last decades. Yet the growth of the industry has led to the spreading of this power production infrastructure around the globe, reaching seismically active regions. That is why, it becomes of interest to refer to the State-of-Practice. In the course of this subchapter, the prescriptions related to earthquakes provided by the principal codes and standards applied for design of wind turbines in the Europe and United States are intended to be brought to light. Firstly, the most widely known IEC 61400-1 [20] standard is considered. This standard is of common application in Europe as well as in US. It includes the requirement for seismic design of the wind turbines if the location is classified as earthquake prone area. The design procedure specified in this code is based on the FDA, where the modal response spectra is the base of the 9
design process (local codes referred). Through this method, the aerodynamic loads are afterwards superposed to those derived from the modal analysis, where the earthquake loads are recommended to be combined with normal operational wind conditions. Additionally, it offers the option to run TDA when the number of simulations is broad enough. Despite all that, there is no minimum requirement for earthquake design criteria, since there is assumed that seismic events are design driving for very limited regions on the globe. Apart from that, there is another major standard applied for design of wind turbines foundation structures in Europe, DNVGL-ST-0126 [21]. Since the union of Det Norske Veritas & Germanischer Lloyd in September 2013, harmonization of two especially relevant guidelines has given rise to this code. Along this last standard, several mentions are done to the seismic hazard regarding the need to determine the geological and soil characteristics of the site in order to predict the influence that an earthquake might have. Soil structure interaction is also mentioned as a necessary aspect to be analysed. Nevertheless, as far as the design procedure is concerned, it does not provide details in order to proceed with the design of an earthquake resistant wind turbine, it refers to EC 1998-1 [22] and ISO 19901-2 [23] insted. However, the first reference is mainly oriented to building structures and the second is defined to be applied for offshore structures within oil and gas industry. Therefore, due to the complexity of wind turbine’s response as a consequence of the interaction of the control system, aerodynamic & hydrodynamic loads (offshore) as well as the structural flexibility, it makes the structure unsuitable for frequency domain analysis presented by the previous EC 1998-1 [22] or ISO 19901-2 [23]. Additionally, frequency domain analysis is dependent on specific values such as the behaviour factor (𝑞), importance factor (𝛾𝐼 ) or damping ratio (𝜁) which are critical for response spectra definition purposes. Then, since the wind turbines are not within the prescribed types of structures, the determination of these parameters turns out to be unclear. The suitability of the FDA is found to be even worse when the damping is considered, which is a function of the operational state, wind velocity and the quality of the control system. That means that even considering two cases of 0.5-1.0% of critical damping for the parked and 5% for the operating case, the accuracy might be questioned. As for the United States of America, the lack of specific codes for wind turbines besides the IEC 61400-1 [20] is considerable. Owing to that, the design of a wind turbine has to be based on European codes or codes prescribed for building structures. In this context, 2006 IBC [24] code represents the basis for the main part of the guidelines. Despite that, it does not include
10
specifically prescriptions to assess earthquake resistant structures and thus it refers to ASCE 710 [25] for that purpose. Nevertheless, this second standard requires different parameters related to the structural properties and usage or application in order to proceed with the earthquake design. Some of those parameters are the occupancy category, site characteristics (soil type), response modification factor (𝑅) which is determined based on the ductility level, damping and over-strength features. Therefore, the application of this codes for wind turbines is not a straight forward approach since they are not included as usual structures considered by these codes, similarly to the limitations found in EC 1998-1 [22]. Hence, aforementioned provides sustenance to the fact that the current growing wind energy industry demands for more accurate and all-inclusive time domain based design prescriptions and guidance since regions located in earthquake prone areas are investing in wind energy. So that the methodology to assess the demand derived on an Offshore Wind Turbine (OWT) subjected simultaneously to wind, waves and earthquake ground motions is becoming an actual need.
The state-of-industry review shown in chapter 1.1 has demonstrated that wind energy industry is growing and thus reaching earthquake prone areas. That is why the assessment of the impact that an earthquake has on this green energy harnessing infrastructure is becoming of high interest. Consequently, as the literature review presented above in chapter 1.2 reveals, some studies have been conducted with the aim to quantify or predict the damage that an earthquake could cause on a wind turbine. Different approaches have been studied for this purpose, however, a limited amount of those studies include the earthquake in a fully coupled multihazard environment and even less are based within a probabilistic framework. Furthermore, along the state-of-practice reciew it is identified that codes and standards implemented for wind turbine design purposes prescribe the need for seismic load consideration if the project is designed for seismically active areas. Nevertheless, accurate standardized procedures that avoid too far conservative approaches (pseudo-static superposition of loads) are not defined. Hence, due to all the aforementioned, this MSc thesis intends to assess the demands exerted on an OWT as a consequence of a fully coupled action of environmental loads (wind, waves and currents) and earthquake hazards within a probabilistic framework based on the fragility assessment. In order to provide a generic and widely relevant study, a common monopile based 11
offshore wind turbine of 5MW is selected to conduct this analysis. In addition, the wind turbine is subjected to the aforementioned actions (wind, waves, currents and earthquakes) through time domain analysis performed by means of an aero-servo-elastic code. The aeroelastic code implemented along this study is HAWC2, whose application for seismic assessment elsewhere has been found to be inexistent. Through this code (HAWC2), the simulation of wind turbines for and scenario where all coupled loads as well as the structural elasticity and the advanced control system of contemporary wind turbines is combined in time domain has been made possible. Consequently, by implementation of HAWC2, wind turbine’s response is aimed to be determined, which accordingly data derived from response analysis is applied in order to derive Probabilistic Seismic Demand Models (PSDM). This leads to the probability estimate of exceeding a deliberately specified demand threshold with regard to the Intensity Measure of the earthquake (IM). Herein, all required steps to build up this Master thesis are described in detail.
12
In this chapter of the MSc thesis report, the horizontal axis offshore wind turbine (HAOWT), adopted herein for research purposes, is comprehensively described, while the basis of the aeroelasticity principles is presented along with the relevant aero-elastic code, HAWC2, utilized in the current study to model the wind turbine. Moreover, this chapter includes the definition of the environmental loads, i.e., wind and waves, as well as the earthquake loads and the related seismic analysis, conducted within the HAWC2 framework. Various modelling techniques are also described, since they are necessary for the time-domain simulations.
Figure 2-1: Coord. systems for HAWC2 (left) [26] and 5MW OWT model studied (right)
13
The wind turbine utilized for this analysis is the NREL 5-MW (National Renewable Energy Laboratory) wind turbine, which is the reference wind turbine model for several studies and related research work. The main characteristics of this turbine are specified in Table 2-1 while more details about the particular wind turbine can be readily found elsewhere (e.g., [11]). Table 2-1: NREL 5-MW specifications [11] Rating
5 MW
Rotor orientation, configuration
Upwind, 3 blades (61.5 m long)
Control
Variable speed, collective pitch
Drivetrain Cut-in, rated and cut-out wind speeds Cut-in and rated rotor speed Rated tip speed
High speed, multiple-stage gearbox 3 m/s, 11.4 m/s, 25 m/s 6.9 rpm, 12.1 rpm (≈ 0.72 rad/s, 1.27 rad/s) 80 m/s
Rotor mass
110,000 kg
Nacelle mass
240,000 kg
Tower mass
347,500 kg
Co. location of overall centre of mass (CM)
(0.0 m, -0.2m, 64.0 m)
As mentioned above, the specific wind turbine is based on a monopile (MP) foundation of 30 meters installed in 20 meters water depth. In addition, the tower installed on top of the monopile is 78 meters high and comprises a linearly tapering cross section contrarily to the constant section of the MP. The structural properties of the monopile as well as the tower are described in Table 2-2. As far as the wind turbine implemented for this study is concerned, it is fully described in by Jonkman et al. [27], where all details about the nacelle, shaft, hub and the blades can be found. 14
Table 2-2: Structural properties of the monopile and tower [11] Monopile Foundation
Tower (base/top)
6.0 m
6.0 m / 3.87 m
60.0 mm
27.0 mm / 19.0 mm
Location of base
20.0 m below the mudline
10.0 m above the mudline
Location of top
10.0 m above the mudline
87.60 m above the sea level
Diameter (ext.) Thickness
Young’s modulus
210,000 MPa
Shear modulus
80,800 MPa
Eff. density of steel Overall support mass Centre of gravity (C.G.)
8,500 kg/m3 * (Monopile plus tower) 522,617 kg 37.172 m (from mudline in vertical centreline)
* In order to account for bolts, welds, paint and not considered flanges
The analysis of the response of the wind turbine, which is subjected to the multi-hazard environment consisting of the wind and wave loads as well the earthquake excitation, will be performed in the time domain. Both the modelling of the wind turbine and the consequent time history analysis will be performed using the aero-servo-elastic code called HAWC2 (Horizontal Axis Wind turbine Code, 2nd generation). HAWC2 has been developed by DTU Wind Energy Department and it required intensive effort from its first development until the current and more advanced version.
HAWC2 code is based on the aeroelasticity theory and was built to conduct time domain analysis of the response of wind turbines. It was mainly developed from 2003 to 2007 within the framework of the Aeroelastic Design Reasearch Program at DTU Wind Energy located in DTU Risø Campus in Denmark. 15
HAWC2 comprises the interaction between the elastic and inertial forces (related to the stiffness and motions of the body) with the aerodynamic forces derived from the exposure to the wind flow. Additionally, the inclusion of the control theory in this interactive system derives the denomination of aero-servo-elastic code which can be applied to define HAWC2. Hence, by means of this code, the time domain analysis where the elasticity of the structure (deflections), the structural dynamics (motion) and the aerodynamic loads can be simulated in a fully coupled environment. Furthermore, additional hazards encountered in offshore environments, turbulent turbine wakes encountered in wind farms, turbine control operations or seldom extreme occurrences like earthquakes lead to an even more complex coupled system. HAWC2 is designed to be capable of assessing the turbine response in time domain for such a complex scenario in a refined manner.
Using HAWC2, the structure is modelled based on flexible multi-body dynamics (MBD) approach. Along these lines, several independent bodies are coupled together and the global response can be analysed with increasing accuracy. The combination of the multiple bodies is performed by applying joints, which are algebraically constraints and may either represent a fixed relative position for the linked bodies or restrain certain degrees of freedom (DOF-s). Through this procedure, the overall motion of the entire flexible entity is divided into two separate parts. One regarding the rigid body motions, which is represented by the motions referred to the floating frame and the other corresponding to elastic motions which are superimposed. By means of the floating reference method the large motions of some elements such as the rotation of blades will be easily combined to the elastic deformations originated from loads [28]. In the multi-body system applied in order to model the NREL 5MW OWT, the following six main bodies have been defined. The monopile which represents the foundation. The tower, which is the column installed on top of the monopile.
16
Tower-top body, which intends to model the properties of the joint between the supporting media and the horizontal drivetrain. Connected to the tower-top is the shaft which aims to model the nacelle and the connection to the rotor. The hub stands as the last connection between the shaft and blades. Finally, three blades are connected to the hub what concludes the multi-body system. Additionally, in order to take flexibility into account, each main body of the system is split into several segments, which are modelled through the Timoshenko beam elements. Timoshenko beam elements consist of six degrees of freedom and thus, they enable capturing three translations (along the three main axes, x-y-z) and three rotations (around the aforementioned elements axes). Numerous Timoshenko beam elements were defined based on the geometrical and material properties of the structural members of the wind turbine. However, elements with relatively constant geometrical characteristics, like the monopile foundation and the tower, can be defined with a limited number of Timoshenko beam elements while, on the other hand, several of them are needed to model the turbine’s blades because of the constantly changing geometry along blade’s length. Regarding the supporting conditions that were adopted herein, the monopile foundation is assumed as fixed-base at its last node. Hence, the soil flexibility and the effect of the dynamic soil-structure-interaction is disregarded within the framework of the current MSc thesis. It is also noted that the HAWC2-based structural model, used herein, has been fully developed by the DTU Wind Energy Department. The accuracy of the particular model is ensured by tis extensive validation undertaken by the experienced researchers of the DTU Wind Energy Department. Furthermore, this model is the reference model for several relevant studies (e.g., [11]).
Since the structural model of the 5 MW offshore wind turbine has been established, HAWC2 enables the calculation of the dynamic characteristics of the structure through the eingenvalue (or modal) analysis. Especially, the natural frequencies of the first side to side and fore-aft. modes are listed in Table 2-3. As mentioned above, a fixed-base model of the 5MW offshore wind turbine was considered herein.
17
Table 2-3: Modal natural frequencies and modal damping ratio of the MBD system 𝑓𝑛 [𝐻𝑧]
𝑇𝑛 [𝑠]
𝜁 [%]
1st Tower fore-aft. mode (y dir.)
2.7568E-01
3.627
0.209
1st Towerside to side mode (x dir.)
2.7976E-01
3.575
0.212
1st Drivetrain torsion
6.0825E-01
1.644
0.414
1st Blade collective flap
6.4017E-01
1.562
0.463
1st Blade asymmetric flapwise pitch
6.7489E-01
1.482
0.498
1st Blade asymmetric flapwise yaw
9.9690E-01
1.003
0.466
1st Blade asymmetric edgewise pitch
1.0118E+00
0.988
0.473
1st Blade asymmetric edgewise yaw
1.5771E+00
0.634
0.996
In the study conducted by Jonkman et al. [11] several aeroelastic codes were compared for the NREL 5MW wind turbine and the modal damping was found to be controversial for codes based on MBD theory like HAWC2 when compared to modal theory based codes (i.e., FAST or Bladed). In that study, the modal damping could not be matched between MBD and modal codes without leading to numerical problems consequence of too low Rayleigh coefficients. For the case of HAWC2 (based on MBD) the Rayleigh damping is implemented. Consequently, Rayleigh coefficients are determined in a three-dimensional basis for each of the main bodies defined in the multi-body system. Then, in order to derive the modal damping ratios of the whole infrastructure, the contribution from each body is summed. Nevertheless, as can be found in the manual for the use of this HAWC2 [26], when a multi-body system is considered the mass proportional damping is found to be limited (see eq. (2.1)). 𝑪 = 𝛼𝑴 + 𝛽𝑲
(2.1)
where 𝑪 is the Rayleigh damping, 𝑴 refers to the mass matrix, 𝛼 is the mass proportional Rayleigh coefficient, 𝑲 is the stiffness matrix and hence 𝛽 is the stiffness proportional Rayleigh coefficient.
18
In order to provide some insight for the applied Rayleigh damping in the NREL 5MW model for time domain analysis in HAWC2, mass and stiffness proportional coefficients are presented in Table 2-4 for each body. Note that these values are provided incorporated in the model developed and provided by the DTU Wind Energy Department. Table 2-4: Rayleigh damping coefficients for the main bodies of the multi-body-system 𝛼𝑥
𝛼𝑦
𝛼𝑧
𝛽𝑥
𝛽𝑦
𝛽𝑧
Monopile
4.50E-2
4.50E-2
8.00E-1
1.20E-3
1.20E-3
4.50E-4
Tower
6.45E-4
6.45E-4
1.25E-3
1.40E-3
1.40E-3
1.25E-3
Tower-top
2.50E-4
1.40E-4
2.00E-3
3.00E-5
3.00E-5
2.00E-4
Shaft
0.00
0.00
7.0725E-3
4.65E-4
4.65E-4
7.0725E-3
Hub
2.00E-5
2.00E-5
2.00E-4
3.00E-6
3.00E-6
2.00E-5
0.00
0.00
0.00
1.41E-3
2.39E-3
4.5E-5
Body
Blade
In this section of the document relevant information is provided in order to explain the environmental and earthquake load definition and implementation techniques in HAWC2
Apart from the structural model, the accurate simulation of the loads subjected to the wind turbine is essential in order to generate precise and meaningful output. Through this section, the environmental loads, i.e., the wind and wave loads, which are relevant for the offshore environment, are thoroughly described. Moreover, modelling techniques are also presented, which are related to the loads definition within the HAWC2 framework for the time domain analysis.
19
For any time-domain simulation of a wind turbine using the aero-elastic code HAWC2, wind conditions can be defined on the basis of three main parameters: the wind profile, the mean wind speed and the turbulence intensity. More specifically, the wind profile determines how the mean wind speed varies in height and the turbulence intensity is a magnitude characterizing the turbulence of the wind field. The latest is defined as the standard deviation of the wind speed with respect to its mean value. Based on the aforementioned characteristics, the wind action can be fully determined for any height within the atmospheric boundary layer, i.e., from zero wind speed found on the crust of the Earth up to the geostrophic wind velocity. In order to define these wind field characteristics, IEC 61400-3 [29] is considered as the reference code, where design requirements for offshore wind turbines are also prescribed. With regard to the wind profile, it is assumed to follow the power law which determines the mean wind speed in function of the height as defines the eq. (2.2), 𝑧 𝛼 𝑉(𝑧) = 𝑉ℎ𝑢𝑏 ∙ ( ) 𝑧ℎ𝑢𝑏
(2.2)
where 𝑉ℎ𝑢𝑏 is the mean wind speed at hub height (averaged 10 minutes wind speed), 𝑧ℎ𝑢𝑏 is the height from the mean water level (MWL) to the hub, 𝑉(𝑧) is the mean wind speed at any 𝑧 height and 𝛼 is the power law exponent, which is defined for offshore normal conditions equal to 0.14. Regarding to the turbulence intensity, it is defined based on eq. (2.3) following the prescription of IEC 61400-3 [29].
𝑇𝑖𝑛𝑡 =
𝜎1 𝑉ℎ𝑢𝑏
(2.3)
where 𝜎1 is the standard deviation of the longitudinal (incident) wind speed component, which can be determined applying the following expression shown in eq. (2.4).
𝜎1 =
𝑉ℎ𝑢𝑏 + 1.28 ∙ 1.44 ∙ 𝐼15 𝑧 𝑙𝑛 ( ℎ𝑢𝑏 ) 𝑧0
(2.4)
20
Where 𝐼15 is the average value of the turbulence intensity at hub height which is determined for 𝑉ℎ𝑢𝑏 = 15 𝑚/𝑠. In accordance with IEC 61400-1 [20], the I15 parameter is equal to the reference turbulence intensity (𝐼𝑟𝑒𝑓 ), and assuming that for open sea conditions the reference turbulence intensity can be the minimum for the onshore case, the Iref is eventually equal to 0.12. Furthermore, the parameter 𝑧0 is the roughness length that represents the extrapolated height, where zero wind speed is found for a logarithmic variation along the height. This parameter is obtained by an iterative process leading to the minimum value of Charnock expression, as shown in eq. (2.5), 2
𝐴𝑐 𝜅 · 𝑉ℎ𝑢𝑏 𝑧0 = [ ] 𝑔 𝑙𝑛 (𝑧ℎ𝑢𝑏 ) 𝑧0
(2.5)
where 𝑔=9.81 m/s2 is the acceleration originated by gravity, 𝜅=0.40 represents the von Karman coefficient and 𝐴𝑐 is the Charnock´s constant which for this case where open sea conditions are assumed leads to 𝐴𝑐 = 0.011. Once the necessary wind conditions are defined, the wind field can be generated within the framework of HAWC2 aero-elastic code based on the Mann turbulence model, which fulfils the IEC 61400-3 [29] requirements. Particularly, using the HAWC2 aero-elastic code, the turbulent atmospheric wind can be simulated by the Mann model. Along these lines, the atmospheric wind state is modelled with three-dimensional spectral tensors, whom spectral energy distribution is based on the Von Karman spectrum. The latter is is defined in eq. (2.6),
𝐸(𝑘) =
2 5 𝛼𝜀 3 𝐿3
(𝐿𝑘)4
(2.6)
17
(1 + (𝐿𝑘)2 ) 6
2𝜋𝑓 ̅ 𝑈
where 𝐸(𝑘) is the energy concentrated at a specific frequency, 𝑘 =
is the wave number,
̅ and the frequency, f of the energy spectrum, 𝐿 is the length which relates the mean wind speed 𝑈 scale, which is recommended to be calculated as 𝐿 = 0.7Λ1 (based on codes prescriptions [20], 2
Λ1 = 42 𝑚 for heights above 60 m. Finally, the 𝛼𝜀 3 input parameter enables the analytical 2
expression of the longitudinal variance of the wind speed (𝛼𝜀 3 = 1).
21
Finally, in order to ensure a proper analytical relationship for the longitudinal variance of the wind, the scaling factor shown in eq. (2.7) is multiplied to every vector in the Mann turbulence box in a three-dimensional field (𝑢, 𝑣, 𝑤), see Figure 2-1. The definition of the scaling factor includes the turbulence intensity parameter. The latter along with the mean wind speed constitute the two variables, which are necessary to define the turbulence box according to the Mann model.
𝑆𝐹 = √
2 𝜎𝑡𝑎𝑟𝑔𝑒𝑡 𝜎2
(2.7)
Where 𝜎𝑡𝑎𝑟𝑔𝑒𝑡 is derived from the requested turbulence intensity and 𝜎 is the actual standard deviation of the wind turbulence calculated at the middle of the turbulent box.
The definition of the wind field and its critical parameters (see above), enable the calculation of the aerodynamic loads that are imposed to any structural member being subjected to the wind flow. For a wind turbine, the blades are mostly exposed to the aerodynamic loads, while the tower is marginally influenced by them. Next, the theoretical background for the calculation of the aerodynamics loads is shortly provided. Based on the actuator disc conception, the air flow approaching the turbine experiences a reduction in the axial speed, which causes the downstream flow to slow down. However, it does not get compressed and it leads to the conclusion that based on the fact that the mass flow rate must remain constant, the volume of the mass crossing the rotor section increases. This change of the velocity is quantified by the axial flow induction factor shown in eq. (2.8), where the wind speed at the disc (𝑈𝐷 ) and upstream (𝑈∞ ) are related.
𝑈𝐷 = 𝑈∞ (1 − 𝑎)
(2.8)
From eq. (2.8), the simple momentum theory is derived, where the change in momentum is related to the variation of velocity the pressure difference across the disc actuator, widely known as thrust. On the other hand, as a consequence of the axial momentum decrement, the drop in the energy content on the air flow is transferred to the rotational kinematic energy of the blades which make the shaft rotate. In this case, the shaft experiences two torques in 22
opposite directions. One originated from the air flow and the other by the generator. As for the first one, it requires an identical and contrary torque to be foisted on the air. Hence the air particles undergo a change in angular momentum. Hence, apart from the axial velocity change, there is a tangential velocity change, which is quantified by the tangential flow induction factor (𝑎′), as expressed in eq. (2.9)
𝑎(1 − 𝑎) = 𝜆2𝑟 · 𝑎′
(2.9)
where 𝜆𝑟 is the local speed ratio obtained as the ratio between the multiplication of tangential velocity of rotating circular ring (𝑟) and the rotational velocity of the shaft (Ω) over the upstream air flow velocity. Based on the simple momentum theory and the angular momentum theory explained before, the resultant relative velocity acting on the blade can be derived in the form shown in eq. (2.10), which is required in order to determine the angle of attach of the wind.
2 (1 − 𝑎)2 + 𝑟 2 Ω2 (1 + 𝑎 ′ )2 𝑊 = √𝑈∞
(2.10)
Once the angle of attach as well as the magnitude of the acting wind are determined, the lift and drag forces on a span-wise basis (𝛿𝑟) can be defined for each blade (see eq.s (2.11) & (2.12)) following the blade element theory’s assumption. This states that the aerodynamic forces (lift and drag), of two identical but isolated elements with the same angle of attach are identical in a two-dimensional flow. 1 𝜌𝑊 2 𝑐𝐶𝑙 𝛿𝑟 2
(2.11)
1 𝛿𝐷 = 𝜌𝑊 2 𝑐𝐶𝑑 𝛿𝑟 2
(2.12)
𝛿𝐿 =
where 𝛿𝐿 and 𝛿𝐷 are the lift and drag forces on span-wise, 𝜌 is the air density, 𝑐 is the chord length of the blade, while 𝐶𝑙 and 𝐶𝑑 are the lift and drag coefficients.
23
Hence, the axial thrust and the torque generated at each annular ring can be defined in a spanwise length which leads to the overall lift and drag loads [30]. Nevertheless, in order to solve the parameters leading to the determination of aerodynamic loads, pertinent assumptions proposed by the Blade-Element/Momentum (BEM) theory are necessary. More specifically, following the BEM theory, it is assumed that the flow affecting to each annuli is independent of the contiguous annulus flow. Therefore, this radially independent interaction enables the discretization of the blades into different annulus and the combination of the loads derived from each of the annuli considered. In reality, the axial flow induction factor does not remain constant in the span-wise direction (radially), which would make the BEM theory main assumption incorrect. However, the literature has already demonstrated that this assumption is acceptable and provides accurate enough results. It is notable that HAWC2 utilizes the BEM theory in order to derive the aerodynamic loads on the blades.
For wind turbines installed offshore, the hydrodynamic loads are of relevant importance. These loads are derived from current and wave actions and in HAWC2 are determined through the Morison’s equation. Additionally, in order to simulate an entire sea-state for the time window of the simulation, JONSWAP spectrum is utilized, which generates a number of linear waves. The superposition of the waves leads eventually to the calculation of the actual sea-state made up of irregular waves. On the other hand, for simplification reasons and due to its marginal influence to the dynamic response of the offshore wind turbine, the current is simulated as constant in depth instead of using a more realistic shear flow distribution.
With the purpose to define sea conditions, waves and current parameters have to be determined and applied within the HAWC2 framework. Regarding the currents definition, the current velocity has to be determined while the waves are defined on the basis of the significant wave height (𝐻𝑠 ) and the peak period (𝑇𝑝 ). Even though the direct relation between the significant wave height and the peak period is typically non-existent, according to the Guide to Wave Analysis and Forecasting [31], the peak wave period (or frequency) can be related, via eq. (2.13), to a specific significant wave height and the mean wind speed, u, calculated 10 m above the mean sea level
24
𝑓𝑝 = 0.148 ∙ 𝐻𝑚0−0.6 ∙ 𝑢0.2
(2.13)
where 𝑓𝑝 is the frequency of the peak period, 𝑇𝑝 = 1/𝑓𝑝 and 𝐻𝑚0 is the theoretical significant wave height related to zero moment of the variance spectrum distribution (𝑚0 ). The zero moment represents the area under the curve defining the spectrum which determines periods and amplitudes of the linear waves included on the irregular wave. In addition, commonly the spectrum width (𝜀) is considered to be 0.4-0.5, which leads to an accurate expression for 𝐻𝑠 = 3.7√𝑚0 [32]. Based on a deterministic 𝐻𝑠 value, the 𝑚0 can be calculated and inserted into 𝐻𝑚0 = 4√𝑚0 which enables the derivation of the accurate value for 𝑓𝑝 . The definition of the both the significant wave height (𝐻𝑠 ) and the peak period (𝑇𝑝 ) enable the simulation of the sea-state based on the JONSWAP spectrum. In accordance with the IEC 61400-3 [29], the wave modelling technique is based on the Joint North Sea Wave Project (JONSWAP) spectrum. This spectrum is usually applied to model waves in the growing phase. In the following equation the spectrum is defined.
−5
𝑆𝐽𝑆 (𝑗) =
0.3125𝐻𝑠2 𝑇𝑝
𝑓 ( ) 𝑓𝑝
4
𝑓 exp (−1.25 ( ) ) (1 − 0.287 ln 𝜆)𝜆 𝑓𝑝
2 𝑓 −1 𝑓𝑝 ) exp −0.5( 𝜎
(
(2.14) )
where σ is defined as follows:
𝜎={
0.07 0.09
𝑓𝑜𝑟 𝑓𝑜𝑟
𝑓 ≤ 𝑓𝑝 𝑓 > 𝑓𝑝
(2.15)
Furthermore, an important feature of the JONSWAP spectrum is its dependence on the peak enhancement factor (𝜆), which determines the sharpness of the spectrum (see eq. (2.16)). 5 𝑒𝑥𝑝 (5.75 − 1.15
𝜆=
{
1
𝑓𝑜𝑟 𝑇𝑝 √𝐻𝑠
)
𝑇𝑝
< 3.6 √𝐻𝑠 𝑇𝑝 𝑓𝑜𝑟 3.6 ≤ ≤5 √𝐻𝑠 𝑇𝑝 𝑓𝑜𝑟 >5 √𝐻𝑠
(2.16)
25
Based on the spectral distribution of the energy concentrated for each frequency related to a regular wave, the actual sea-sate can be calculated by superposition of uniform waves.
The definition of the wave field and its critical parameters (see above), enable the calculation of the hydrodynamic loads that are imposed to any structural member that is located underwater. Next, the theoretical background for the calculation of the hydrodynamic loads is shortly described. In order to quantify the hydrodynamic loads on the vertical cylindrical body (e.g., the tower of the wind turbine), the Morison’s equation can be applied [33]. More specifically, for steady current conditions, the inline directional force acting on the cylinder can be calculated from the following equation
𝐹=
1 𝜌𝐶 𝐷𝑈|𝑈| 2 𝐷
(2.17)
where 𝐹 denotes the force per unit length of the cylinder, 𝐶𝐷 is the drag coefficient, 𝜌 represents sea water density (1027 kg/m3) and 𝑈 defines the flow velocity. Nevertheless, when the steady state drag is combined with the oscillatory flow scenario derived from wave motions, the total inline force has to be coupled with two other contributions, which are the hydrodynamic mass and Froude-Krylov force. The hydrodynamic mass effect is originated by the equivalent principle of a moving body on a steady fluid. In such a case, what the body experiences is that the acceleration of the fluid in the vicinity of the body originated from the body’s motion causes the hydrodynamic mass forces, as defined in eq. (2.18). This is equivalent of the accelerations that the water would experience with a oscillatory flow and steady body in absolute terms.
𝐹ℎ.𝑚𝑎𝑠𝑠 = 𝜌𝐶𝑚 𝐴𝑈̇
(2.18)
where 𝐶𝑚 is the so called hydrodynamic mass coefficient (circular cross section 𝐶𝑚 = 1), 𝐴 is the cross sectional area of the cylinder and 𝑈̇ denotes the acceleration as first partial derivate of the velocity, 𝑑𝑈/𝑑𝑡. 26
On the other hand, the Froude-Krylov force accounts for the effect derived from the outer-flow region. There, the accelerated motion of the fluid causes a pressure gradient, which can be expressed in the form of eq. (2.19).
𝐹𝑝 = 𝜌𝐴𝑈̇
(2.19)
Therefore, the Morison force can be derived combining eq.s (2.17), (2.18) & (2.19). However, the equation is further refined in order to account for the motions of the structure due to loading actions (see eq. (2.20)). That is done by calculating the relative velocity of the motion in fluidbody basis (𝑈 − 𝑈𝑏 ).
𝐹=
1 𝜌𝐶 𝐷(𝑈 − 𝑈𝑏 )|𝑈 − 𝑈𝑏 | + 𝜌𝐶𝑚 𝐴(𝑈̇ − 𝑈̇𝑏 ) + 𝜌𝐴𝑈̇ 2 𝐷
(2.20)
Regarding the seismic motions, their inclusion into the aero-elastic code has been performed through advanced dynamic linked libraries (dll files). These dll files are compiled scripts coded in FORTRAN 90 and allow the implementation of the external, earthquake-imposed actions at the base of the structure. More specifically, within the HAWC2 framework, the earthquake motions are defined in the form of ground accelerations that can be potentially applied along the two main horizontal directions and the vertical one. Hence, due to the acceleration input in the last node of the wind turbine (i.e., the last node of the monopile), the structure will experience the motion of the earthquake and its displacements that will be applied for cross sectional forces determination. The equation of motion, considering only the earthquake excitation and disregarding the complexity of the turbine as well as the other environmental loads, is expressed in eq. (2.21), where 𝑀𝑒 𝑢̈ 𝑔 represents the force in form of the accelerated equivalent mass.
𝑀𝑥̈ + 𝐶𝑥̇ + 𝐾𝑥 = −𝑀𝑒 𝑢̈ 𝑔 (𝑡)
(2.21)
27
This chapter aims to provide the necessary validation related to the use of HAWC2 aero-elastic code for the seismic analysis, through real strong ground motions, of a structural model. Furthermore, sensitivity analysis was conducted in order to define with increased accuracy critical parameters for the time-domain simulation of the offshore wind turbine when subjected to a multi-hazard environment consisting of wind and waves loads and the earthquake-imposed ground accelerations.
The HAWC2 aero-elastic code is primarily dedicated to run time-domain analysis for a wind turbine model when subjected to wind and wave loads. However, the development and utilization of the aforementioned dll files, related to the application of the external forces at the base of HAWC2 model, and their interaction with the core of the HAWC2 code enable the incorporation of seismic forces into a time-domain analysis. In other words, HAWC2 can be additionally used for seismic analysis of a wind turbine, when the earthquake excitations are provided as time series of strong ground motions (i.e., accelerations). Still, there is no existing documented effort to run seismic simulations within the framework of HAWC2. Hence, it is necessary to materialize a validation scheme in order to be sure that the current code can run reliably and with accuracy seismic simulations. Based on the above, two real earthquake records (i.e., time-series of strong ground motions recorded during past seismic events) have been selected and seismic simulations for two different structural models were performed using both HAWC2 and SAP2000, a widely used finite element code especially in earthquake engineering. The comparison from the results, which were obtained from these analyses (i.e., two seismic records, two structural models, two
28
finite element codes), will enable to understand if the HAWC2-based simulation can be reliably used for the seismic analysis of structural systems. Along these lines, two different structural models were created using both HAWC2 and SAP2000: (a) a long mast (high-rise tower) and (b) a shorter tower. The geometrical characteristics and the material properties of these two structural systems, which were modelled identically in HAWC2 and SAP2000, are listed in Table 3-1. In addition, the perfect matching of their dynamic characteristics (derived after eigenvalue analysis both in HAWC2 and SAP2000 and presented in Table 3-2), verifies that the models created using these two finite element codes are almost identical in terms of geometrical, material and mass properties. Table 3-1: Geometrical and material properties of the validation models
Height [m]
Mast
Short column
66
8.0
Elasticity modulus [N/m2]
2.10E11
Shear modulus [N/m2]
8.0768231E10
Density [kg/m3]
7850
Boundary cond.
Fixed base
Mesh, FEM Diameter, solid C.S. [m]
44 FE, 1.50 m each
32 FE, 0.25 m each
0.60
0.25
Furthermore, the modal damping corresponding to the first 40 modes was harmonized for the both HAWC2 and SAP2000 models. This harmonization of the structural damping for both models was performed based on the modal logarithmic decrement obtained through HAWC2. Considering this logarithmic decrement as the basis, the modal damping ratio was derived following the expression shown in eq. (3.1) and provided by [34].
𝛿 ≃ 2𝜋𝜁
(3.1)
29
Hence, by implementing these modal damping ratios in the short as well as the high-rise models built up in Sap2000 the global harmonization of the dynamic characteristics of the structures was achieved. Evidence of the harmonization are the highly similar first eight natural frequencies presented in Table 3-2. Table 3-2: Natural Eigen frequencies of validation models, [Hz] Mast (high-rise tower)
Short tower
HAWC2
Sap2000
HAWC2
Sap2000
𝑓1
0.0965
0.0996
2.824
2.8229
𝑓2
0.0965
0.0996
2.824
2.8229
𝑓3
0.6242
0.6238
17.604
17.586
𝑓4
0.6242
0.6238
17.604
17.586
𝑓5
1.7465
1.7445
48.892
48.811
𝑓6
1.7465
1.7445
48.892
48.811
𝑓7
3.4189
3.4134
94.695
94.489
𝑓8
3.4189
3.4134
94.695
94.489
Based on the above, since the models (from HAWC2 and SAP2000) are identical, they should also lead to common response results if they are subjected to the same loading conditions (i.e., common earthquake excitation). Table 3-3 provides the details of the two seismic motions that were selected to excite the structural models described above. Additionally, the time history of the selected strong ground motions, recorded in two horizontal directions (x and y directions) as well as their corresponding response spectra (𝜁 = 5%) are depicted in Figure 3-1.
30
Table 3-3: Characteristics of Kobe and Manjil earthquakes applied for validation purposes Earthquake
Magnitude
Epic. (km)
Soil (NEHRP)
PGA (x,y), g
Record st.
Kobe, Japan
6.90
25.40
B
0.290/0.310
Kobe
(1995) Manjil, Iran
Univeristy 7.37
77.84
D
0.132/0.209
Abhar
(1990)
The seismic, time domain analysis of both the mast and the short column, modelled by HAWC2 and SAP2000, enabled calculating the seismic response in terms of 𝑥, 𝑥̇ , 𝑥̈ , 𝐹𝑥,𝑦 and 𝑀𝑥,𝑦 . Both the deformations (displacement, velocity and acceleration) and the forces were calculated for several different locations along the height of the structures.
Figure 3-1: Time history of strong ground motions and response spectra of Kobe and Manjil earthquakes As can be seen in Figure 3-2, the bending moments, calculated at the base of both structural systems when they subjected to both the seismic motions, show an almost perfect matching. The same trend was also observed for the entire set of response results that were investigated. Hence, it is profound that the HAWC2 code can reliably predict the seismic response of different structural systems, since it provides common response results with a widely used finite element code (SAP2000) for earthquake engineering purposes.
31
Figure 3-2: Moment demand at base frame, HAWC2 vs. Sap2000, Manjil & Kobe earthquakes. (a) refers to mast and (b) to short column response.
The influence of the initial conditions of the wind turbine (rotational speed of the rotor’s shaft, Ω [rad/s], and the initial pitch angle of the blades 𝜃 [rad]) as well as the time step applied for the time-domain simulations have been thoroughly studied. From this study, the initial conditions of the wind turbine have been found to be of high influence. The further are the initial angular velocity of the shaft and the pitch angle from the correspondent to the simulated wind speed, the larger is the motion pulse generated by the wind on the structure in the first time-step which leads to a detrimental transient with regard to the simulation’s accuracy. This is because wind acts like an impulse when the simulation starts, at time zero the turbine has no loads while for the first time-step the air flow is hitting the structure. So, the transient part is product of this initial instability, which after some time finds equilibrium.
32
Therefore, if the transient part is too large, a small modification of the time-step leads to a drift at the peak of the transient part, which derives into different equilibrium conditions for different time-steps. These conditions mainly concern the mean displacement generated by the wind at the tower top as well as the moment at the monopile base. Hence, in order to avoid such unsustainable conditions, the analysis based on the wind ramp has been implemented. With this way, the specific initial conditions depending on the mean wind speed applied for the simulation are known and hence they can be implemented for each simulation, preventing the drift of the response results once after the transient part.
Figure 3-3: Wind ramp (𝑈𝑚 ), pitch angle (𝜃) & shaft rotational velocity (𝛺) In Figure 3-3, the incremental procedure of applying the wind speed without no turbulence is shown in conjunction with the turbine performance. It is clearly visible, that when the wind speed exceeds the rated wind speed the controller starts pitching in order to keep constant the angular velocity of the shaft and thus the power production. From that analysis, the response of the wind turbine, once the equilibrium is reached for each wind speed, is adopted as the initial condition for further simulations at that specific mean wind speed. Applying the initial conditions obtained from Figure 3-3, the problem with the dramatic transient is solved. However, the choice of an acceptable time-sept for the analysis have been answered by a further sensitivity analysis. This analysis intends to compare the relative accuracy of the results with regard to the applied time-step and the corresponding computational time cost. This is done based on the assumption that for the lowest time step applied (𝑑𝑡 = 0.0005 𝑠) the relative error is negligible. The simulations for this purpose have been run 33
with the most numerically instable conditions that might be intended to analyse through this thesis. That is why, the significant wave height is 𝐻𝑠 = 6.5𝑚, the mean wind speed 𝑈𝑚 = 25𝑚/𝑠 and the time step of the seismic motion were applied is 𝑑𝑡 = 0.005𝑠. Based on those environmental loads, the moment at the base of the monopile has been considered and the mean value after the transient part was calculated for each of the applied time-steps. The 𝑑𝑡 has been doubled seven times covering the range of 𝑑𝑡 = 0.005 − 0.064 𝑠.
Figure 3-4: Time-step influence on HAWC2, relative error vs time cost It can be seen from Figure 3-4 that the Newmark time integration procedure becomes instable after a value of approximately 0.03s for the time step. Apart from that, considering that the lowest time step for the earthquakes that are intended to be analysed is 0.005s, and based on the fact that this time-step reflects a good relationship between accuracy and time cost, it has been selected to be applied for further simulations.
34
Within this chapter, the response of a wind turbine when subjected to an earthquake is the theme of interest. The magnitude of the demand induced by the earthquake is studied and the critical locations on the structure, in terms of maximum response, are intended to be determined. Nonetheless, the turbine is not simply an object but a subject. A complex structure designed to interact with the wind in order to optimize the power production. For this reason, the response is not only dependant on the environmental load, but also the interaction of the turbine state and its controller.
In order to estimate the response of a wind turbine when it is excited by a seismic event, the operational state of the turbine is essential. In this document three different configurations are considered for the preliminary study.
This operational state is defined by the no opportunity for the wind turbine to shutdown for any reason. In this way, the turbine keeps operating in spite of the earthquakes.
On the other hand, the main feature of this second operational state main feature is that the acceleration sensor at the hub states a threshold of 1m/s2, which in case of being exceeded a shutdown mechanism is activated. This shutdown consists of active blade pitching with an angular speed of up to 4 deg/s in the final phase. In This way, the torque originated by the wind action on the blades almost disappears smoothly but the electrical torque acts against the remaining rotation, which leads to a soft decay of the angular velocity of the rotor. When the 35
shutdown protocol is activated, the exposed area of blades to the wind is decreased leading to a progressive transition from operational mean displacement downstream of the hub to the zero mean of the parked condition.
In order to model the parked condition, the fixed parked condition has been specified, where not rotation of the shaft is allowed. The initial pitch angle reaches the maximum, which leads to the minimum exposed area of the blades perpendicular to wind stream.
With the purpose of analysing the performance of a wind turbine when it is excited by the earthquake-imposed strong ground motions, a specific loading scenario is determined for each of the operational states described above. Especially, three different reference case scenarios are considered in order to get more insight into the influence of the turbine state on the overall response. These three cases are defined in Table 4-1 as well as five earthquakes in Table 4-2. Table 4-1: Environmental load specification for reference case scenarios Operating
Operating & shutdown
Parked
𝑈𝑚 = 13𝑚/𝑠
𝑈𝑚 = 13𝑚/𝑠
𝑈𝑚 = 35𝑚/𝑠
𝐻𝑠 = 6.5𝑚
𝐻𝑠 = 6.5𝑚
𝐻𝑠 = 6.5𝑚
Table 4-2: Five earthquakes applied for seismic response study of a wind turbine Earthquake
Magnitude
Epic.(km)
Soil NEHRP
PGA (x,y), g
Record st.
Nahanni, Canada
6.76
6.80
C
0.978/1.096
Site 1
Northridge-01
6.69
20.27
D
0.583/0.590
Newhall – Fire sta
Smart, Taiwan
7.30
76.21
D
0.122/0.153
Smart 1 C00
Superstitionhills-02
6.54
15.99
D
0.455/0.377
Parachutetest site
Tabas, Iran
7.35
55.24
B
0.836/0.852
Tabas
36
Figure 4-1: Time domain response analysis of a wind turbine subjected to an earthquake
37
Figure 4-1 depicts the time history of the response (deflection tower-top) of the OWT subjected to five different earthquakes. From this time domain analysis conducted for five different earthquakes, it is noticable the gobernance of the earthquake on the reponse of the wind turbine for any of the three wind turbine case scenarios studied. In addition, the effect of the wind is appreciable to be much larger on the operating and the operating & shutdown cases than for the parked wind turbine. This is given by the smaller pressure that the blades of the wind turbine are subjected to when the parked condition is adopted and hence the exposure area of the blades becomes the minimum. Furthermore, as can be identified mainly from the third column of the diagrams in Figure 4-1, the maximum response is generally found for the operational case. Nonetheless, for those earthquakes that excite the wind turbine more severely (e.g., containing long predominant period) the parked case seem to reach even exceed the maximum response derived from an operational wind turbine. This is the case for the Tabas earthqueake shown in Figure 4-1. In order to study further the maximum response of a wind turbine excited by an earthquake, a singular case which subjectst the wind turbine to substantial degree of demand has been selected deliberately. Manjil earthquake, utilized above in the validation process in chapter 3.1 has chosen to represent the effects of the earthquake along the height of a wind turbine.
Figure 4-2: Tower-top deflection in X and Y & moment at the base of MP around X and Y (operating, shut-down and parked states)
38
For the case of Manjil earthquake depicted in Figure 4-2, the parked wind turbine shows, as expected and in agreement with Figure 4-1, much lower damping than the two operational cases, where the aerodynamic damping has an important role in the dynamic response of the wind turbine as can be observed by the smooth decrement of the response along the Y direction. Furthermore, the response in X direction of the parked wind turbine is limited compared to the other two cases. That represents how the motion of the structure excited mainly by the earthquake is modified by the wind in its parallel direction forcing the oscillation direction to turn onto the other axis. This reflects how different is the motion in amplitude and direction of the turbine depending on the operational state conditions. On the other hand, it can be perceived that even though the shutdown mechanism attenuates the response, the maximum response of the turbine is very similar to the constant operational condition. Nevertheless, the shut-down case shows smaller maximum response for the Manjil earthquake. Therefore, from Figure 4-2 it can be concluded that the most detrimental condition happens to be the parked one while the difference in terms of maximum response of the operational and operational & shut-down is almost negligible. Furthermore, Figure 4-3 shows the response along the height of the wind turbine structure in maximum terms when it is subjected to the Manjil earthquake. Hence, the critical spots with regard to different demands can be identified. With regard to the maximum shear (𝑉𝑥&𝑦 ), moment (𝑀𝑥&𝑦 ) and displacement (𝐷𝑖𝑠𝑝.𝑥&𝑦 ) induced by the seismic hazard, the response can be predictable based on the previous Figure 4-2. However, the acceleration in both directions (𝐴𝑐𝑐.𝑥&𝑦 ) show that second vibration mode is activated by the fully coupled load scenario together with the first vibration mode as represented by the displacement plots. For Y direction, the influence of the second vibration mode is significantly larger than for the vibrations in x direction, where even higher vibration modes can be identified. This results are aligned with the conclusions made by Alati et al. [13] who found significant accelerations at approximately two thirds (2/3) of the total height. It can be seen that in [13] the stiffness of the foundation (tripod and jacket) was larger than in this analysis (monopile) and that makes the shape of the maximum acceleration response in height smoother for this more flexible structure.
39
Figure 4-3: Maximum and minimum response along the height of the wind turbine to a single seismic motion, operating (Op.), shut-down (Sh.d) & parked (Pa.) Despite the fact that three different wind turbine states have been analysed along this section, it has to be considered that the shut-down scenario is not a straightforward solution. The characteristics of the shut-down triggered by the earthquake depends significantly on the design of the controller and the strategy followed to proceed with this protocol. In addition, for the shut-down scenario applied above the maximum response was always found for the constantly operating case.
40
In this chapter of the report, the fragility that the horizontal axis offshore wind turbine studied herein experiences due to its exposure to a multi-hazard environment that consists of wind, waves and earthquake excitations. In order to quantify the fragility, the probability of exceeding a certain limit state threshold or damage condition has to be determined. With this aim, fragility curves are derived, where the capacity of the turbine to withstand the demand related to a specific loading scenario and varying intensity measure is estimated. Besides the probabilistic background which is necessary for the fragility calculations, this chapter includes also the definition of: (a) the loading parameters that were used to compose the multi-hazard environment adopted herein and (b) the capacity limit states based on different failure modes that were investigated in this study. The convolution of all the aforementioned parts enables the calculation of the fragility curves for the specific wind turbine.
In this study the vulnerability related to the wind turbine foundation structure is analysed. In this case, the monopile is mainly considered to be seriously threatened by the multi-hazard environment adopted in this study. Due to the lack of redundancy on the structural system, one single failure mode within the monopile is enough to result in collapse or non-operative condition consequences. Moreover, as a result of the fact that the monopile is the only structural member connecting the tower to the supporting media, retrofitting is rarely a feasible option when the MP experiences certain degree of damages that would require replacement. That makes even more dramatic the need to assess the probability of failure of this element in detail. Additionally, the fragility of the nacelle, where highly sensitive devices are located, is also assessed on the basis of acceleration-based capacity limit states. No fragility analysis was performed for the blades within the framework of this study. In structural risk assessment, different structures such as bridges [35] & [36] or buildings [37] have been addressed based on fragility curves. The theory behind fragility curves is derived as
41
a brunch from structural reliability. The structural vulnerability is assessed with regard to some parameter defining the magnitude and characteristic of the external hazard. Within the framework of this report, fragility functions have been defined to determine the probability of the demand (D) reaching or exceeding a certain capacity limit or damage state (C) for a given intensity measure (IM) of the seismic hazard as is represented in eq. (5.1).
𝐹𝑟𝑎𝑔𝑖𝑙𝑖𝑡𝑦 = 𝑃[𝐷 ≥ 𝐶|𝐼𝑀] = 𝑃[𝐶 − 𝐷 ≤ 0.0|𝐼𝑀]
(5.1)
Following this definition, the probability distribution of the demand, (D), has to be related to the selected intensity measure, IM. For the case of earthquakes hazard, this is named probabilistic seismic demand model, referred as PSDM from this point. Cornell et al. [38] showed that the approximate median demand (𝑆𝑑 ) can be derived from a power model expressed in eq. (5.2).
𝑆𝑑 = 𝑎 · 𝐼𝑀𝑏
(5.2)
where 𝑆𝑑 is the median estimate of the selected demand parameter, 𝑎 and 𝑏 are coefficients that are obtained from regression analysis for the response results (i.e., demand parameters like displacements, accelerations, forces, stresses) that have been calculated after running several time-domain analysis with seismic excitations of varying intensity, quantified through appropriate 𝐼𝑀𝑠. In the following chapters, the regression analysis that was applied herein and the related PSDM are presented in figures. When the power model is estimated and the regression coefficients defined (𝑎 &𝑏), the distribution of the demand about its median can be determined assuming a two parametric probabilistic lognormal distribution [38]. However, a third parameter is required and this is the logarithmic dispersion, 𝛽𝐷|𝐼𝑀 , of the demand with regard to its median. This dispersion is obtained from eq. (5.3). 𝑛
𝛽𝐷|𝐼𝑀
= √∑ 1=1
(𝑙𝑛(𝑆𝑑 𝑖 ) − 𝑙𝑛 (𝑆𝑑 𝑟𝑒𝑔. ))
(5.3)
𝑛−1
42
where 𝑆𝑑 𝑖 is the actual demand obtained from numerical simulations, 𝑆𝑑 𝑟𝑒𝑔. is the demand obtained from the regression analysis. From this point, the conditional seismic demand is modelled through a lognormal probability distribution as was applied in [39] & [40] and is shown in eq. (5.4). ln(𝑑) − 𝑙𝑛(𝑎𝐼𝑀𝑏 ) 𝑃[𝐷 ≥ 𝑑|𝐼𝑀] = 1 − Φ ( ) 𝛽𝐷|𝐼𝑀
(5.4)
where 𝐷 is the demand, which is related to a probability value of occurrence, 𝑑 is the actual demand related to an intensity measure (reference for threshold), 𝛽𝐷|𝐼𝑀 is the logarithmic standard deviation or dispersion about the median estimate and Φ(●) is the standard normal cumulative distribution function. Nevertheless, when structural reliability is the matter of interest, it is important to consider that the capacity of the structure is not one fixed value but has a probability distribution. In order to account for this probabilistic distribution of the capacity of the structure, these limit state capacities are determined related to the demand that the structure will be exposed during the earthquake hazard. Hence, the fragility can be eventually represented in the form shown in eq. (5.5). 𝑙𝑛(𝑆𝑑 /𝑆𝑐 )
𝑃[𝐷 > 𝐶|𝐼𝑀] = Φ
2 √𝛽𝐷|𝐼𝑀
[
+ 𝛽𝑐2
(5.5) ]
where 𝑆𝑐 is the median of the capacity calculated from its lognormal distribution and 𝛽𝑐 is the logarithmic standard deviation of the capacity. It is noted that the aforementioned probabilistic procedure to estimate the fragility of a structurral system, which is exposed to seismic actions, has been widely used in numerous studies from the earthquake engineering field.
The second component that is required for the fragility analysis is the definition of the capacity limit states along with the corresponding demand estimates. Within the framework of this MSc thesis, two main capacity limit states were considered: (a) the Ultimate Limit State (ULS) and (b) the Serviceability Limit State (SLS). The ULS concerns the overall structural integrity of 43
the wind turbine and more specifically, the monopile. This capacity limit state incorporates herein two main failure modes that are expected for the monopile based foundation. More specifically, the plastic failure mode and the local buckling failure mode. The latest has a significant importance due to the fact that as a consequence of the geometry of the monopile (cross section class 4 according to EN 1993-1-1 [41]), it has to be treated as a shell. In addition, within the scope of the plastic failure three different failure modes have been studied. More accurately, the bending moment failure mode, the shear failure mode and their combination through the widely known von-Mises criterion were considered within the framework of the ULS. Regarding the second limit state that was considered, SLS, it was evaluated through this study on the basis of acceleration-related thresholds for the top of the tower (nacelle), where the damage in the highly-sensitive devices located there can cause serious operational problems and malfunctions for the wind turbines. Next, the capacity limit states, which were adopted in this study, will be thoroughly described.
The distribution of the loads derived from seismic events decrease along the height of the monopile and hence, the wind turbine’s base (i.e., bottom of the monopile) experiences the maximum demand (at least in terms of sectional forces like the bending moments and the shear forces). The cross section of the monopile is usually constant until the transition piece, which might be conical in order to reduce settlements. Therefore, a potential failure is expected to take place at the bottom of the MP due to seismic loads along with the wind and wave loads. In this study, the wind turbine model is considered fixed-based at the sea-bed, thus the critical crosssection is considered to be the base of the monopile. However, it must be noted that the lack of redundancy for such a high-rise and slender structure may affect adversely the structural integrity of the entire infrastructure even in case of damage(s) located only in the monopile. Based on the above, it is necessary to calculate the cross sectional capacity of the monopile. The cross section of a wind turbine monopile is a thin wall cylinder of big diameter. If the cross section is classified according to EN 1993-1-1 [41], for tubular sections the limit to be within the first three cross section classes is 𝑑 ≤ 90𝜀 2 𝑡
(5.6)
44
where 𝜀 = √235/𝑓𝑦 , 𝑑 is the external diameter of the tubular section and 𝑡 represents the thickness of the wall. Therefore, for the dimensions specified in the subchapter 2.1 the cross section clearly falls within the class 4. The latter means that it is necessary to treat the monopile as a shell rather than a beam cross section and especially the provisions of EN 1993-1-6 [42] must be followed. Particularly, this Eurocode standard prescribes four different ULS scenarios that shall be taken into consideration:
LS1 - Plastic Limit. The failure if due to yielding of the material. The capacity of the cross section is reached. Hence the demand on a cross section 4 has to be derived only by elastic stress distribution as depicted in Figure 5-1.
Figure 5-1: Elastic stress distribution in a tubular cross section
LS2 - Cyclic Plasticity. The failure in this case is derived from repeated cyclic loadingunloading reaching the plastic limit in a specific point, which leads to excessive plastic work.
LS3 - Buckling. The failure mode in this case is driven by sudden displacements normal to the shell surface. The whole or part of the structure reaches its buckling capacity.
LS4 - Fatigue. Failure consists on propagation of small cracks generated by long term repeated cyclic loading, named fatigue cracks.
Despite the code considers four Limit States regarding the capacity of shell cross-sections, the maxima demand results, derived from the concurrent exposure of the wind turbine to the windwaves loads and the earthquake excitations, are of relevance for this study. Hence, LS1 and LS3 are adopted herein as the capacity ultimate limit states that will be compared with the hazards-driven demand. LS1 is related to the yielding capacity of the cross section. However, since the monopile’s cross section is classified into the forth category, exceeding the yielding point means plastic failure, since there is no further redundancy of the cross section to resist the imposed stresses after the yielding point. The latter explains why LS1 is associated with the Plastic Limit, according to the EN 1993-1-6 [42] prescriptions. Furthermore, the second 45
ultimate limit state, LS3, that was considered in this study is related to the stability of the monopile and more specifically to the shell buckling of the cross section. Several refined studies have already highlighted the relevance of the buckling failure mode and its complicated mechanism for the integrity and reliability of wind turbines [43]. According to the definition of cross section class 4 provided by EN 1993-1-1, buckling in terms of local instability is expected to happen before yielding. Hence, these two ultimate limit states (LS1 and LS3) for the capacity that were considered in this study correspond to two different damage states or better, different failure types, i.e., buckling failure mode, which is expected to happen before the yielding point, and the plastic limit, which is related to the exceedance of the yielding strength. In the following, the ultimate limit sates and the related failure modes are comprehensively described.
Before describing the limit states, it is necessary to mention that the different loading conditions, derived from the multi-hazard environment that the wind turbine is exposed to, create a complex stress field for the shell cross-sections. This stress field is presented in Figure 5-2.
Figure 5-2: Membrane stresses in shells, from [42] Regarding the LS1 evaluated in case of a shell cross-section, the membrane theory has to be applied for stress calculation. According to Annex A in EN 1993-1-6 [42], the following expressions can be used to define the stresses in unstiffened cylindrical shells:
46
Stress derived from uniform axial load (𝐹𝑥 ):
𝜎𝑋,𝑁 = −
𝐹𝑥 2𝜋𝑟𝑡
(5.7)
𝑀 𝜋𝑟 2 𝑡
(5.8)
Stress derived from bending moment (𝑀):
𝜎𝑋,𝑀 = ±
Stress derived from torsional moment (𝑀𝑡 ):
𝜏𝑡 =
𝑀𝑡 2𝜋𝑟 2 𝑡
(5.9)
Stress derived from shearing transverse force (𝑉):
𝜏𝑉 = ±
𝑉 𝜋𝑟𝑡
(5.10)
where 𝑟 is the radius of the monopile measured from the centreline and 𝑡 is the wall thickness. Furthermore, the aforementioned bending moment (𝑀) and shear force (𝑉) are defined as the square root of sum of squares in order to combine two directional results derived in the analysis. Meaning that the stresses are calculated based on resultant loads. Apart from that, the internal pressure and external pressure of the monopile’s cross section are assumed to be identical since the monopile is considered flooded. Therefore, no circumferential stresses are introduced into the section (𝜎𝜃 = 0). When different loads act in the cross section at the same time, EN 1993-1-6 [42] considers of good practice to represent the resultant two dimensional primary stresses field by equivalent von Mises stresses. This approach is shown in eq. (5.11).
2 2 2 𝜎𝑒𝑞 = √𝜎𝑋2 + 𝜎𝜃2 − 𝜎𝑋 · 𝜎𝜃 + 3(𝜏𝑋,𝜃 + 𝜏𝑋,𝑛 + 𝜏𝜃𝑛 )
(5.11)
47
where, 𝜎𝑋 = 𝜎𝑋,𝑁 + 𝜎𝑋,𝑀 and 𝜏𝑋.𝜃 = 𝜏𝑡 + 𝜏𝑉 for the current assessment. Therefore, the resulting plastic stress capacity can be rewritten in the following form.
2 𝜎𝑒𝑞 = √𝜎𝑋2 + 3 · 𝜏𝑋,𝜃
(5.12)
It should be noted that the expression derived from von Mises stress principle and shown in eq. (5.12) is partly conservative. This is due to the fact that the maximum stress produced by the shear and the maximum stress generated by the moment occur with 90 deg. phase lag within the tubular cross section. This can be seen from the elastic stress distribution in a tubular section shown in Figure 5-1, where the maximum stresses derived from moments coincide with zero tangential shear stresses and vice versa. Hence, in the current study, the LS1 was quantified through three different failure modes.
The demand in the cross section is derived for the locations of maximum tangential shear. At these points, the shear force is combined with the normal force and torsional moment.
As for the moment failure, it refers to the yielding of the cross section, where maximum bending moment takes place. In this case, the bending moment is superposed to the axial load and the torsional moment.
This last failure mode assumes that the stresses derived from tangential shear and bending moment are occurring at the same location. The latter is a conservative approach, still it is prescribed by the codes. Finally, regarding the LS1 analysis, the Utilization Ratio, shown in eq. (5.13), is considered as a representative parameter in order to evaluate the relation between the structural demand derived from the multi-hazard environment and the capacity of the cross section.
𝑼𝑹𝑳𝑺𝟏 =
𝝈𝒆𝒒 𝒇𝒚
(5.13)
48
where fy is the yield strength of the steel, used for the monopile. This representative value is obtained for each time step of the simulation. The resultant moment and shear are calculated following Pythagoras theorem and by means of the calculation of the membrane stresses the derivation of 𝑈𝑅 becomes possible.
The cross sectional stability is a key issue in case for shells and hence it must be thoroughly considered. The European Standard for Shell Strength and Stability offers different criteria to assess this problem. From the linear analysis (LA), which is evaluated based on simple stress calculations, to the complex geometrically and materially nonlinear analysis (GMNIA), which demands high time and computational cost [44]. The linear-bifurcation analysis and the materially nonlinear analysis (LBA/MNA) stand as an intermediate point. Along these lines, the conservativeness of the LA can be overcome and the complex solution, offered by GMNIA analysis, is avoided (only feasible to apply for really singular scenarios) achieving a final rather accurate result. This LBA&MNA method is described in the subchapter 8.6 of [42] and its application within the framework of this study is explained below. Buckling resistance is calculated from the amplification factor, 𝑟𝑅𝑑 which represents the capacity of the shell to withstand a specific loading condition. This factor is applied to the load case scenario (𝐹𝐸𝑑 ) that the cross section is facing and hence, the buckling resistance is represented by 𝐹𝑅 = 𝑟𝑅𝑑 · 𝐹𝐸 . Therefore, 𝑟𝑅𝑑 depicts the proportion that the load acting on the shell can be amplified before encountering buckling. Apart from that, this amplification factor has to be derived by means of the reference plastic and elastic resistances, 𝑟𝑅𝑝𝑙 and 𝑟𝑅𝑐𝑟 consequently. The first determines the utilization ratio of the cross sectional capacity as described in eq. (5.14), 𝑟𝑅𝑝𝑙 =
𝑡 · 𝑓𝑦 2 √𝑛𝑥2 − 𝑛𝑥 · 𝑛𝜃 + 𝑛𝜃2 + 𝑛𝑥𝜃
(5.14)
where the 𝑛𝑥 is the normal membrane stress, 𝑛𝜃 is the hoop membrane stress and 𝑛𝑥𝜃 is the shear membrane stress in the direction shown Figure 5-1. While the elastic resistance ratio, 𝑟𝑅𝑐𝑟 , is related to the geometry and it is determined from a Linear Bifurcation Analysis. The basis of this LBA procedure is an eigenvalue analysis which enables the determination of the lowest
49
eigenvalue, named bifurcation load factor. This lowest eigenvalue, dictates the factor that can be multiplied to the loads implemented for the LBA analysis before reaching the first mode shape or linear bifurcation scenario, commonly known as local buckling. Therefore, since the predominant load that may cause local buckling within the length of the OWT’s monopile is the bending moment, the LBA has been performed only considering the governing action (same applies for the determination of 𝑟𝑅𝑝𝑙 ). This assumption has been made since at the instant of maximum response of the OWT, the normal force derived from self-weight contributes to second order effects while shear stress does not coincide with bending moment stress (see Figure 5-1). Furthermore, the eigenvalue analysis has been conducted by means of the recognized FEM software Abaqus, where the applied load has been a unit moment. Consequently, the result from the LBA analysis has reflected the elastic critical moment, 𝑀𝑏𝑢𝑐𝑘𝑙𝑖𝑛𝑔 which has been used to determine the elastic critical resistance ratio as expressed in eq. (5.15).
𝑟𝑅𝑐𝑟,𝑖 =
𝑀𝑏𝑢𝑐𝑘𝑙𝑖𝑛𝑔 𝑀𝐸,𝑖
(5.15)
There, the 𝑀𝑏𝑢𝑐𝑘𝑙𝑖𝑛𝑔 refers to the critical elastic moment found from the LBA analysis for different can heights, 𝑀𝐸,𝑖 is the resultant moment for each time step of the simulations at the corresponding can height and 𝑟𝑅𝑐𝑟,𝑖 represents the elastic critical resistant ratio for a specific time step and a height along the cylinder. From the calculation of 𝑟𝑅𝑝𝑙 and 𝑟𝑅𝑐𝑟,𝑖 the derivation of the overall relative slenderness factor becomes possible. This way, the demand from the stresses considered by 𝑟𝑅𝑝𝑙 and the buckling demand reflected on 𝑟𝑅𝑐𝑟,𝑖 combined in eq. (5.16).
𝜆̅𝑜𝑣 = √𝐹𝑅𝑝𝑙 /𝐹𝑅𝑐𝑟 = √𝑟𝑅𝑝𝑙 /𝑟𝑅𝑐𝑟
(5.16)
As a result of the definition of the overall relative slenderness (𝜆̅𝑜𝑣 ), the concluding overall buckling reduction factor is determined as described by eq. (5.17),
50
1
𝑤ℎ𝑒𝑛
𝜆̅𝑜𝑣 ≤ 𝜆̅0𝑣,0
𝑤ℎ𝑒𝑛
𝜆̅𝑜𝑣,0 < 𝜆̅𝑜𝑣 < 𝜆̅𝑝
𝑤ℎ𝑒𝑛
𝜆̅𝑝 ≤ 𝜆̅𝑜𝑣
𝜂𝑜𝑣
𝜒𝑜𝑣
𝜆̅𝑜𝑣 − 𝜆̅0𝑣,0 1 − 𝛽 ( ) 𝑜𝑣 = 𝜆̅𝑝 − 𝜆̅0𝑣,0 𝛼𝑜𝑣 { 𝜆2̅𝑜𝑣
(5.17)
where 𝛼𝑜𝑣 is the overall elastic imperfection, 𝛽𝑜𝑣 is the plastic range factor, 𝜂𝑜𝑣 refers to the interaction exponent, 𝜆̅0𝑣,0 represents the squash limit relative slenderness and 𝜆̅𝑝 defines the 𝛼 plastic limit relative slenderness and is determined as 𝜆̅𝑝 = √1−𝛽𝑜𝑣 . These parameters are 𝑜𝑣
defined based on the Appendix D from EN 1993-1-6 [42] where the case of meridional compression represents the current buckling case. Additionally, the fabrication tolerance quality is assumed to be “normal” and the considered length of the cylinder that will be involved in the local buckling is assumed to be within the first six (6) meters. Moreover, the overall elastic imperfection has been derived in eq. (5.18),
𝛼𝑜𝑣 =
1
0.62 Δ𝑤 1.44 1 + 1.91 · ( 𝑡 𝑘 )
(5.18)
𝑟
where, Δ𝑤𝑘 = 𝑄 √𝑡 · 𝑡 and the value of 𝑄 is related to the fabrication quality which in this case has been conservatively assumed to be “normal”, 𝑄 = 16. Consequently, once the overall buckling reduction factor 𝜒𝑜𝑣 is defined, the overall buckling resistance ratio is derived in eq. (5.19).
(5.19)
𝑟𝑅𝑘 = 𝜒𝑜𝑣 · 𝑟𝑅𝑝𝑙
Therefore, buckling resistance verification is represented in eq. (5.20),
𝐹𝐸 ≤ 𝐹𝑅 = 𝑟𝑅 · 𝐹𝐸
𝑜𝑟
𝑟𝑅 ≥ 1
(5.20)
meaning that the utilization ratio in terms of buckling can be defined in the form of eq. (5.21),
51
𝑼𝑹𝑳𝑺𝟑
2 √𝑛𝑥2 − 𝑛𝑥 · 𝑛𝜃 + 𝑛𝜃2 + 𝑛𝑥𝜃 𝟏 = = 𝒓𝑹 𝜒0𝑣 · 𝑓𝑦 · 𝑡
(5.21)
where the capacity estimate for fragility assessment is driven by the yielding strength, 𝑓𝑦 . In addition, the location along the 30 meters of height of the monopile where the buckling demand should be considered in order to assess the cross section stability has been analysed in detail. For that purpose, the linear bifurcation analysis has been performed for different heights, covering the total length of the cylinder. This way, the buckling capacity along the height has been derived and considering that the maximum moment demand increases in depth, the critical region where buckling is most likely to happen has been identified.
Figure 5-3: Buckling capacity along the monopile (MP) height The result of this study is shown in Figure 5-3 where the buckling capacity is reflected as function of the height. Based on this figure, it can be concluded that the critical location where buckling is most likely to happen is around 5 m height. At this point the capacity converges with any further length meaning that the lowest capacity is reached. An example of the linear bifurcation analysis is shown in Figure 5-4. There, the first eigenvalue and mode shape are depicted. Additionally, boundary conditions have been always defined as a fixed base and a free end of the cylinder and the moment applied has been a unit [1N·mm]. Therefore, the bifurcation value (eigenvalue) represents the factor that can be multiplied to the applied load before buckling happens.
52
Eigen Value = 4.47684e+12
Figure 5-4: Linear Bifurcation Analysis (LBA) result for 5m high can and loading M=1Nmm Additionally, following the good practice recommendation given by EN 1993-1-6 [42], buckling assessment should be performed considering a minimum distance from the closest boundary condition defined in eq. (5.22). 0.1 · 𝐿 𝑙𝑅 = 𝑚𝑖𝑛 {
𝑟 𝑡
0.16 · 𝑟 · √
(5.22)
Where 𝐿 is the total length of the cylinder, 𝑟 is the radius of the shell from its centreline and 𝑡 represents the thickness of the shell wall. For the current monopile analysis, 𝑙𝑅 = 3.343𝑚. Hence, based on the sensitivity analysis shown in Figure 5-3 and code prescriptions [42], buckling assessment has been performed every meter from 3 to 6 m height of monopile as represented in Figure 6-4 with the purpose of defining the most critical spot.
Regarding the serviceability limit state (SLS), it was evaluated through this study on the basis of acceleration-related thresholds for the top of the tower (nacelle), where the damage in the highly-sensitive devices located there can cause serious operational problems and malfunctions for the wind turbines. Devises, like electrical generators, inverters and transformers as well as electrical controllers are sensitive to excessive accelerations, and hence it was deemed to investigate the fragility of these components to different acceleration thresholds, that have been identified as critical for these type of machinery. More specifically, the accelerations-based capacity of these devices have been already defined by Porter & Kiremidjian [45] as well as 53
Swan & Kassawara [46], who based their studies on post-earthquake assessment of the condition of different non-structural elements and components that can be found in industrial or residential buildings. The exact values for these acceleration thresholds and the related dispersion are given in the next section.
The two previous sections describe thoroughly all the demand-capacity considerations that were adopted in this study in order to be used for the fragility calculations. For ULS-LS1, the demand is calculated based on the stress field that is developed at the base of the monopile due to the concurrent application of the different loading actions (wind, waves, earthquakes) and the capacity is evaluated on the basis of yield strength of the steel (see the corresponding Utilization ratio shown in eq. (5.13)). For the ULS-LS3, the demand is based on the bending moment while the related capacity is defined according to the yield strength of the steel multiplied by the buckling reduction factor. Hence, it is required to define the median yield strength of the steel since this study is based on the fragility assessment power low for regression analysis which is recognized to be accurate for median demand estimate predictions as shown by Cornell et al. [38]. In the following, the way followed in this study to calculate the median value is thoroughly presented. Firstly, the mean value of the steel’s yield strength must be calculated. According to the Joint Committee for Structural Safety [47], the yield strength of the steel follows a lognormal distribution and its mean value can be estimated via the following expression:
𝑓̅𝑦 = 𝑓𝑦𝑠𝑝 · 𝛼 · 𝑒𝑥𝑝(−𝑢 · 𝑣) − 𝐶
(5.23)
where 𝑓𝑦𝑠𝑝 is the nominal value prescribed by the codes for specific plate thickness, 𝑣 is the related covariance, 𝑣 = 𝐶𝑂𝑉 = 0.07, 𝛼 is the spatial position factor (𝛼 = 1), 𝐶 is the reduction yield strength constant to get statistical data from tests results and 𝑢 is the factor, which relates the statistical proximity between the nominal value obtained from codes and the mean calculated based on eq. (5.23). This parameter always has to remain within the range (−1.5 𝑡𝑜 − 2.0). The value for this parameter has been chosen based on the literature, where experimental data has been the foundation for the decision. For this part of study, two different classes of steel were considered: S235 and S355. For the S235, Melcher et al. [48] showed that 54
the most appropriate value for u is equal to −2.0. For S355 steel, the research conducted by Zdeněk et al. [49] has been used as reference. In this second case, it has been found that the best match between the statistical results derived through experimental data and the calculated value following JCSS was achieved for 𝑢 = −2.0. Therefore, for both steel classes same value for the u factor has been applied. Consequently, the material data (steel’s yield strength characteristics) adopted herein are summarized in Table 5-1. Table 5-1: Nominal and mean yield strength for S235 and S355 𝑓𝑦𝑠𝑝 [𝑀𝑃𝑎]
𝑓̅𝑦 [𝑀𝑃𝑎]
S235
215
227.31
S355
335
365.34
Since the mean value of the steel’s yield strength has been already calculated, its median value follows. As mentioned above, in order to proceed with the calculation of the fragility or the probability of exceeding a certain damage state as it is defined in eq. (5.5), the median yield strength and its logarithmic standard deviation have to be defined. The median yield strength of a lognormally distributed case is calculated based on [50], as it is shown in eq. (5.24).
𝑓𝑦.𝑚𝑒𝑑𝑖𝑎𝑛 = 𝑒 𝜆
(5.24)
where 𝜆 is the standard deviation of the distribution natural logarithm, and can be obtained from the eq. (5.25), where the relationship between the mean and the mean & standard deviation of the distribution natural lograithm, 𝜁 and 𝜆 is shown.
𝑓̅𝑦 = 𝜇 = exp (𝜆 +
𝜁2 ) 2
(5.25)
where 𝜁 is obtained from the next expression, where the logarithmic standard deviation of the capacity is also defined.
55
(5.26)
𝜁 = 𝛽𝑐 = √𝑙𝑛(1 + 𝐶𝑂𝑉 2 ) = 0.0699
Hence, the capacity terms, assumed to be lognormally distributed, can defined by means of the following parameters as shown in Table 5-2. Table 5-2: Median yield strength and logarithmic standard deviation for S235 and S355 𝑓𝑦.𝑚𝑒𝑑𝑖𝑎𝑛 [𝑀𝑃𝑎]
𝛽𝑐
S235
226.75
0.0699
S355
364.45
0.0699
What has been described above concerns the capacity that is related to the ULS. On the other hand, the capacity for the SLS was defined according to acceleration thresholds for the devices that are located in the nacelle (see 5.2.2). More specifically, the thresholds for the acceleration and the related dispersion that were adopted based on the literature, are the following: electrical generators – 𝑎𝑐𝑐𝑆𝐿𝑆.𝑡ℎ𝑟 = 0.9𝑔 and 𝛽𝑐 = 0.50, inverters and transformers – 𝑎𝑐𝑐𝑆𝐿𝑆.𝑡ℎ𝑟 = 1.2𝑔 and 𝛽𝑐 = 0.60 as well as electrical controllers – 𝑎𝑐𝑐𝑆𝐿𝑆.𝑡ℎ𝑟 = 1.6𝑔 and 𝛽𝑐 = 0.80. To summarize all the above, Table 5-3 lists all the capacity-related parameters that were used for the fragility analysis of the offshore wind turbine. Table 5-3: Capacity estimate for structural and component reliability 𝑆𝑐
Limit state
𝛽𝑐
ULS1, S235
𝑓𝑦.𝑚𝑒𝑑.𝑆235
226.75
𝑀𝑃𝑎
0.0699
ULS1, S355
𝑓𝑦.𝑚𝑒𝑑.𝑆355
364.45
𝑀𝑃𝑎
0.0699
ULS3, S235
𝑓𝑦.𝑚𝑒𝑑.𝑆235
226.75
𝑀𝑃𝑎
0.0699
ULS3, S355
𝑓𝑦.𝑚𝑒𝑑.𝑆355
364.45
𝑀𝑃𝑎
0.0699
SLS, electrical generators
0.9 · 𝑔
8.82
𝑚/𝑠 2
0.50
SLS, inverters and transf.s
1.2 · 𝑔
11.76
𝑚/𝑠 2
0.60
SLS, electrical controllers
1.6 · 𝑔
15.68
𝑚/𝑠 2
0.80
56
The theoretical and probabilistic background for the fragility calculations and the demandcapacity considerations have been already described in the first two sections of the 5th chapter. The last component, which must also be described in order to close the loop for the necessary ingredients of a fragility analysis, is the loading environment that will be applied to the structural system studied. In this study, the 5MW offshore wind turbine is assumed to be exposed, at the same time, to three different hazardous sources: wind, waves and earthquake excitations. This is multi-hazard environment that several offshore wind turbines experience when thay have been installed in earthquake prone areas. This study intends to provide a generic assessment of the fragility for an offshore wind turbine. Hence, the selection of the loading components (wind, waves and erarthquake motions) do not correspond to a specific location of the world but it was intended to create such an exposures set, being representative of various sites. Therefore, some reasonable values for the wind and waves loading parameters were assumed, which are common for typical offshore conditions. Furthermore, special attention was paid to create a set of earthquake motions that includes strong ground motions with varying characteristics in terms of amplitude, frequency content, earthquake magnitude and epicentral distance. In this way, the fragility results will indirectly cover different cases of hazardous environment. Next, the parameters related to wind, waves are presented as well as a detailed description of the earthquake motions selected is provided.
Regarding the wind loads, they have been determined based on the mean wind speed, 𝑈𝑚 . Consequently, the turbulence intensity, corresponding to each 𝑈𝑚 , has been derived as explained in subchapter 2.3.1.1 according to IEC 61400-3 [29]. With the purpose to define the entire range of the operating window as well as the several wind speeds exceeding the operating range, different mean wind speeds were adopted for the time domain simuations and the consequent fragilty analysis (see Table 5-4).
57
Table 5-4: Mean wind speeds simulated Mean wind speed, 𝑈𝑚 [m/s] Operating
4
7
10
13
16
19
22
25
Parked
26
29
32
35
38
41
44
47
In this way, the operational case as well as the parked condition of the wind turbine can be studied.
Regarding the wave loads (sea-states) considered for the time-domain simulations and the consequent fragility analysis, different wave heights have been selected for the study. The basic assumption that an extreme hazard like the earthquake has a limited probability of happening jointly with a large waves derived from extreme storms. That is why the range of waves considered covers from 0.5 m wave height until a maximum of 6.5 m stepping every two meters leading to four different sea-sates. Furthermore, the current considered for simulation purposes, it has been stacked to 1m/s velocity. Ocean currents can vary from few centimetres to 4 meters per second. However, a common range would be 5m/s-50cm/s [51]. Therefore, as a conservative value it has been chosen to simulate 1m/s current as mentioned. It is also noted that both the wind and the wave loads were assumed to be aligned, meaning that they act in the same direction.
Time series of strong ground motions (in terms of accelerations) that have been recorded during past seismic events were used herein to represent the earthquake loads. The motions have been obtained by the PEER NGA Database, that currently includes more than 3551 publicly available, three-component strong ground motions (i.e., approximately 10,650 individual seismic records) that were recorded during 173 shallow crustal earthquakes from seismic prone regions around the world. As it is shown in Figure 5-5, the majority of the embedded seismic
58
events have occurred in California, US. Their magnitude, M, ranges from 4.2 to 7.9 and their epicentral distance, R, ranges from 0.20 to 600 km (Figure 5-6).
Figure 5-5: World map presenting the distribution of the 173 seismic events stored in the PEER-NGA database (distribution computed based on the information available online in http://peer.berkeley.edu/peer_ground_motion_database) Furthermore, the strong motion database provides information regarding the seismic source, including the date and time of the seismic event, the hypocentre’s location, the fault mechanism and additional, characteristics that reflect the overall seismotectonic environment. Data about 1600 recording stations (i.e., site characterizations and geological properties, shallow subsurface conditions and the exact location of the recording instrument) is also included in the PEER NGA-Database. Figure 5-5 presents the geographical distribution of 173 seismic events stored in the PEER-NGA database. Contrarily to the common practice, according to which the selection of earthquake motions is either based on a site-specific seismic hazard scenario or includes a limited number of earthquake records (close to seven motions) [52], the scope of the specific study and more specifically the fragility analysis purposes dictated the formation of a large set of motions with wealth of different characteristics. Especially, it is intended through this study to derive generic results for the fragility of the 5MW wind turbine and this requires using strong ground motions with a wide spectrum of characteristics in order to reflect the seismotectonic environment of different locations around the world. Therefore, a set of 40 pairs of horizontal components of strong ground motions were carefully selected in order to fulfil a multi-criteria framework, which is summarized in the next:
59
20 different seismic events recorded in different areas around the world were utilized to avoid dominating the set of motions with records from regions with common seismotectonic characteristics, the maximum number of seismic motions that are related to the same earthquake event was kept limited to 12.5% of the total, the preliminary selection criteria associated with seismological characteristics (i.e., earthquake magnitude, fault mechanism, directivity of seismic waves), amplitude (peak ground acceleration), frequency content and soil profile were deliberately kept wide enough, 20 pairs of horizontal components of strong motions were recorded during earthquake events of moderate magnitude, 6 25 km.
Figure 5-6: Earthquake magnitude as function of epicentral distance and peak ground acceleration of the strong ground motion records available in PEER-NGA database (distribution computed based on the information available online in http://peer.berkeley.edu/peer_ground_motion_database)
60
In accordance with aforementioned goals for earthquake motions selection, the set of earthquakes is presented in details in the Appendix A of this document. As an additional consideration, the directionality of the real ground motion has been studied for all 40 cases. By real ground motion, the combination of the X and Y direction of the motion pair in the horizontal plane is meant. Through this approach, the actual natural Peak Ground Acceleration (PGA) can be found from a trigonometric calculation where at each time step a resultant ground motion and the angle of attach is derived. The resultant picture where the maximum natural PGA in the corresponding direction for the 40 motions is shown in Figure 5-7.
Figure 5-7: Natural PGA magnitude and direction (around Z axis of the tower) As it has been described by means of eq. (5.1), the assessment for the structural reliability is performed regarding the magnitude which characterises the externally acting hazard. This magnitude is defined as Intensity Measure (IM) and its influence can be detrimental if the wrong IM is implemented for the probability study. Most commonly applied IM-s to define the magnitude and attributes of seismic hazards are two, Peak Ground Acceleration (𝑃𝐺𝐴) and Spectral Acceleration (𝑆𝑎 ). However, it is crucial to provide an IM that ensures a high correlation with the demand. Otherwise, the probabilistic model built to predict the fragility will not provide acceptable preditions. This is the reason why, the correlation of the Peak Ground Acceleration and the Spectral Acceleration with regard to the demand is studied at the beginning of results Chapter 6. In the following the definition and correlation of each of the potential IM-s are explained.
61
𝑃𝐺𝐴 There is more than one procedure when it comes to define the PGA. Motions in the horizontal plane of the earthquakes are depicted in two dimensional accelerograms. Pairs of motions represent the actual motion generated by the seismic event. That is the reason why the PGA can be defined independently of the actual motion (see eq. (5.27)), combining unidirectional peak ground accelerations or it can be determined based on the actual ground motion (see eq. (5.28)).
𝑃𝐺𝐴𝑔𝑚 = √𝑃𝐺𝐴𝑋 · 𝑃𝐺𝐴𝑌
𝑃𝐺𝐴𝑟𝑒𝑎𝑙 = 𝑚𝑎𝑥 (√𝐺𝐴2𝑥.𝑖 + 𝐺𝐴2𝑦.𝑖 )
(5.27)
(5.28)
In this two equations, 𝑃𝐺𝐴𝑔𝑚 is the geometrical mean between two independent 𝑃𝐺𝐴-s, 𝑃𝐺𝐴𝑋 & 𝑃𝐺𝐴𝑌 are two independent 𝑃𝐺𝐴-s obtained for each part of the pair of motions, 𝑃𝐺𝐴𝑟𝑒𝑎𝑙 is the actual Peak Ground Acceleration found for the combined accelerograms which is calculated as the maximum ground acceleration found along the duration of the earthquake and 𝐺𝐴𝑥.𝑖 & 𝐺𝐴𝑦.𝑖 are the independent ground accelerations found in each component of the pair of motions which are combined for each 𝑖 time-step. Some studies recommend to get rid of the influence of the directionality of the strong ground motions by means of the geometrical mean as described by Boore et al. [53]. This way, orientation-independent ground motions would be applied as Intensity Measure. In spite of that, based on the fact that the response of the wind turbine is assessed in terms of resultant response which combines two directional responses into one resultant vector, it is concluded that the most realistic way to consider PGA as an intensity measure would be by means of eq. (5.28). Through this procedure, the magnitude as well as the direction of the maximum or actual Peak Ground Motion can be derived. 𝑆𝑎 As far as the Spectral Acceleration is concerned for determination of the Intensity Measure, the first advantage is that it involves features of both, the earthquake as well as the structure. That makes it easier the find a higher correlation between the IM and the demand, especially when the structural predominant frequency is of long period. 62
The Spectral Acceleration is determined based on the elastic response spectra of each of the motions, where the spectral acceleration corresponding to the specific natural period of the first mode of the turbine can be obtained. For the case of the wind turbine, the difference of damping between the operating and the parked state is large, therefore, the 𝑆𝑎 has to be defined separately for each state (i.e., operational and parked conditions). In Figure 5-8, x and y response spectra for operating as well as parked turbine are shown where the 𝑆𝑎 -s are identified.
Figure 5-8: Elastic response spectra for operating (𝜁 = 5%) and parked (𝜁 = 1%) wind turbine [54]& [55], x & y directions Once the 𝑆𝑎 is determined for every case, the x and y values have to be combined. In this case the geometrical mean [53] is applied as shown in eq. (5.29).
2 · 𝑆2 𝑆𝑎 = √𝑆𝑎.𝑥 𝑎.𝑦
(5.29)
At the end, two sets of Spectral Accelerations are obtained, one for parked state and the other for operating state.
In order to determine the most appropriate Intensity Measure for the fragility assessment of the offshore wind turbine, the best correlated to the demand is chosen. In order to quantify the 63
correlation, the coefficient of determination is considered. This coefficient is based on the correlation coefficient which is defined in eq. (5.30) obtained as determined in [56],
𝑟𝑋𝑌 =
1 ∑𝑛𝑖=1(𝑥𝑖 − 𝑥̅ ) · (𝑦𝑖 − 𝑦̅) · 𝑛 𝑆𝑥 · 𝑆𝑦
(5.30)
where 𝑛 is the number of data points within each sample, 𝑥𝑖 & 𝑦𝑖 are the 𝐼𝑀 and the 𝑆𝑑 for each of the 𝑖 earthquakes, 𝑥̅ & 𝑦̅ refer to the mean value of the all the 𝐼𝑀 and the 𝑆𝑑 and 𝑆𝑥 & 𝑆𝑦 1
represent the sample standard deviation, 𝑆𝑥 = √ · ∑𝑛𝑖=1(𝑥𝑖 − 𝑥̅ )2. 𝑛
From eq. (5.30) the coefficient of determination is defined as expressed in eq. (5.31).
2 𝑅 2 = 𝑟𝑋𝑌
(5.31)
Hence, the comparison of the correlation in terms of coefficient of determination between the 𝑃𝐺𝐴 and 𝑆𝑎 is performed in Chapter 6, since data from the simulations is required for such a comparison.
Besides the selected 40 real ground motions, each of them has been scaled with a unique scaling factor, ranging from 2 up to 5, in order to get an homogeneously distributed sample of motions. Hence, all in all 16 different mean wind speeds are combined to 4 sea states and 80 earthquakes are simulated leading to an overall number of simulations of 5120. In the following the results derived from the fragility analysis of the OWT subjected to a fully coupled load scenario (comprising wind, waves, current and earthquakes) are brought to light.
64
Along this section of the report, the results related to the fragility assessment of a wind turbine subjected to aerodynamic, hydrodynamic and seismic loads is presented herein. It is noteworthy to mentioned that two approaches are presented. In the beginning three wind turbine states are studied according to the reference environmental load case scenarios presented in chapter 4.2, specifically in Table 4-1. Hence, for the first part of the analysis the sea state defined by a 6.5m high significant wave height, the suite of 80 motions and rated mean wind speed (for operating and operating & shutdown) as well as the 35 m/s wind speed (parked condition) have been studied in detail. On the other hand, the second part of this chapter considers all load case scenarios aforementioned in chapter 5.4, leading to a wide analysis which derives fragility surfaces.
Within this subchapter, different failure modes are studied according to the reference environmental load scenarios explained above. As defined by chapter 4.2, three different states of the wind turbine are analysed, parked condition, operating (no option for shutdown) and operating with a threshold of 1m/s2 at the nacelle which in case being exceeded a shutdown protocol is activated. The fragility analysis is divided into two parts. The first is related to the ULS that comprises the LS1 and LS3 according to EN 1993-1-6 [42]. While the second is regarding the SLS defined by the acceleration limits at the hub which have been defined in Table 5-3. Nevertheless, as a preliminary step, the most appropriate 𝐼𝑀 has been determined. For this purpose, parked wind turbine has been selected and based on the its response to the 80 strong ground motions, the determination coefficient has been derived for different failure modes with regard to two intensity measure (𝑃𝐺𝐴 & 𝑆𝑎 ). The results of this correlation comparison are shown in Table 6-1. From these results, it can be perceived that the shear failure mode is the only demand which has a higher correlation with the 𝑃𝐺𝐴. However, shear failure mode is of no influence in the integrity of the structure as it is shown below in subchapter 6.1.1. That is 65
why, according to the correlation comparison, it can be concluded that the most appropriate 𝐼𝑀 for an earthquake when the demand is related to a slender structure like the horizontal axis wind turbine is the Spectral Acceleration, 𝑆𝑎 . Table 6-1: Correlation comparison {𝐼𝑀(𝑃𝐺𝐴 𝑣𝑠. 𝑆𝑎 ) − 𝑆𝑑 } by means of the Coefficient of determination, 𝑅 2 𝑅 2 (coefficient of determination) IM
Shear
Moment
vonMises
Buck.S235
Buck.S235
Acc. hub hei.
PGA
0.83
0.17
0.18
0.11
0.11
0.28
Sa
0.19
0.67
0.67
0.72
0.72
0.56
Hence, the spectral acceleration described in chapter 5.3.3.2 is applied herein as the Intensity Measure within the framework of this study.
When the structural integrity of the wind turbine is considered, it refers to the reliability of the monopile which has to withstand the larges forces in terms of moment and shear. This part of the whole structure is crucial since is the only mean of support and connection to the ground. Therefore, in this study, the relation between the capacity of the monopile and the demand derived from the environmental and the seismic hazard is intended to be demonstrated.
As far as the ULS of the monopile is concerned, the EN 1993-1-6 is considered the reference for capacity and demand determination. Nevertheless, even though this code prescribes the stress calculation in a combined manner where normal, bending and shear loads are combined by means of vonMises criteria, it has been found of interest to perform at the same time a separate analysis in order to quantify the influence of each of the actions. That is why, in the following, shear failure mode, moment failure mode and the vonMises criteria failure mode are analysed separately.
66
The only difference between them is that when the resultant stress is derived, the contributing actions vary depending on the location of the cross section where the resultant stress is wanted to be calculated (see section 5.2.1). From Figure 5-1 it can be seen that the maximum stresses derived from moment and those derived from the shear occur at opposite locations within the cross section. That is why, it is possible a separate analysis of each failure mode.
Figure 6-1: Shear failure mode, PSDM In the Figure 6-1, the power model for shear mode failure is represented. It can be deduced that the correlation between shear (𝜎𝜏 ) and spectral acceleration (𝑆𝑎 ) is not very high due to the distribution of the scatter data of the demand. Additionally, the range of stresses where the maximum response of the eighty (80) motions lay are located is very low compared to the yielding limit, meaning that the structural integrity is not threatened by the shearing action.
Figure 6-2: Moment failure mode, PSDM However, for the case of the moment failure mode the correlation between the demand (𝜎𝑀 ) and the spectral acceleration (𝑆𝑎 ) is found higher than for the shear demand (see Figure 6-2). 67
Additionally, the peak stress level experienced by the cross section is considerable for this second case. In this occasion where the moment failure mode is studied, the most interesting case is the active shutdown simulation. It reflects that for some of the low intensity earthquakes the shutdown protocol is an effective method leading to a lower demand while for strong motions the performance of the turbine in maximum terms is very similar to the fully operational case. Additionally, the parked turbine shows how the difference in damping does not strictly impose larger peak response. From Figure 6-2 it can be seen that the parked wind turbine performance demand is lower than for the operating case for low intensity earthquakes while for higher intensity strong ground motions they are subjected to similar demand levels.
Figure 6-3: vonMises criteria’s failure mode, PSDM Hence, with regard to Figure 4-2 and based on Figure 6-2 it can come to an end that the response is significantly dependent on the frequency content and intensity of the seismic motion. However, the main different between parked and operational states resides on the damping ratio which is not significant with regard to peak response analysis. In order to derive higher responses due to lower damping, a number of oscillations is required where the motion keeps increasing during the earthquake excitation. This way the initial larger bending moment imposed by the wind to the operating turbine´s monopile (see chapter 4.2) can be exceeded by the parked case as a consequence of lower energy dissipated by the damping at every cycle (i.e., Manjil earthquake case in Figure 4-2). When all load actions in the cross section are combined by means of vonMises criteria as shows Figure 6-3, the conclusion is that the moment is the governing load since no significant difference is found between the moment failure mode and vonMises criteria failure mode. Furthermore, the deviation with regard to the median estimate increases slightly compared to 68
the moment failure mode. Meaning that the shearing action is contributing conservatively with regard to the prediction of the probability of failure. This conservativeness is inherent in the vonMises criteria, nevertheless, the actual difference can be assumed negligible as shows Figure 6-5.
When the buckling assessment of such a long cylinder is pursued, the location along the height where the local buckling or cross sectional instability will occur has to be identified. From the Linear Bifurcation Analysis (LBA) results presented in Figure 5-3, four spots were selected to be analysed. The outcome of such analysis is shown in Figure 6-4.
Figure 6-4: Buckling failure mode, PSDM (critical height identification) It can be seen from Figure 6-4 that no matter the increase on the buckling capacity for the shorter cans that the higher the demand, the higher will be the buckling fragility for these cases. Therefore, and based on eq. (5.22) where the minimum required distance from boundary conditions for buckling is defined (𝑙𝑅 = 3.343𝑚), the height of 4 meters is determined as critical spot and hence selected for analysis purposes. Note that the characteristics of the PSDM for buckling scenario at 4m height have been defined Table 6-2.
Finally, fragility curves for structural reliability have been derived based on the PSDM-s presented above and applying the eq. (5.5). In these CDF plots all three states of the wind turbine are presented for two different design options which are related to the applied material for the monopile production. 69
Figure 6-5: Fragility curves for structural reliability, CDF (ULS) It is clearly visible that for steel grade S235 the probability of exceeding the capacity is much higher for any of the analysed failure modes. On the other hand, shearing action shows relatively lower probability of failure which is product of the low correlation between the intensity measure and demand as well as the limited impact of the shear on the structural integrity in terms of stresses. Apart from that, von Mises approach is found to be almost equal to the case where moment demand is considered. Thus, vonMises approach does not represent a conservative approach for wind turbines monopile seismic assessment. On the other hand, buckling has been found to be the most demanding damage consideration. This is in accordance with the definition that is offered in EN 1993-1-1 [41] for the cross section class for which is expected to experience local buckling before yielding as it is shown by Figure 6-5. In addition, in order to interpret Figure 6-5, it has to be mentioned that the as a consideration of the different damping for the derivation of the spectral acceleration, the 𝐼𝑀 shown in the horizontal axis is not the same for the two operational cases and the parked case. Hence, as explained in chapter 5.3.3.2, different spectral accelerations are derived for each seismic event leading to a higher spectral acceleration for the lower damping of the parked wind turbine. That is why, since PSDM-s show a similar demand (𝑆𝑑 ) for all three turbine conditions, if a different 𝑆𝑎 is applied for regression analysis the fragility curve will be displaced to a region of higher
70
spectral acceleration values for the parked case. This means that the consideration of the fragility for parked and operating wind turbines has to be assessed separately.
In this section, the results regarding the serviceability limit state defined by accelerations at the hub are represented. This consideration is pertinent due to the sensitivity of different devices located at the hub (see chapter 5.2.2). In order to assess the probability of reaching the acceleration thresholds defined for these devices, probabilistic seismic demand models have been derived and they are plotted in Figure 6-6.
Figure 6-6: Acceleration at hub height, PSDM From Figure 6-6 the very little difference for different operational states is the predominant conclusion that can be derived. The response of the turbine derived from high instensity motions seem no to be affected by the shutdown mechanism while those in the lower magnitude area are more susceptible to experience modifications. However, due to the soft pitching applied during a long transient which enables the transition from operational to idling, the resulting acceleration and thus the motions at the hub do not suffer any dramatic alteration. Furthermore, regarding the relevance of the design of the applied shutdown protocol, it is worth to consider that the shutdown is usually designed for accelerations derived from gust wind impulses or too large irregularities in the torque. Nevertheless, when an earthquake hits the foundation of the turbine and the threshold acceleration is reached at hub height, the seismic wave will still be active leading to an undesirable scenario. The impulse generated by an emergency shutdown is combined to the quick loss of aerodynamic damping while the ground
71
motion still excites the structure. Therefore, the shutdown method for earthquake scenarios has to be progressively applied in order to avoid dramatic excitement of the wind turbine. Figure 6-7 reflects the CDF for the acceleration prediction at the top of the tower, which as mentioned before shows a similar estimate for the two operational conditions.
Figure 6-7: Fragility curves for component level reliability, CDF (SLS)
In this last subchapter of the Fragility Analysis, the results of the PSDM-s are summarized in tabulated format in the following Table 6-2. There, the fully operational, operational & shutdown plus the parked wind turbines PSDM-s key parameters are determined with the purpose to provide full information about the power model. However, the analysis is extended to different wind velocities as well as significant wave heights, leading to a wider scenario where the actual influence of environmental loads besides earthquake can be studied. This is shown in section 6.2 of the same chapter.
72
Table 6-2: Summary of PSDM results, power model definition for operation (OP), parked (PK) and operating & shutdown (SH) 𝑆𝑑 = 𝑎 · 𝐼𝑀𝑏 Limit State
Turbine state
Failure mode
𝑎
𝑏
𝛽𝐷|𝐼𝑀
𝑅2
18.629
0.507
0.610
0.187
21.154
0.565
0.580
0.270
SH
20.984
0.561
0.587
0.263
PK
98.588
0.666
0.389
0.675
128.589
0.673
0.356
0.682
SH
122.293
0.686
0.382
0.670
PK
99.714
0.664
0.393
0.666
129.813
0.671
0.361
0.673
SH
123.462
0.684
0.387
0.661
PK
0.404
0.710
0.379
0.719
0.545
0.710
0.347
0.721
SH
0.515
0.723
0.371
0.712
PK
0.275
0.710
0.379
0.719
0.371
0.710
0.347
0.721
SH
0.350
0.723
0.371
0.712
PK
3.595
0.686
0.460
0.558
4.586
0.704
0.534
0.543
4.722
0.683
0.512
0.533
PK OP
ULS -yielding-
OP
OP
OP ULS3 -buckling-
OP
SLS -hub acc.-
OP SH
Shear
Moment
vonMises
Buckling S235
Buckling S355
Acc. Hub height
73
Pursuing the initial objective of this MSc thesis, the seismic analysis of the offshore wind turbine subjected to the multi-hazard environemnt has been broadened to a large number of environmental load combinations. Consequently, a wide set of PSDMs have been derived that can be applied as preliminary reference cases when the seismic assessment of an OWT of similar characteristics is required. The purpose of the selection of the load case scenarios is to cover the whole range of operational mean wind speeds as well as the same extension beyond the maximum operational wind speed. For that reason, the load cases specified in section 5.3 have been studied. As a result, eight different operational wind speeds, eight parked wind states and four distinct sea-states have been simulated for the suite of 80 earthquakes. Nevertheless, unlike in the previous section 6.1, the shutdown simulation has not been included in the vast analysis. The main reasons why the shutdown simulation has not been taken into consideration is based on the fact that the maximum response has been found for the operational condition and hence does not represent a more detrimental operational condition than the constantly operational wind turbine. Moreover, the shutdown protocol can vary depending on the design which mainly considers the pitching velocities and time delays as well as transient periods when the pitching is applied. Based on the aforementioned reasons, the assessment of the shutdown for the probabilistic framework of this study has been omitted in this final part where reliable PSDM-s are aimed to be offered for further application. On the other hand, in the current offshore market, seldom structure is constructed whose primary steel is made of S235 grade steel. So this material is also disregarded for the main analysis. In addition, in order to avoid duplicated similar results, the main focus will be placed on one sea state (𝐻𝑠 = 6.5𝑚). However, all results are attached in Appendix B, where the required parameters to determine the PSDM for each load case are provided. Therefore, along the following subchapters the results derived by means of the same procedure developed in Chapter 5 are shown. The fragility analysis of the structure with regard to the Ultimate Limit State (yield & buckling) and Serviceability Limit State (acceleration at the hub) is shown herein.
74
Within the Ultimate Limit Estimate consideration, the results derived for fragility assessment for yielding of the cross section (ULS1) and local buckling of the shell (ULS3) are represented in the next chapters.
By means of the von Mises approach specified in EN 1993-1-6 [42] for cylindrical shells, the maximum stress calculation at the support of the monopile has been derived for each of the seismic events and specific load cases. Figure 6-8 provides the picture of the demand distribution in terms of maximum stresses for the whole set of Intensity Measures. From the general picture, it is evident the little affection of the wind speed on the demand for the operational OWT and the nearly inexistent influence on the parked case. In addition, it proves to be difficult to state a pattern relating the response of the parked and the operating wind turbines.
Figure 6-8: Stress demand in the cross section for wind speeds ranging from 4 to 47 m/s and Hs=6.5m This can be due to the fact that the frequency content, intensity and directionality of the seismic motions are influencing the response of the operating wind turbine which is experiencing an almost constant thrust from wind action. Consequently, due to the thrust and its inherent aerodynamic damping in one direction, the impact of the earthquake on the operating structure becomes hard to be related to the nearly independent effect of the earthquake on the parked turbine. 75
Additionally, significantly different probability prediction can be appreciated from Figure 6-9 where the fragility curves obtained based on the demand exposed in Figure 6-8 are visible. This dissimilarity is the result of the difference in the magnitude of the IM obtained from response spectrum derived considering the difference in damping between parked and operational turbines. The 𝑆𝑎 for the parked wind turbine has been obtained applying a damping ratio of 𝜁 = 1% while for the operational state 𝜁 = 5% has been utilized.
Figure 6-9: Cumulative Distribution for the probability of exceeding the yielding limit for wind speeds ranging from 4 to 47 m/s and Hs=6.5m Moreover, in line with the aforementioned predominant role of the earthquake, the fragility surface happens to be considerably independent from the influence of the mean wind speed. Despite that, a small plunge in the fragility surface around the cut-out wind speed can be perceived from Figure 6-9, which has been found for all sea-sates analysed (see Appendix B). Thus, the influence of the increased pitch angle for these higher operational wind speeds seem to be a reasonable explanation.
Comparably to the procedure followed for the yield limit state, the fragility surface for the buckling limit state has been derived. In this case the demand has been quantified based on the LBA approach prescribed by EN 1993-1-6 [42]. From that analysis, the Utilization Ratio has been determined for the location of 4 meters from the support which has been found to be the critical spot for cross section instability. The same line of actions has been followed for every single time domain simulation and the demand in terms of the UR for buckling as a function of the IM of the earthquake have been obtained as Figure 6-10 shows. 76
Figure 6-10: Buckling demand in terms of Utilization Ratio at a height of 4 m from the support for wind speeds ranging from 4 to 47 m/s and Hs=6.5m The distribution of the demand found in Figure 6-10 is very similar to the demand distribution shown above for the yielding limit. Hence, same conclusions apply for the case of the buckling as well. Nonetheless, as far as the fragility surface is concerned, it can be appreciated a higher probability estimate for the buckling limit state in Figure 6-11 compared to Figure 6-9 which is in agreement with the results shown in the previous section 6.1, where the prescriptions regarding local buckling of a cross section class 4 by the codes [42] has been corroborated.
Figure 6-11: Cumulative Distribution for the probability of exceeding the buckling limit for wind speeds ranging from 4 to 47 m/s and Hs=6.5m With the purpose of comparing the fragility estimate for the failure by reaching yielding and by local buckling of the monopile the most significant wind states have been selected. These notable cases are the cut-in, rated and cut-out wind speeds which are plotted in Figure 6-12 for both fragility estimates (LS1 & LS3). 77
Figure 6-12: Comparison of fragility curves for yielding limit state (LS1) and buckling limit state (LS3) according to the cut-in, rated and cut-out wind speeds The most perceptible idea that can be derived from Figure 6-12 is that the failure in local buckling will be precipitated earlier than the yield limit is reached within the cross section of the monopile. Therefore, the statement for the definition of cross section class 4 [41] is fulfilled in this analysis where the buckling failure mode happens to be the governing demand for the monopile. In addition to that, the rated wind speed can be pointed out as the most critical while the cut-out appears to be the least demanding. Therefore, from the perspective of the structural integrity it can be concluded that the worst case scenario for an operating wind turbine results at the rated wind speed in terms of local buckling demand of the cross section. Nevertheless, the probability of reaching the yield limit is considerable and hence must not be disregarded. In addition, the fragility surface for the parked condition has a very homogeneous distribution which is almost constant for varying mean wind speeds. Thus the comparison among the wind speeds results of limited relevance, nevertheless the overall fragility has to be considered with the same priority as for the operational case. It has to be considered that as it has been explained in section 6.1.1.3 of the report, the spectral acceleration specified as 𝐼𝑀 has a large influence on the difference between the shape of the fragility curve and hence surfaces for parked and operating wind turbines as a consequence of the damping implemented for the spectral acceleration derivation.
78
Within the probabilistic framework, the Serviceability Limit State determined based on component level fragility has been assessed with relation to the acceleration at the hub. This fragility estimate is based on the affection of the earthquake on the serviceability state of the infrastructure. As explained in the previous chapter 5.2.2 of this document, the sensitivity to accelerations of electrical generators, inverters & transformers as well as electrical controllers are intended to be considered from the serviceability perspective. In Figure 6-13 the demand of the accelerations at hub height are presented as function of the mean wind speed. Compared to the previous cases where the structural reliability has been analysed, the study of the accelerations shows a more homogeneous demand. Meaning that, the accelerations at the tower top are even more independent from the wind speed than the yield and buckling demands.
Figure 6-13: Acceleration demand at the hub for wind speeds ranging from 4 to 47 m/s and Hs=6.5m This higher independence with respect to the wind can be seen in Figure 6-14 where the fragility surface appears to be more homogeneous without any relevant sink or drop like has been identified for the LS1 and LS3 close to cut-out wind speed.
79
Figure 6-14: Cumulative Distribution for the probability of exceeding three significant acceleration thresholds at hub height for 4 to 47 m/s and Hs=6.5m Finally, in order to conclude with the extended fragility analysis from the serviceability perspective, the most important wind states have been compared in Figure 6-15. From this comparison, it can be identified that the lower is the wind speed the higher probability to exceed the acceleration capacity is derived. Therefore, the cut-in wind speed represents the worst case from the perspective of the SLS.
Figure 6-15: Comparison of fragility curves for acceleration limits at hub height (SLS) according to the cut-in, rated and cut-out wind speeds (SLS1 – Elec. generators, SLS2 – Inv.&trans. and SLS3 – elec. control.) As has been mentioned at the beginning of this chapter, four different sea states have been analysed within the framework of this MSc thesis. Details about the PSDM derived from the analysis are included in the Appendix B. Additionally, in the following section of this 80
document the influence of these four different sea-states on the response of the OWT analysed is studied.
With the purpose of getting some insight into the relevance of the sea-state considered for seismic assessment of wind turbines installed offshore, a comparative study is presented in the following. This subchapter intends to assist clarifying the influence that different sea-states determined according to the significant wave height have in demand terms when the OWT is subjected to an earthquake.
Figure 6-16: Influence of the sea states in the fragility estimate for ULS in yielding The fragility estimate for the yielding failure mode in ULS is shown in Figure 6-16. From this figure the extremely limited influence of sea condition on the response of the OWT is visible.
Figure 6-17: Influence of sea states in the fragility estimate for ULS in buckling
81
In total agreement with the findings regarding the influence of the significant wave height on the probability to failure in yielding at the base of the monopile, Figure 6-17 reflects the marginal influence for buckling failure.
Figure 6-18: Influence of sea states on the fragility estimate for SLS (ele. gen. 𝑆𝑐 = 0.9𝑔) Finally, the relevance of the sea-state for the accelerations found at the hub of an OWT when it is subjected to strong ground motions is represented by Figure 6-18. As it shows clearly, the influence is negligible for the SLS assessment as well. Hence, within the framework of this MSc thesis four different sea-states have been analysed in detail. Nevertheless, the influence of these sea states on the peak response (both in ULS as well as for SLS) of the OWT when it is subjected to an earthquake have been found to be marginal and thus not considerable.
82
In order to conclude this MSc thesis which has shown the possibility to run time domain analysis for earthquake assessment by means of HAWC2, final conclusions and potential considerations for further research purposes are explained in the following.
This vast probabilistic seismic assessment of a monopile based offshore wind turbine by means of time domain analysis performed through an aero-servo-elastic code (HAWC2) provides a wide set of probabilistic seismic demand models derived for different environmental load case scenarios based on 80 seismic motion pairs. In addition, a refined procedure has been implemented for the selection of the forty strong ground motions applied during the research which increases the reliability of the study by reducing the inherent uncertainty of earthquake’s seismic characteristics. Hence, this MSc thesis makes available pertinent PSDMs which can be utilized as generic reference for the seismic response study of comparable OWTs. Throughout the research conducted for the fragility assessment of an OWT, several concluding statements have been derived, which are concisely summarized in the upcoming. The earthquake has been identified as the governing hazard for the wind turbine response both in terms of stresses at the monopile base and accelerations at the hub. This leads to the fact that the seismic hazard may represents the predominant design load case if the OWT is planned to be installed in an earthquake prone area. The influence of the mean wind speed for the response of the wind turbine subjected at the same time to the seismic action has been found to be of limited significance. Nevertheless, the rated wind speed was found to be the most critical in terms of the fragility for the monopile while the cut-in wind speed was the most detrimental for the accelerations at hub height. Furthermore, the relevance of the significant wave height is determined to be marginal for the fragility assessment. Meaning that the influence of the sea-state for fragility
83
analysis of an offshore wind turbine subjected at the same time to an earthquake is negligible based on the findings of this study. On the other hand, considerably high probability (i.e., 0.5 to 0.7) for the demand to exceed the capacity has been found through the study even for moderate earthquake intensity (𝑆𝑎 ⋍ 0.5𝑔). This means, that the earthquake hazard should be considered for the analysis, design and/or assessment of OWT for seismically active areas. For the vast majority of the earthquake strong ground motions used herein, the operational state of the wind turbine has been found to be the case where the highest demands were calculated for the bottom of the monopile. This means that it is quite important to be able to capture the dynamic response of an OWT during its operational status. This can be only done using aero-elastic codes like HAWC2, which at the same time allow for earthquake analysis in a fully coupled environment with the rest concurrent load conditions (wind, waves & currents).
As a consequence of the time constrain and the need to orient the study in a specific direction, some probably interesting aspects within the framework of the fragility analysis of offshore wind turbines were disregarded. Nevertheless, in this last section of the document, some further steps for future implementation are brought to light. Along this study performed as a MSc thesis, the structure has been assumed to be fixed at the bottom of the monopile. Hence, soil structure interaction has been neglected and consequently the influence of the soil flexibility on the total dynamic response of the OWT has not been captured within this work. Apart from that, determination of the effect of the directionality of the seismic waves with regard to the orientation of the wind turbine can be assumed to be relevant. Even more when the study is design oriented and the location of the turbines is known. At the same time, further analysis aiming to define the influence of the frequency content of the seismic motions on the dynamic response of a wind turbine could derive interesting results. In addition, it has to be considered that the central point of this MSc thesis is the prediction of the wind turbine’s performance when it is exposed to multi-hazard environment consisting of wind, waves and earthquakes. The set of the loading conditions was chosen deliberately wide 84
in order to derive fragilities for a broad spectrum of different wind, wave and earthquake conditions. However, it could be interesting to consider the definition of the loading parameters based on a probabilistic approach and to specify accordingly the loading parameters for a particular region considering a life-time framework. This means that site-specific data about the seismicity and the wind/wave environment can be convoluted in a probabilistic way with demand/capacity models and hence the risk for a specific wind turbine located in a particular region could be estimated. Moreover, aging and deterioration of the infrastructure as well as the fatigue consideration could be involved in a study if the lifetime reliability of the structure was the objective.
85
[1]
U.S. Energy Information Administration, "International Energy Statistics - EIA," [Online].
Available:
https://www.eia.gov/cfapps/ipdbproject/IEDIndex3.cfm?tid=2
&pid=2&aid=12. [Accessed February 2016]. [2]
G. Bredehoeft, "AEO2014 - Issues in Focus Articles - U.S. Energy Information Administration," 29 April 2014. [Online]. Available: http://www.eia.gov/forecasts/ archive/aeo14/section_issues.cfm#elec_proj. [Accessed March 2016].
[3]
J. MacDonald, “Bloomberg New Energy Finance,” 14 January 2016. [Online]. Available:
http://about.bnef.com/press-releases/clean-energy-defies-fossil-fuel-price-
crash-to-attract-record-329bn-global-investment-in-2015/. [Accessed 2016 February]. [4]
Global Wind Energy Council (GWEC), “Global Wind Report - Annual Market Update 2015,” [Online]. Available: http://www.gwec.net/wp-content/uploads/vip/GWECGlobal-Wind-2015-Report_April-2016_22_04.pdf. [Accessed June 2016].
[5]
N. Bazeos, G. Hatzigeorgiou, I. Hondros , H. Karamaneas , D. Karabalis and D. Beskos, "Statis, seismic and stability analysis of a prototype wind turbine steel tower," Engineering Structures, pp. 1015-1025, 2002.
[6]
U. Ritschel, I. Warnke, J. Kirchner and B. Meussen, "Wind turbines and earthquakes," in Proceedings of the 2nd World Wind Energy Conference, Cape Town, South Africa, 2003.
[7]
D. Witcher, “Seismic analysis of wind turbines in the time domain,” Wind Energy, pp. 81-91, 2005.
[8]
L. Wang and Y. Zhang , “Influence of simplified models on seismic response analysis of wind turbine towers,” Applied Mechanics and Materials, pp. 369-374, 2011.
86
[9]
B. Song, Y. Yi and J. Wu, “Study on seismic dynamic response of offshore wind turbine tower with monopile foundation based on M method,” Advanced Materials Research, pp. 686-691, 2013.
[10] G. He and J. Li, “Seismic analysis of wind turbine system including soil-structure interaction,” in Proceedings of the 14th World Conference on Earthquake Engineering, Beijin, China, 2008. [11] J. Jonkman and W. Musial, “Offshore Code Comparison Collaboration (OC3) for IEA Task 23 Offshore Wind Technology and Deployment,” NREL, National Renewable Energy Laboratory, Golden, Colorado, 2010. [12] M.-A. Asareh and J. Volz, “Evaluation of aerodynamic and seismic coupling for turbines using finite element approach,” in Proceedings of the ASME 2013 International Mechanical Engineering Congress and Exposition, San Diego, CA, USA, 2013. [13] N. Alati, G. Failla and F. Arena, “Seismic analysis of offshore wind turbines on bottomfixed support structures,” The Royal Society, University of Reggio Calabria, 2015. [14] I. Prowell, “An experimental and numerical study of wind turbine seismic behaviour,” San Diego, CA, USA, 2011. [15] X. Jin , H. Liu and W. Ju, “Wind turbine seismic load analysis based on numerical calculation,” Slovenian Journal of Mechanical Engineering, pp. 638-648, 2014. [16] O. Kiyomiya, T. Rikiji and P. van Gelder, “Dynamic response analysis of onshore wind energy power units during earthquakes and wind,” in Proceedings of the 12th International Offshore and Polar Engineering Conference, Kitakyushu, Japan, 2002. [17] E. Nuta, C. Christopoulos and J. Packer, “Methodology for seismic risk assessment for tubular steel wind turbine towers: application to Canadian seismic environmnet,” Canadian Journal of Civil Engineering, pp. 293-304, 2011. [18] A. Mensah, L. Duenas-Osorio, I. Prowell and M. Asareh, “Probabilistic combination of earthquake and operational loads for wind turbines,” in Proceedings of the 15th World Conference of Earthquake Engineering, Lisboa, Portugal, 2012.
87
[19] M. Mardfekri and P. Gardoni, “Multi-hazard reliability assessment of offshore wind turbines,” Wind Energy, pp. 1433-1450, 2014. [20] IEC - International Electrotechnical Commission, “Wind turbines - Part 1: Design requirements,” Geneva, Switzerland, 2005. [21] Det Norske Veritas & Germanischer Lloyd (DNV-GL), Standard DNVGL-ST-0126: Support structures for wind turbines, 2016. [22] European Commitee for Standards, CEN, Eurocode 8: Design of structures for earthquake resistance - Part 1: General rules, seismic action and rules, 2004. [23] International Organization for Standarization (ISO), Petroleum and natural gas industries: Specific requirements for offshore structures - Part 2: Seismic design procedures and criteria, 2015. [24] International Building Code (IBC), International Code Council, Country Club Hills, IL, USA, 2006. [25] American Society of Civil Engineers, ASCE/SEI 7-10 Minimum design loads for buildings and other structures, Reston, VA, USA, 2010. [26] T. J. Larsen and A. H. Hansen, “How 2 HAWC2, the user's manual,” Risø National Laboratory, Technical University of Denmark, Roskilde, Denmark, 2015. [27] J. Jonkman, S. Butterfield, W. Musial and G. Scott, “Definition of a 5-MW Reference Wind Turbine for Offshore System Development,” National Renewable Energy Laboratory (NREL), Golden, Colorado, USA, 2009. [28] O. A. Bauchau, Flexible Multibody Dynamics, Atlanta, Georgia: Springer, 2011. [29] IEC - International Electrotechnical Commission, “Wind turbines - Part 3: Design requirements for offshore wind turbines,” in IEC 61400-3, Geneva, Switzerland, 2009. [30] T. Burton, N. Jenkins, D. Sharpe and E. Bossanyi, Wind Energy Handbook, chapter 3, Chichester, United Kingdom: John Wiley & Sons, Ltd, 2011.
88
[31] WMO - World Meteorological Organization, “Guide to Wave Analysis and Forecasting,” 1998, pp. 6-14. [32] Z. Liu and P. Frigaard, Generation and Analysis of Random Waves, Aalborg: Lab oratoriet for Hydraulik og Havnebygning; Instituttet for Vand, Jord og Miljøteknik; Aalb org Unive rsitet, 2001. [33] B. M. Sumaer and J. Fredsøe, Advanced Series on Ocean Engineering - Volume 26. Hydrodynamics around Cylindrical Structures, Singapore: World Scientific, 2006. [34] A. K. Chopra, Dynamics of Structures, Theory and Applications to Earthquake Engineering, University of California at Berkeley: Prentice Hall, 2012. [35] O. Taskari and A. Sextos, “Multi-angle, multi-damage fragility curves for sesimic assessment of bridges,” Earthquake Engineering and Structuctural Dynamics, pp. 22812301, 2015. [36] O. Taskari and A. Sextos, “Probabilistic Assessment of Abutment-Embankment Stiffness and Implications in the Predicted Performeance of Short Bridges,” Journal of Earthquake Engineering, pp. 822-846, 2015. [37] L. E. Yamins, A. I. Hurtado, A. H. Barbat and O. D. Cardona, “Seismic and wind vulnerability assessment for the GAR-13 global risk assessment,” International Journal of Disaster Risk Reduction, pp. 452-460, 2014. [38] A. Cornell, F. Jalayer and R. Hamburger, “Probabilistic basis for 2000 SAC federal emergency management agency steel moment frame guidelines,” Journal of Structural Engineering, no. 128, pp. 128: 526-532, 2002. [39] B. G. Nielson and R. DesRoches, “Seismic fragility methodology for highway bridges using a component level approach,” Earthquake Engineering and Structuctural Dynamics, pp. 823-839, 2006. [40] J. E. Padgett and R. DesRoches, “Methodology for development of analytical fragility curves for retrofitted bridges,” Earthquake Engineering and Structuctural Dynamics, pp. 1157-1174, 2008.
89
[41] European Commitee for Standards, CEN, Eurocode 3: Design of steel structures - Part 1-1: General Rules and Rules for Buildings, 2005. [42] European Commitee for Standards, CEN, Eurocode 3: Design of steel structures - Part 1-6: Strength and Stability of Shell Structures, 2007. [43] S. Thöns, M. H. Faber and W. Rücker, “Ultimate Limit State Model Basis for Assessment of Offshore Wind Energy Converters,” Journal of Offshore Mechanics and Arctic Engineering, vol. 134, no. 031904, 2012. [44] M. J. Rotter, “Shell buckling design and assessment and the LBA-MNA methodology,” Ernst & Sohn, pp. 791-803, 2011. [45] K. Porter and A. Kiremidjian, “Assembly-based vulnerability of buildings and its uses in seismic performance evaluation and risk management decision-making,” The John A, Blume Earthquake Engineering Center, Stanford University, Stanford, 2001. [46] S. W. Swan and R. Kassawara, “The Use of Earthquake Experience Data for Estimates of the Seismic Fragility of Standard Industrial Equipment,” in ATC-29-1, Proc., Seminar on Seismic Design, Retrofit, and Performance of Nonstructural Components, Redwood City CA: Applied Technology Council, 1998. [47] Joint Comitee on Structural Safety, Technical University of Denamrk, “JCSS Probabilistic Model Code, Part 3: Material Properties,” 2000. [48] . J. Melcher, Z. Kala, M. Holický, M. Fajkus and L. Rozlívka, “Design characteristics of structural steels based on statistical analysis of metallurgical products,” Journal of Constructional Steel Research, pp. 795-808, 2004. [49] K. Zdeněk, J. Melcher and L. Puklický, “Material and Geometrical Characteristics of Structural Steel,” Journal of Civil Engineering and Management, pp. 299-307, 2009. [50] D. Kececioglu, Reliability Engineering Handbook, vol. 1, Englewood Cliffs, New Jersey: Prentice Hall, Inc., 1991.
90
[51] Enciclopædia Britanica, School and Library Subscribers, “ocean current | Britannica,” [Online]. Available: http://global.britannica.com/science/ocean-current. [Accessed 03 2016]. [52] E. Katsanos, A. Sextos and G. Manolis, “Selection of earthquake ground motion records: A state-of-the-art review from a structural engineering perspective,” Soil Dynamics and Earthquake Engineering, pp. 157-169, 2010. [53] D. M. Boore, J. Watson-Lamprey and N. A. Abrahamson, “Orientation-Independent Measures of Ground Motion,” Bulletin of Seismological Society of America, vol. 96, no. 4A, pp. 1502-1511, 2006. [54] G. Stamatopoulus, “Response of a wind turbine subjected to near-fault excitation and comparison with Greek seismic code provisions,” Soil Dynamic and Earthquake Engineering, pp. 77-84, 2013. [55] I. Prowell, M. Veletzos, A. Elgamal and J. Restrepo, “Experimental and numerical seismic response of a 65 kW wind turbine,” Journal of Earthquake Engineering, pp. 1172-1190, 2009. [56] M. H. Faber, Statistics and Probability Theory: In Persuit of Engineering Decision Support, Springer, 2012.
91
92
1979
1979
1979
1980
1985
1986
1986
1987
1987
1987
1989
1989
1992
1992
1992
Imperial Valley-06
Imperial Valley-06
Imperial Valley-06
Irpinia, Italy01
Nahanni, Canada
Taiwan SMART1(45)
Taiwan SMART1(45)
Superstition Hills-02
Superstition Hills-02
Superstition Hills-02
Loma Prieta
Loma Prieta
Erzican, Turkey
Cape Mendocino
Cape Mendocino
Landers
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
Northridge01
1979
Imperial Valley-06
3
20
1978
Tabas, Iran
2
1994
1992
1976
Gazli,USSR
1
Year
Earth. name
Id#
117
628
425
425
313
1018
1018
1124
1124
1124
1114
1114
1223
1123
1015
1015
1015
1015
916
517
MODY
6.69
7.28
7.01
7.01
6.69
6.93
6.93
6.54
6.54
6.54
7.3
7.3
6.76
6.9
6.53
6.53
6.53
6.53
7.35
6.8
Magnit .
27.47 28.09 19.81
El Centro Array #6 El Centro Array #8 Holtville Post Office
Jensen Filter Plant Generator
Lucerne
Petrolia
Cape Mendocino
13
44.02
4.51
10.36
8.97
42.13
Sunnyvale Colton Ave. Erzincan
18.46
29.41
Wildlife Liquef. Array LGPC
11.2
15.99
Parachute Test Site Poe Road (temp)
76.01
76.21
6.8
SMART1 I07
SMART1 C00
Site 1
30.35
19.44
EC Meloland Overpass FF
Sturno
55.24
12.82
EpiD [km]
Tabas
Karakyr
Record. station
C
C
C
C
D
D
C
D
D
D
D
D
C
B
D
D
D
D
B
C
NEHRP
525.8
684.9
712.8
513.7
274.5
267.7
477.7
207.5
207.5
348.7
274.5
274.5
659.6
1000
202.9
206.1
203.2
186.2
766.8
659.6
Vs30[m/s]
0.571/1.024
0.727/0.789
0.590/0.662
1.497/1.039
0.515/0.496
0.207/0.209
0.966/0.587
0.181/0.202
0.446/0.300
0.455/0.377
0.118/0.122
0.122/0.153
0.978/1.096
0.251/0.358
0.253/0.221
0.602/0.145
0.410/0.366
0.317/0.296
0.836/0.852
0.609/0.716
PGA (x/y) [g]
0.2440.080
0.327/0.112
0.070/0.142
0.192/0.078
0.167/0.176
0.168/0.183
0.424/0.093
0.094/0.198
0.064/0.080
0.2320.120
0.110/0.091
0.106/0.096
0.089/0.107
0.092/0.209
0.160/0.168
0.113/0.003
0.204/0.098
0.148/0.298
0.203/0.334
0.191/0.169
Sa,3.6s(x/y)[g]
0.791/0.575
0.314/0.165
0.540/0.677
0.363/0.352
1.384/0.839
1.539/1.399
0.802/0.564
0.896/1.159
0.484/0.478
1.110/0.663
1.102/1.095
0.943/0.933
0.283/0.357
0.676/0.862
0.618/0.704
0.612/0.364
0.939/0.421
1.048/1.768
0.483/0.456
0.431/0.403
Tm (x/y) [s] [g]
93
Year 1994 1994 1994 1995 1995 1995 1999 1999 1999 1999 1999 1999 1999 1999 1979 1979 1990 1990 1999 1999
Earth. name
Northridge01
Northridge01
Northridge01
Kobe, Japan
Kobe, Japan
Dinar, Turkey
Kocaeli, Turkey
Kocaeli, Turkey
Chi-Chi, Taiwan
Chi-Chi, Taiwan
Chi-Chi, Taiwan
Chi-Chi, Taiwan
Duzce, Turkey
Duzce, Turkey
St Elias, Alaska
St Elias, Alaska
Manjil, Iran
Manjil, Iran
Hector Mine
Chi-Chi, Taiwan
Id#
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
920
1016
620
620
228
228
1112
1112
920
920
920
920
817
817
1001
116
116
117
117
117
MODY
7.62
7.13
7.37
7.37
7.54
7.54
7.14
7.14
7.62
7.62
7.62
7.62
7.51
7.51
6.4
6.9
6.9
6.69
6.69
6.69
Magnit .
TCU065
Amboy
Abhar
Abbar
Yakutat
Icy Bay
Duzce
Bolu
TCU075
TCU068
TCU052
CHY101
Yarimca
Duzce
Dinar
Takatori
KJMA
SylmarConverter Sta
Newhall - W Pico Canyon Rd.
Newhall - Fire Sta
Record. station
26.67
47.97
77.84
40.43
159.99
74.84
1.61
41.27
20.67
47.86
39.58
31.96
19.3
98.22
0.44
13.12
18.27
13.11
21.55
20.27
EpiD [km]
D
D
D
C
D
D
D
D
C
C
C
D
D
D
D
D
D
D
D
D
NEHRP
305.9
271.4
274.5
724
274.5
274.5
276
326
573
487.3
579.1
258.9
297
276
219.8
256
312
251.2
285.9
269.1
Vs30[m/s]
0.814/0.603
0.182/0.150
0.132/0.209
0.515/0.496
0.083/0.065
0.099/0.176
0.348/0.535
0.728/0.822
0.333/0.264
0.566/0.462
0.348/0.419
0.353/0.440
0.268/0.349
0.312/0.358
0.352/0.282
0.611/0.616
0.821/0.599
0.612/0.897
0.455/0.325
0.583/0.590
PGA (x/y) [g]
0.383/0.208
0.117/0.069
0.146/0.479
0.093/0.112
0.158/0.166
0.049/0.093
0.236/0.221
0.106/0.064
0.303/0.135
0.438/0.247
0.366/0.252
0.311/0.310
0.310/0.296
0.236/0.085
0.088/0.067
0.217/0.226
0.105/0.117
0.345/0.119
0.258/0.103
0.114/0.165
Sa,3.6s(x/y)[g]
1.130/1.029
0.781/0.629
0.744/1.514
0.316/0.321
1.904/2.017
1.097/1.436
0.693/0.818
0.547/0.785
0.776/0.664
1.511/1.286
1.584/1.528
1.054/0.964
1.242/1.341
0.991/0.876
0.892/0.683
1.131/0.986
0.646/0.656
1.158/0.994
1.586/1.249
0.525/0.697
Tm (x/y) [s] [g]
Along this appendix, the parameters required to define the Probabilistic Seismic Demand Models are provided for all environmental load case scenarios considered in this seismic analysis of wind turbines. In addition, the determination coefficient is offered for each of the PSDM derived from the 80 motions. Different load cases are classified based on the sea-state or significant wave height applied for the simulation.
94
𝑯𝒔 = 𝟎. 𝟓𝒎
Cumulative distribution of probability of failure in LS1, Hs=0.5m 𝑺𝒅,𝒚𝒊𝒆𝒍𝒅 = 𝒂 · 𝑰𝑴𝒃
𝐻𝑠 = 0.5𝑚 𝑈𝑚 [m/s]
𝒂
𝒃
𝜷
𝑹𝟐
4
117.895
0.717
0.394
0.682
7
114.807
0.732
0.370
0.700
10
130.486
0.665
0.347
0.691
13
128.201
0.675
0.356
0.682
16
123.948
0.684
0.369
0.672
19
120.849
0.692
0.374
0.671
22
113.296
0.694
0.400
0.650
25
127.170
0.629
0.369
0.638
26
96.326
0.685
0.397
0.681
29
96.599
0.680
0.395
0.676
32
97.028
0.675
0.392
0.676
35
97.770
0.669
0.390
0.675
38
98.840
0.662
0.386
0.673
41
100.096
0.655
0.383
0.669
44
101.610
0.647
0.379
0.667
47
103.456
0.638
0.376
0.663 95
Cumulative distribution of probability of failure in LS3, Hs=0.5m
𝑺𝒅,𝒃𝒖𝒄𝒌𝒍𝒆 = 𝒂 · 𝑰𝑴𝒃
𝐻𝑠 = 0.5𝑚 𝑈𝑚 [m/s]
𝒂
𝒃
𝜷
𝑹𝟐
4
0.336
0.758
0.385
0.721
7
0.329
0.771
0.357
0.740
10
0.378
0.698
0.334
0.732
13
0.370
0.709
0.343
0.723
16
0.356
0.720
0.356
0.714
19
0.345
0.731
0.363
0.713
22
0.321
0.735
0.388
0.695
25
0.372
0.651
0.356
0.681
26
0.268
0.730
0.387
0.725
29
0.269
0.725
0.384
0.721
32
0.271
0.719
0.381
0.720
35
0.274
0.712
0.378
0.719
38
0.278
0.703
0.374
0.717
41
0.282
0.695
0.370
0.715
44
0.287
0.685
0.366
0.714
47
0.293
0.674
0.363
0.711
96
Cumulative distribution function of probability of exceeding three different acceleration thresholds, Hs=0.5m
𝑺𝒅,𝒂𝒄𝒄𝒆𝒍. = 𝒂 · 𝑰𝑴𝒃
𝐻𝑠 = 0.5𝑚 𝑈𝑚 [m/s]
𝒂
𝒃
𝜷
𝑹𝟐
4
4.643
0.706
0.490
0.542
7
4.477
0.724
0.472
0.554
10
4.592
0.703
0.486
0.544
13
4.574
0.706
0.492
0.544
16
4.559
0.702
0.496
0.534
19
4.515
0.705
0.498
0.533
22
4.446
0.690
0.501
0.515
25
4.570
0.676
0.491
0.512
26
3.623
0.692
0.455
0.568
29
3.613
0.690
0.458
0.561
32
3.606
0.688
0.458
0.560
35
3.599
0.686
0.459
0.558
38
3.592
0.683
0.459
0.560
41
3.589
0.682
0.462
0.550
44
3.586
0.678
0.460
0.554
47
3.589
0.676
0.460
0.552
97
𝑯𝒔 = 𝟐. 𝟓𝒎
Cumulative distribution of probability of failure in LS1, Hs=2.5m 𝑺𝒅,𝒚𝒊𝒆𝒍𝒅 = 𝒂 · 𝑰𝑴𝒃
𝐻𝑠 = 2.5𝑚 𝑈𝑚 [m/s]
𝒂
𝒃
𝜷
𝑹𝟐
4
118.074
0.716
0.394
0.682
7
115.017
0.731
0.370
0.701
10
130.794
0.664
0.346
0.690
13
128.229
0.674
0.355
0.682
16
123.938
0.683
0.369
0.672
19
120.948
0.692
0.373
0.671
22
112.998
0.696
0.400
0.652
25
127.669
0.627
0.366
0.638
26
96.719
0.684
0.396
0.682
29
96.484
0.680
0.397
0.675
32
96.936
0.676
0.393
0.676
35
98.210
0.668
0.388
0.675
38
96.865
0.671
0.390
0.674
41
101.001
0.652
0.379
0.670
44
102.095
0.645
0.376
0.667
47
103.319
0.639
0.378
0.663 98
Cumulative distribution of probability of failure in LS3, Hs=2.5m
𝑺𝒅,𝒃𝒖𝒄𝒌𝒍𝒆 = 𝒂 · 𝑰𝑴𝒃
𝐻𝑠 = 2.5𝑚 𝑈𝑚 [m/s]
𝒂
𝒃
𝜷
𝑹𝟐
4
0.336
0.757
0.385
0.722
7
0.330
0.769
0.358
0.741
10
0.378
0.698
0.334
0.732
13
0.370
0.709
0.343
0.723
16
0.356
0.719
0.357
0.713
19
0.345
0.731
0.360
0.713
22
0.321
0.736
0.390
0.697
25
0.373
0.651
0.354
0.681
26
0.269
0.729
0.386
0.726
29
0.269
0.725
0.385
0.720
32
0.271
0.720
0.382
0.720
35
0.275
0.710
0.376
0.719
38
0.270
0.716
0.380
0.718
41
0.284
0.691
0.367
0.715
44
0.289
0.683
0.363
0.714
47
0.293
0.675
0.366
0.711
99
Cumulative distribution function of probability of exceeding three different acceleration thresholds, Hs=2.5m
𝑺𝒅,𝒂𝒄𝒄𝒆𝒍. = 𝒂 · 𝑰𝑴𝒃
𝐻𝑠 = 2.5𝑚 𝑈𝑚 [m/s]
𝒂
𝒃
𝜷
𝑹𝟐
4
4.640
0.706
0.490
0.542
7
4.489
0.722
0.473
0.555
10
4.597
0.702
0.486
0.544
13
4.571
0.706
0.492
0.544
16
4.560
0.702
0.496
0.535
19
4.522
0.704
0.498
0.532
22
4.464
0.689
0.501
0.512
25
4.578
0.675
0.491
0.511
26
3.623
0.693
0.455
0.568
29
3.618
0.690
0.457
0.562
32
3.594
0.689
0.461
0.559
35
3.610
0.684
0.457
0.559
38
3.604
0.682
0.459
0.559
41
3.600
0.680
0.459
0.550
44
3.599
0.677
0.458
0.555
47
3.578
0.677
0.462
0.551
100
𝑯𝒔 = 𝟒. 𝟓𝒎
Cumulative distribution of probability of failure in LS1, Hs=4.5m 𝑺𝒅,𝒚𝒊𝒆𝒍𝒅 = 𝒂 · 𝑰𝑴𝒃
𝐻𝑠 = 4.5𝑚 𝑈𝑚 [m/s]
𝒂
𝒃
𝜷
𝑹𝟐
4
118.079
0.716
0.394
0.681
7
114.863
0.732
0.370
0.701
10
130.976
0.664
0.343
0.692
13
128.514
0.673
0.356
0.682
16
124.214
0.681
0.369
0.671
19
118.953
0.700
0.386
0.670
22
114.072
0.691
0.399
0.650
25
125.454
0.638
0.369
0.640
26
96.838
0.682
0.397
0.680
29
97.184
0.678
0.394
0.677
32
97.629
0.673
0.392
0.675
35
98.607
0.667
0.387
0.676
38
98.911
0.661
0.386
0.672
41
100.375
0.653
0.383
0.668
44
101.993
0.646
0.379
0.667
47
103.299
0.639
0.378
0.660 101
Cumulative distribution of probability of failure in LS3, Hs=4.5m
𝑺𝒅,𝒃𝒖𝒄𝒌𝒍𝒆 = 𝒂 · 𝑰𝑴𝒃
𝐻𝑠 = 4.5𝑚 𝑈𝑚 [m/s]
𝒂
𝒃
𝜷
𝑹𝟐
4
0.336
0.758
0.385
0.721
7
0.329
0.770
0.357
0.741
10
0.381
0.693
0.331
0.733
13
0.369
0.710
0.343
0.723
16
0.357
0.718
0.357
0.712
19
0.340
0.739
0.375
0.712
22
0.323
0.732
0.386
0.696
25
0.368
0.659
0.355
0.683
26
0.269
0.728
0.388
0.724
29
0.271
0.721
0.383
0.721
32
0.272
0.718
0.381
0.720
35
0.275
0.710
0.376
0.720
38
0.277
0.703
0.374
0.717
41
0.282
0.694
0.370
0.714
44
0.288
0.684
0.365
0.714
47
0.292
0.675
0.365
0.708
102
Cumulative distribution function of probability of exceeding three different acceleration thresholds, Hs=4.5m
𝑺𝒅,𝒂𝒄𝒄𝒆𝒍. = 𝒂 · 𝑰𝑴𝒃
𝐻𝑠 = 4.5𝑚 𝑈𝑚 [m/s]
𝒂
𝒃
𝜷
𝑹𝟐
4
4.640
0.706
0.490
0.542
7
4.462
0.725
0.476
0.554
10
4.598
0.702
0.485
0.544
13
4.586
0.704
0.491
0.543
16
4.564
0.701
0.495
0.535
19
4.544
0.701
0.496
0.532
22
4.453
0.689
0.500
0.517
25
4.583
0.675
0.490
0.511
26
3.627
0.692
0.456
0.568
29
3.611
0.690
0.458
0.562
32
3.605
0.688
0.459
0.559
35
3.602
0.685
0.459
0.559
38
3.584
0.684
0.460
0.559
41
3.587
0.682
0.462
0.550
44
3.594
0.678
0.458
0.555
47
3.574
0.678
0.464
0.550
103
𝑯𝒔 = 𝟔. 𝟓𝒎
Cumulative distribution of probability of failure in LS1, Hs=6.5m 𝑺𝒅,𝒚𝒊𝒆𝒍𝒅 = 𝒂 · 𝑰𝑴𝒃
𝐻𝑠 = 6.5𝑚 𝑈𝑚 [m/s]
𝒂
𝒃
𝜷
𝑹𝟐
4
117.825
0.717
0.394
0.683
7
115.385
0.730
0.368
0.700
10
131.492
0.661
0.344
0.691
13
128.590
0.673
0.356
0.682
16
124.797
0.681
0.366
0.671
19
120.850
0.693
0.375
0.670
22
113.339
0.696
0.396
0.653
25
129.651
0.617
0.367
0.636
26
97.536
0.678
0.395
0.681
29
96.792
0.678
0.399
0.673
32
98.249
0.672
0.387
0.677
35
98.588
0.666
0.389
0.675
38
99.134
0.660
0.390
0.670
41
100.568
0.652
0.382
0.668
44
102.418
0.645
0.377
0.668
47
104.305
0.636
0.374
0.663 104
Cumulative distribution of probability of failure in LS3, Hs=6.5m
𝑺𝒅,𝒃𝒖𝒄𝒌𝒍𝒆 = 𝒂 · 𝑰𝑴𝒃
𝐻𝑠 = 6.5𝑚 𝑈𝑚 [m/s]
𝒂
𝒃
𝜷
𝑹𝟐
4
0.336
0.758
0.386
0.722
7
0.329
0.770
0.357
0.741
10
0.380
0.695
0.332
0.733
13
0.371
0.710
0.347
0.722
16
0.357
0.720
0.353
0.714
19
0.345
0.731
0.363
0.712
22
0.320
0.737
0.385
0.699
25
0.378
0.641
0.355
0.679
26
0.272
0.723
0.385
0.725
29
0.270
0.722
0.389
0.718
32
0.273
0.717
0.377
0.722
35
0.275
0.710
0.379
0.719
38
0.277
0.704
0.379
0.715
41
0.282
0.694
0.370
0.714
44
0.289
0.683
0.364
0.715
47
0.295
0.672
0.361
0.711
105
Cumulative distribution function of probability of exceeding three different acceleration thresholds, Hs=6.5m
𝑺𝒅,𝒂𝒄𝒄𝒆𝒍. = 𝒂 · 𝑰𝑴𝒃
𝐻𝑠 = 6.5𝑚 𝑈𝑚 [m/s]
𝒂
𝒃
𝜷
𝑹𝟐
4
4.640
0.706
0.490
0.542
7
4.481
0.723
0.472
0.555
10
4.593
0.702
0.487
0.544
13
4.586
0.704
0.492
0.543
16
4.562
0.702
0.495
0.534
19
4.509
0.705
0.499
0.533
22
4.452
0.690
0.497
0.519
25
4.578
0.674
0.492
0.512
26
3.635
0.691
0.454
0.569
29
3.609
0.690
0.461
0.560
32
3.614
0.687
0.458
0.560
35
3.596
0.686
0.460
0.557
38
3.575
0.684
0.463
0.558
41
3.594
0.681
0.462
0.547
44
3.586
0.678
0.460
0.554
47
3.592
0.676
0.460
0.551
106
DTU Civil Engineering Department of Civil Engineering Technical University of Denmark
Brovej, Building 118 2800 Kgs. Lyngby Telephone 45 25 17 00
107
