Comparative Study of Time-Domain Fatigue Assessments for an Offshore Wind Turbine Jacket Substructure by Using Conventional Grid-Based and Monte Carlo Sampling Methods

Currently, in the design standards for environmental sampling to assess long-term fatigue damage, the grid-based sampling method is used to scan a rectangular grid of meteorological inputs. However, the required simulation cost increases exponentially with the number of environmental parameters, and considerable time and effort are required to characterise the statistical uncertainty of offshore wind turbine (OWT) systems. In this study, a K-type jacket substructure of an OWT was modelled numerically. Time rather than frequency-domain analyses were conducted because of the high nonlinearity of the OWT system. The Monte Carlo (MC) sampling method is well known for its theoretical convergence, which is independent of dimensionality. Conventional grid-based and MC sampling methods were applied for sampling simulation conditions from the probability distributions of four environmental variables. Approximately 10,000 simulations were conducted to compare the computational efficiencies of the two sampling methods, and the statistical uncertainty of the distribution of fatigue damage was assessed. The uncertainty due to the stochastic processes of the wave and wind loads presented considerable influence on the hot-spot stress of welded tubular joints of the jacket-type substructure. This implies that more simulations for each representative short-term environmental condition are required to derive the characteristic fatigue damage. The characteristic fatigue damage results revealed that the MC sampling method yielded the same error level for Grids 1 and 2 (2443 iterations required for both) after 1437 and 516 iterations for Kand KK-joint cases, respectively. This result indicated that the MC method has the potential for a high convergence rate.


Introduction
Pollution due to energy production is harmful for our living environment; therefore, an alternative energy source is increasingly required.Offshore wind turbine (OWT) technology offers an opportunity in this field [1].According to the Global Wind Energy Council (GWEC), the global wind turbine installation capacity was 539.1 GW by the end of 2017 and is estimated to reach 840.9 GW in 2022 [2].Jacket-type substructures are commonly used for the OWT systems for water depths of 20 to 50 m.The tubular joint connections in jacket substructures are subject to approximately 10 9 load cycles during their designed lifetime of 20 years [3] and thus are the critical locations that may undergo fatigue failure.Time-domain rather than frequency-domain analyses for OWTs are typically used to handle Energies 2018, 11, 3112 2 of 17 the nonlinearities and interaction problems of wave and wind loads.The primary problem with the time-domain fatigue analyses is that the process is time-consuming because it employs the conventional grid-based approach.Thousands of combinations of environmental parameters, such as wind speed U w , wave height H s , wave period T p , and load direction θ [4], should be considered for long-term fatigue.The uncertainty of fatigue damage distribution in time analysis should be considered because of the stochastic processes of the wind and wave simulations [5].Moreover, the partial safety factor approach [6,7] is often used in wind turbine fatigue design, and the uncertainty can be minimised by performing more simulations [8].Dong [9] reduced the relative statistical uncertainties by conducting 20 simulations each for 400 short-term environmental conditions; that is, 8000 simulations in total.Both environmental parameters and uncertainties can cause 'high dimensionality'; thus, the number of required simulations would render the problem intractable.
Monte Carlo (MC) sampling method is widely used to address high-dimensionality problems in finance [10] and is applied in OWT analysis under probabilistic frameworks.Bilionis [11] studied the cross-sectional fatigue at the foundation of a monopile design by conducting a MC simulation of the wind (velocity) and wave (height and period) characteristics.Vahdatirad [12] quantified the uncertainties due to the stochastic dynamic stiffness of the foundation of large OWTs; the uncertainties were related to the first natural frequency of a turbine that was supported by surface footing.Soil elastic modulus and layer depth were considered as random variables with lognormal distributions; thus, the MC process was used to conduct a probabilistic simulation.Dong [13] evaluated the uncertainty due to local strain and fatigue crack initiation life of welded joints and used an approach based on MC simulation to evaluate the uncertainty due to the randomness of the geometrical and material parameters.Jensen [14] combined the MC sampling approach and the first-order reliability method to obtain a favourable estimation of the tail in the distribution of the wave-induced fatigue damage.Compared with the MC sampling approach, the combination method proposed by Jensen provided better estimation results.Subsequently, Horn [15] verified this method by including wind loads.Graf [16] compared the computational efficiencies of the grid-based and MC sampling methods by calculating millions of short-term damage equivalent loads of floating OWTs with the help of the supercomputer of the National Renewable Energy Laboratory.The MC method presents significantly reduced computational cost at high dimensions because the theoretical convergence rate is independent of dimensions.In general, the MC approach is more efficient when dimensions increase, compared with the grid-based method.
As an alternative to the time-domain methods, frequency-domain approaches can be highly efficient for fatigue analysis.Although inevitably limited to linear systems [17], the frequency-domain method is still useful for the initial design assessments [18].For comprehensive assessments, design calculations for OWTs are typically carried out in the time domain, suggested by the standards of International Electrotechnical Commission IEC 61400-3 [19] and the guidelines developed by Det Norske Veritas (Norway) and Germanischer Lloyd (Germany) Classification Societies DNVGL-ST-0126 [20].Because the environmental parameters and uncertainties can cause 'high dimensionality' for lifetime assessments of the OWTs, we aim to apply the MC approach to the time-domain analysis for reducing the computational cost.The purpose of this study was also to compare the number of simulations required to obtain the same level of fatigue damage at the welded tubular joints of the jacket-type substructure by employing the conventional grid-based and MC sampling methods.The uncertainty of the fatigue analysis was also discussed.This study conducted time-domain simulations to not only accurately reproduce the long-term evolution of the wind and wave loads, but also evaluate the corresponding structural response.Time-domain analyses were conducted for calculating the hot-spot stress and the stress range by using the rainflow-counting algorithm [21].The fatigue limit is commonly characterized as a stress-cycle (S-N) curve, which plots the magnitude of an alternating stress versus the number of cycles to failure for a given material and joint type.The S-N curve and Palmgren-Miner's rule were applied to assess the fatigue life.The entire procedure is depicted in Figure 1, and every module is sequentially described in the following sections.The sea load and wind load were generated via in-house code HydroCRest and the unsteady blade element momentum method (UBEM) respectively.The calculations were conducted on the basis of the guidelines of DNV GL [20,22], in compliance with IEC 61400-1 [23].
Energies 2018, 11, x 3 of 17 load and wind load were generated via in-house code HydroCRest and the unsteady blade element momentum method (UBEM) respectively.The calculations were conducted on the basis of the guidelines of DNV GL [20,22], in compliance with IEC 61400-1 [23].
Figure 1.Flowchart of the methodology; UBEM and S-N curve denote the unsteady blade element momentum method and the stress-cycle curve, respectively.

Long-Term Environmental Statistics and Sampling
The long-term environmental statistics were based on a preliminary site survey that was conducted from December 1991 to August 1999 at the THL3 weather station near the Fuhai Offshore Wind Farm in Taiwan (Figure 2).The scatter diagrams and joint probability tables were collated to produce a covariance matrix of the five considered environmental load parameters, namely wave height, period, wave direction, wind speed, and wind direction.In this study, we assumed constant wind and wave direction.The mean wind speed   was in the range of 0 to 30 m/s with an increment of 5 m/s, the wind direction  was in the range of 0° to 360° with an increment of 30°, the significant wave height   was in the range of 0.25 to 4.75 m with an increment of 0.5 m, and the spectral peak period   was in the range of 3.5 to 14.5 s with an increment of 1 s.By using the four parameters, 8640 combinations were formed.Note that some conditions do not exist (i.e., zero probability) based on the joint probability table of the wave scatter   and   for a given   and the joint probability table of the wind scatter   and .Therefore, 2443 weather conditions can be derived from the joint probability density distribution of the characteristic parameters   , ,   , and   as presented in Equation ( 1 This equation was used to generate stochastic weather states for 20 years.The states served as the basis for obtaining statistically accurate wind and wave load calculations to conduct fatigue lifetime evaluations.

Long-Term Environmental Statistics and Sampling
The long-term environmental statistics were based on a preliminary site survey that was conducted from December 1991 to August 1999 at the THL3 weather station near the Fuhai Offshore Wind Farm in Taiwan (Figure 2).The scatter diagrams and joint probability tables were collated to produce a covariance matrix of the five considered environmental load parameters, namely wave height, period, wave direction, wind speed, and wind direction.In this study, we assumed constant wind and wave direction.The mean wind speed U w was in the range of 0 to 30 m/s with an increment of 5 m/s, the wind direction θ was in the range of 0 • to 360 • with an increment of 30 • , the significant wave height H s was in the range of 0.25 to 4.75 m with an increment of 0.5 m, and the spectral peak period T p was in the range of 3.5 to 14.5 s with an increment of 1 s.By using the four parameters, 8640 combinations were formed.Note that some conditions do not exist (i.e., zero probability) based on the joint probability table of the wave scatter H s and T p for a given U w and the joint probability table of the wind scatter U w and θ.Therefore, 2443 weather conditions can be derived from the joint probability density distribution of the characteristic parameters U w , θ, H s , and T p as presented in Equation (1): Energies 2018, 11, x 3 of 17 load and wind load were generated via in-house code HydroCRest and the unsteady blade element momentum method (UBEM) respectively.The calculations were conducted on the basis of the guidelines of DNV GL [20,22], in compliance with IEC 61400-1 [23].
Figure 1.Flowchart of the methodology; UBEM and S-N curve denote the unsteady blade element momentum method and the stress-cycle curve, respectively.

Long-Term Environmental Statistics and Sampling
The long-term environmental statistics were based on a preliminary site survey that was conducted from December 1991 to August 1999 at the THL3 weather station near the Fuhai Offshore Wind Farm in Taiwan (Figure 2).The scatter diagrams and joint probability tables were collated to produce a covariance matrix of the five considered environmental load parameters, namely wave height, period, wave direction, wind speed, and wind direction.In this study, we assumed constant wind and wave direction.The mean wind speed   was in the range of 0 to 30 m/s with an increment of 5 m/s, the wind direction  was in the range of 0° to 360° with an increment of 30°, the significant wave height   was in the range of 0.25 to 4.75 m with an increment of 0.5 m, and the spectral peak period   was in the range of 3.5 to 14.5 s with an increment of 1 s.By using the four parameters, 8640 combinations were formed.Note that some conditions do not exist (i.e., zero probability) based on the joint probability table of the wave scatter   and   for a given   and the joint probability table of the wind scatter   and .Therefore, 2443 weather conditions can be derived from the joint probability density distribution of the characteristic parameters   , ,   , and   as presented in Equation ( 1 This equation was used to generate stochastic weather states for 20 years.The states served as the basis for obtaining statistically accurate wind and wave load calculations to conduct fatigue lifetime evaluations.This equation was used to generate stochastic weather states for 20 years.The states served as the basis for obtaining statistically accurate wind and wave load calculations to conduct fatigue lifetime evaluations.

Grid-Based Method
A 'grid-based method' or 'lumping of load cases' is widely used in the offshore industry.This method is the most straightforward approach for estimating fatigue damage by selecting increments for each environmental variable and assessing all combinations of the inputs.This approach enables rapid calculation of all load cases accumulated from the long-term environmental statistics because only one simulation is required for each grid.However, for high-dimensional variables, the grid-based method becomes intractable when the number of required simulations grows exponentially with the dimension [16].In this study, four variables were considered-wind speed, wind direction, wave height, and wave period.Thus, the study was a 'four-dimensional fatigue analysis'.

MC Method
The MC method is extensively used in high-dimensionality cases to calculate the expected value α of a certain problem.The convergence rate of the MC method is based on the central limit theorem and the law of large numbers as follows: The sample mean α n has approximately the same distribution as the expected value α added to the standard Gaussian distribution N(0, 1) multiplied by the ratio of the standard deviation σ and the square root of the number of samples n.The approximation implies that the convergence rate of the MC method is of the order of 1 √ n and is independent of the dimension.Therefore, the method provides a high-speed estimation even at high dimensions.The error ie n between the sample mean and the true expected value can be written as follows: Equation (3) specifies that the error decreases when n increases.The probability of the error based on the selected intervals is given as follows: Equation ( 4) provides an approximate confidence interval by replacing σ with σ n , thus implying that a reasonable estimate can be obtained even for small computations.In Equation ( 4), c 1 and c 2 are usually selected as ±1.96 to derive the 95% confidence interval.

Substructure Model Description and Load Generation
This study was focused on a fixed K-type jacket substructure of an OWT designed for the Fuhai Offshore Wind Farm for a water depth of 23 m (Figure 3).The wind turbine has a power of 3.6 MW in the preliminary design stage, the diameter of the rotor is 120 m, the cut-in wind speed is 3 m/s, the cut-out wind speed is 25 m/s, the nominal wind speed is 12 m/s, the tower height is 68.39 m, and the jacket height is 39.22 m.The entire model included the tower with the simplified mass of the rotor nacelle assembly applied on the top of the tower and the jacket substructure that was rooted on four piles for which the interaction with the soil was reproduced using nonlinear elastic spring connectors.The aerodynamic forces acting on the wind turbine were calculated, and then, the equivalent load was applied to the top of the wind tower as a concentrated force in the structural finite element model.Only x and y components were considered, and the z component was neglected because the tilt angle is only 5 • .The bending moment at the root of a rotor blade was assumed to be in balance because of the three blades.

Sea Loads
For each sea state, the selected Joint North Sea Wave Project (JONSWAP) sea spectrum [24] and directional spreading function were applied to obtain an irregular, time-varying flow field.The consequent wave loads for the flow field were determined using the Morison equation, as expressed in Equation ( 5) [20], and executed via the in-house HydroCRest code (sea-state generator/hydrodynamic solver).
where   and   are the drag and inertia coefficients, respectively; A is the area of the cylinder; D is the diameter of the cylinder;  is the density of seawater;  is the horizontal wave-induced velocity of seawater; and ̇ is the horizontal wave-induced acceleration of seawater.The coefficients were determined by following the guidelines of DNV GL [20], in which the slamming term   is neglected in normal wave conditions.On the basis of the superposition solution of the potential flow theory, an irregular sea surface was decomposed into an infinite number of regular component waves, which were formulated in terms of amplitude, direction, frequency, and phase.In the directional JONSWAP power spectrum presented in Equation ( 6) [25],   () is the wave spectrum defined by two parameters: the significant wave height (Hs) and the peak period (Tp), D(, ) is the directional spreading function,  is the wave frequency, and  is the direction of the wave.The wave amplitude is calculated using Equation ( 7) [25].The phase of each component wave was set randomly.For each load case the simulation was conducted for 10 minutes with time steps of 0.25 seconds, comprising totally 2400 steps.A single instruction, multiple data (SIMD) parallel computing algorithm was implemented in HydroCRest upon graphics processing unit (GPU) hardware architecture to accelerate the processing of thousands of component waves and element nodes for wave load simulations.Figure 4a presents the nodal forces and the total force as computed in HydroCRest. Figure 4b illustrates the directional distribution of the total force and moment.In Figure 4b, the directions of the total force and

Sea Loads
For each sea state, the selected Joint North Sea Wave Project (JONSWAP) sea spectrum [24] and directional spreading function were applied to obtain an irregular, time-varying flow field.The consequent wave loads for the flow field were determined using the Morison equation, as expressed in Equation ( 5) [20], and executed via the in-house HydroCRest code (sea-state generator/hydrodynamic solver).
where C D and C M are the drag and inertia coefficients, respectively; A is the area of the cylinder; D is the diameter of the cylinder; ρ is the density of seawater; ν is the horizontal wave-induced velocity of seawater; and .
ν is the horizontal wave-induced acceleration of seawater.The coefficients were determined by following the guidelines of DNV GL [20], in which the slamming term F s is neglected in normal wave conditions.
On the basis of the superposition solution of the potential flow theory, an irregular sea surface was decomposed into an infinite number of regular component waves, which were formulated in terms of amplitude, direction, frequency, and phase.In the directional JONSWAP power spectrum presented in Equation ( 6) [25], S J (ω) is the wave spectrum defined by two parameters: the significant wave height (H s ) and the peak period (T p ), D(ω, θ) is the directional spreading function, ω is the wave frequency, and θ is the direction of the wave.The wave amplitude is calculated using Equation ( 7) [25].The phase of each component wave was set randomly.For each load case the simulation was conducted for 10 minutes with time steps of 0.25 seconds, comprising totally 2400 steps.A single instruction, multiple data (SIMD) parallel computing algorithm was implemented in HydroCRest upon graphics processing unit (GPU) hardware architecture to accelerate the processing of thousands of component waves and element nodes for wave load simulations.Figure 4a presents the nodal forces and the total force as computed in HydroCRest. Figure 4b illustrates the directional distribution of the total force and moment.In Figure 4b, the directions of the total force and overturning moment (OTM) at each step were counted, whereas the major wave direction was set to the north direction.The probability distributions over the entire 10 minutes were plot, such that the times of the forces and OTMs in a certain direction were divided by the total number of time steps.Because the theoretical wave dynamics consist of oscillatory motions, the major directions of the total force fall in both the north and south directions, and the total OTM is in the east and west directions.In the simulations the major wave directions were assumed to be in line with the wind directions.
Energies 2018, 11, x 6 of 17 the north direction.The probability distributions over the entire 10 minutes were plot, such that the times of the forces and OTMs in a certain direction were divided by the total number of time steps.
Because the theoretical wave dynamics consist of oscillatory motions, the major directions of the total force fall in both the north and south directions, and the total OTM is in the east and west directions.
In the simulations the major wave directions were assumed to be in line with the wind directions.

Wind Loads
The short-term wind states were modelled on the basis of the normal turbulence model (NTM) in IEC 61400-1 [23] with a reference turbulence intensity Iref of 0.16.This modelling was conducted because the Taiwanese Ministry of Economic Affairs required a pilot wind turbine to be installed in the Fuhai Offshore Wind Farm with IEC 61400-1 Class IA compliance.By following the NTM, the standard deviation σU for a given mean wind speed at a hub height of 10 m (U10) can be calculated as follows [23]: To generate a realistic short-term wind state, the wind speed fluctuations were calibrated against practical onsite data by applying a nonlinear least squares regression analysis to the power spectral densities (PSD) of the wind speed data over consecutive 10-min periods.Then, the signals were reconstructed (Figure 5a).The regression curve amplitudes were scaled to produce a final signal that conforms to the NTM in IEC 61400-1.The wind fluctuations were further considered (Figure 5b) by randomising the phase shifts of the component harmonics.

Wind Loads
The short-term wind states were modelled on the basis of the normal turbulence model (NTM) in IEC 61400-1 [23] with a reference turbulence intensity I ref of 0.16.This modelling was conducted because the Taiwanese Ministry of Economic Affairs required a pilot wind turbine to be installed in the Fuhai Offshore Wind Farm with IEC 61400-1 Class IA compliance.By following the NTM, the standard deviation σ U for a given mean wind speed at a hub height of 10 m (U 10 ) can be calculated as follows [23]: To generate a realistic short-term wind state, the wind speed fluctuations were calibrated against practical onsite data by applying a nonlinear least squares regression analysis to the power spectral densities (PSD) of the wind speed data over consecutive 10-min periods.Then, the signals were reconstructed (Figure 5a).The regression curve amplitudes were scaled to produce a final signal that conforms to the NTM in IEC 61400-1.The wind fluctuations were further considered (Figure 5b) by randomising the phase shifts of the component harmonics.
are multiplied by the number of blades, as presented in Equations ( 9) and (10), to determine the total thrust and rotor torque: where   is the number of blades;  is the density of air;   is the rotor velocity;   and   are the lift and drag coefficients, respectively; and r is the length of the rotor blades.These equations can be easily found in studies pertaining to the blade element momentum theory.The effects of the wind turbine tower on the upstream flow field were modeled by assuming that the flow field resembled a potential flow around a circular cylinder (Figure 6a), such that the radial and angular components of the flow velocity at a considered point are given as follows [27]: where  ∞ is the undisturbed wind velocity and R is the radius of the tower.A wind shear profile was also included such that the wind velocity at height z for a specified hub height H can be expressed as follows: where the power law exponent α for offshore locations is 0.14, in accordance with DNVGL-ST-0126 [20].Because a large number of transient wind load calculations are necessary over the simulated life time, an unsteady blade element momentum method (UBEM) [26] was adopted to calculate the aerodynamic loads on the wind turbine.The UBEM discretises the rotor into a number of two-dimensional (2D) airfoil sections such that the axial and tangential loads on each 2D section may be calculated from the respective airfoil's lift and drag characteristics for the respective local relative flow velocity and angle.These local loads were integrated along the length of the rotor blades and are multiplied by the number of blades, as presented in Equations ( 9) and (10), to determine the total thrust and rotor torque: where n B is the number of blades; ρ is the density of air; U rel is the rotor velocity; C l and C d are the lift and drag coefficients, respectively; and r is the length of the rotor blades.These equations can be easily found in studies pertaining to the blade element momentum theory.The effects of the wind turbine tower on the upstream flow field were modeled by assuming that the flow field resembled a potential flow around a circular cylinder (Figure 6a), such that the radial and angular components of the flow velocity at a considered point are given as follows [27]: Energies 2018, 11, 3112 where U ∞ is the undisturbed wind velocity and R is the radius of the tower.
Energies 2018, 11, x 8 of 17 The UBEM model was validated against the power curve provided by the manufacturer of the 3.6-MW wind turbine [28] for the full range of normal operating conditions by considering the cutin, cut-out, and supra-nominal (pitch control) wind velocities.The model was found to appropriately correlate with the official data (Figure 6b).

Fatigue Lifetime Assessment
Beam elements were selected for the finite element model (Figure 3) analysed in Abaqus 2017.There are several approaches for fatigue analysis, such as the nominal stress approach, the hot-spot stress approach, the notch stress approach, and the crack propagation approach.The weld joints are the critical locations based on the experience of the offshore industry.We therefore implement the hot-spot stress method, recommended by the DNV GL standard [22], for evaluating the fatigue damage at the tubular joints.The stress raising effect can be considered by the hot-spot stress range by using the stress concentration factor (SCF) to modify the nominal stress:

Hot-Spot Stress Evaluation
In DNV GL, 'the hot-spot stress should be evaluated at eight spots around the circumference of the intersection', and the hot-spot stress range is derived by the superposition of the SCF as follows [22]:  A wind shear profile was also included such that the wind velocity at height z for a specified hub height H can be expressed as follows: where the power law exponent α for offshore locations is 0.14, in accordance with DNVGL-ST-0126 [20].
The UBEM model was validated against the power curve provided by the manufacturer of the 3.6-MW wind turbine [28] for the full range of normal operating conditions by considering the cut-in, cut-out, and supra-nominal (pitch control) wind velocities.The model was found to appropriately correlate with the official data (Figure 6b).

Fatigue Lifetime Assessment
Beam elements were selected for the finite element model (Figure 3) analysed in Abaqus 2017.There are several approaches for fatigue analysis, such as the nominal stress approach, the hot-spot stress approach, the notch stress approach, and the crack propagation approach.The weld joints are the critical locations based on the experience of the offshore industry.We therefore implement the hot-spot stress method, recommended by the DNV GL standard [22], for evaluating the fatigue damage at the tubular joints.The stress raising effect can be considered by the hot-spot stress range by using the stress concentration factor (SCF) to modify the nominal stress:

Hot-Spot Stress Evaluation
In DNV GL, 'the hot-spot stress should be evaluated at eight spots around the circumference of the intersection', and the hot-spot stress range is derived by the superposition of the SCF as follows [22]: ) Energies 2018, 11, 3112 9 of 17 In Equations ( 15)-( 22), SCF AC and SCF AS are the stress concentration factors at the crown for the axial load and at the saddle, respectively; SCF MIP and SCF MOP are the stress concentration factors for in-plane and out-of-plane moments, respectively; σ x is the maximum nominal stress due to the axial load; and σ my and σ mz are the maximum nominal stresses due to in-plane and out-of-plane bending moments, respectively.
The SCFs were calculated using a comprehensive set of simple joint parametric equations proposed by Efthymiou [29].The parametric equations depend on the joint configurations (T/Y, X, K, and KT), the applied load cases (axial load and in-plane and out-of-plane moments), the estimated location of the joint (chord side or brace side, and saddle or crown), and the geometry of the joint.To avoid long numerical and experimental procedures, the specified method was used in this study.Moreover, one hot spot at the brace side of the K joint and one hot spot at the brace side of the KK joint (Figure 3) were investigated.

S-N Curve and Welding Effect
The basic design of the stress-cycle (S-N) curve is presented in Equation ( 23) [22], where N is the number of cycles to failure in the stress range; ∆σ is the stress range derived at eight spots from Equations ( 15)-( 22); and a and m are constants that depend on the material and the environmental conditions, respectively.log N = log a − m log ∆σ (23) where m = 3.0 and a = 12.18 if N ≤ 1.8 × 10 6 , and m = 5.0 and a = 16.13 if N > 1.8 × 10 6 .Three types of S-N curves for the tubular joints (known as T-curves) in different environments were presented in [22] for air, seawater with cathodic protection, and seawater in which free corrosion is enabled.Both selected joints in this study were in the splash zone of the substructure; thus, the T-curve for tubular joints in seawater with cathodic protection was selected.All the tubular joints were assumed to be full-penetration welded joints at the design stage; thus, the welding effect was also considered.

Fatigue Damage Calculation
Fatigue damage calculation was based on the widely used Palmgren-Miner linear damage theory [30], such that the cumulative damage can be assumed to be the sum of the particular damage caused by each cycle and is expressed as follows: In Equation ( 24), D g is the accumulated fatigue damage, n i is the number of cycles for the stress range i, and N i is the number of cycles to failure for the stress range i, which can be evaluated using Equation (23).To deal with the long-term stress, the rainflow-counting algorithm was used for the fatigue damage evaluation.The stress range, mean stress and cycles of the irregularity, and the complex long-term stress signal were determined using the algorithm.By using an open-source rainflow-counting algorithm MATLAB code [21,31], the expected lifetime T of the tubular joints can be calculated using Equation (25), where T 0 is the time interval of the stress history.

Results and Discussion
The accumulated fatigue damage of the tubular joints of the jacket substructure of an OWT was mainly due to the wind loads.The contributions of the wave loads, however, were still considerable [3].Moreover, the wind and wave load simulations were stochastic and irregular processes because they are natural phenomena.Both phase angles of the wave and the wind turbulence models were randomly generated and caused variation in the hot-spot stress distribution.The random nature has a considerable influence on the fatigue damage, thus leading to underestimating or overestimating the fatigue lifetime.In the first part of this section, the uncertainty of fatigue damage is discussed for the cases of only wind load and only wave load.In the second part, 5000 simulations were conducted, with each including 10-min hydrodynamic, aerodynamic, and structural analysis based on the MC method, and 4886 simulations were conducted by using the grid-based method (approximately 68.6 simulation days in total).The computational effort of the hydrodynamic and aerodynamic simulations can be summed up to 1.2 days with the parallelisation technology.The structural simulations in Abaqus were completed in 63.9 days by using one octa-core Intel Core i7-4790 CPU operated at 3.60 GHz.

Uncertainty of Fatigue Damage
The influence of the stochastic process was discussed by splitting the load conditions into four wind conditions and three wave conditions.For the wind conditions, 60 short-term (10 min) simulations were conducted.A hot spot of the KK joint presented in Figure 3 was selected for demonstration.The fatigue damage values for the case in which only the wind load is considered are presented in Figure 7a, whereas the mean wind speed was set to 7.5 m/s and the wind direction was 0 • .Fatigue damage varied with the load case for the 60 simulations under identical load conditions; the maximum damage was 5.007 × 10 −12 and the minimum was 1.020 × 10 −12 .The difference between the minimum and maximum values is considerable because of the stochastic process of wind simulations, thus causing simulation deviation for each representative short-term environmental condition if the cycle numbers of simulations were limited.The variation in the characteristic fatigue damage (CFD) for different simulation samples can be calculated using Equation (26), as presented in Figure 7a.Approximately seven simulations must be conducted to minimise the CFD error for reducing the uncertainty effect.A highly representative result can be achieved for each load condition by increasing the samples.
Three other wind conditions were also evaluated with a wind direction of 0 • and mean wind speed set to 12.5, 17.5, and 22.5 m/s.The fatigue damage distributions were obtained using the 60 simulations for each load condition (Figure 7b).The fatigue damage ranges due to the hot-spot stress were from 1.812 × 10 −9 to 7.754 × 10 −9 for the wind speed of 12.5 m/s, from 4.376 × 10 −9 to 1.566 × 10 −8 for the wind speed of 17.5 m/s, and from 1.418 × 10 −8 to 3.663 × 10 −8 for the wind speed of 22.5 m/s.Because the fatigue lifetime of the welded tubular joints of the jacket substructure was shown to be dominated by wind loads [3], the variation in fatigue damage was mainly influenced by wind loads.
The maximum, minimum, average, coefficient of variation (CV), and range between the extreme values of the aforementioned four wind conditions are listed in Table 1.In the worst case, the fatigue damage was overestimated by a value in the range of 154.0 to 224.6% under the four wind conditions if only one simulation was conducted because of the extreme value in the 60 simulations.For the fatigue lifetime in Equation ( 25), the lifetime of a hot spot was underestimated by a value in the range of 44.5 to 64.9%.The fatigue damage was underestimated by a value in the range of 42.1 to 59.6%, thus resulting in a simulation deviation.Moreover, the lifetime was overestimated by a value in the range 167.7 to 237.5% compared with the average lifetime.The deviation of the extreme values was significant; thus, the uncertainty due to the random effect of the load simulation was considerable.
simulations, thus causing simulation deviation for each representative short-term environmental condition if the cycle numbers of simulations were limited.The variation in the characteristic fatigue damage (CFD) for different simulation samples can be calculated using Equation (26), as presented in Figure 7a.Approximately seven simulations must be conducted to minimise the CFD error for reducing the uncertainty effect.A highly representative result can be achieved for each load condition by increasing the samples.The CV, also known as the relative standard deviation, is defined as the ratio of the standard deviation (STD) to the average (AVG) of data, as presented in Equation (27), and is widely used to measure the dispersion of data.The CV of fatigue damage was evaluated using the 60 simulations for the case in which only the wind load was considered and was up to 0.420 for the wind speed of 7.5 m/s.
The simulation deviation or error (%) for the wind-load-alone cases can be calculated using Equation (28).The CFDs of the 60 samples (i.e., the average fatigue damage of 600-min or 10-h simulations) were adopted as the data and termed as the 'long-term average'.The simulation deviation of the wind is presented in Figure 8.
We observed that the error curve converges when the number of simulations increases, thus implying that more simulations should be conducted to obtain a smaller simulation error for each load condition.The error curves based on the 60 short-term simulations lead to the prediction that the simulation errors for CFD can be reduced to less than 5% after 19, 13, 31, and 3 simulations for wind speeds of 7.5, 12.5, 17.5, and 22.5 m/s, respectively.Although the CV of the 17.5 m/s wind load condition was smaller than those of the 7.5 m/s and 12.5 m/s wind load conditions, the former case still requires the highest amount of cases to converge within a 5% error.Next, the variation in the fatigue damage for the wave-load-alone cases was examined.The hotspot stress distribution for welded tubular joints of the jacket substructures indicated that the wave load is a crucial factor for fatigue lifetime estimation in the OWT industry.To assess the uncertainty of wave-induced fatigue damage in welded tubular joints, 60 short-term (10 min) simulations for three wave environmental conditions were performed.The wave heights were set to 0.75, 1.25, and 1.75 m, and the wave direction and period were set to 0° and 6.5 s, respectively, for all three load conditions.
The fatigue damage, presented in Table 2, was overestimated by a value in the range of 225.2 to 325.0% or underestimated by a value in the range of 32.3 to 38.9% because of the maximum and minimum damage in the 60 simulations.The deviations in the extreme values for the wave-loadalone cases are as high as those in the wind-load-alone cases.Next, the variation in the fatigue damage for the wave-load-alone cases was examined.The hot-spot stress distribution for welded tubular joints of the jacket substructures indicated that the wave load is a crucial factor for fatigue lifetime estimation in the OWT industry.To assess the uncertainty of wave-induced fatigue damage in welded tubular joints, 60 short-term (10 min) simulations for three wave environmental conditions were performed.The wave heights were set to 0.75, 1.25, and 1.75 m, and the wave direction and period were set to 0 • and 6.5 s, respectively, for all three load conditions.
The fatigue damage, presented in Table 2, was overestimated by a value in the range of 225.2 to 325.0% or underestimated by a value in the range of 32.3 to 38.9% because of the maximum and minimum damage in the 60 simulations.The deviations in the extreme values for the wave-load-alone cases are as high as those in the wind-load-alone cases.For a Gaussian sea surface, the following simplification of ocean surface wave elevation is used extensively in offshore engineering [9]: x ω j ∆ω j cos ω j t − ε j (29) where S + x is the directional wave spectrum.The phase angle ε j of waves was assumed to be uniformly distributed over (0, 2π) and was randomly generated in the time-domain simulation because of the stochastic and irregular manner of the wave simulation based on the natural phenomena.This random generation affected the distribution of the hot-spot stress range and yielded uncertainties in the wave-induced fatigue damage.The simulation error of the wave condition was obtained using Equation ( 28) and is illustrated in Figure 9.The error curves, such as those of the wind-load-alone cases, converged as the number of simulations increased.Inevitably, when considering the combined wind-wave conditions, the uncertainty of the fatigue damage caused by the stochastic process of the simulations increased further.Therefore, more simulations are required to reduce the simulation deviation.Thus, the number of required evaluations increases and the problem becomes intractable.
Energies 2018, 11, x 13 of 17 obtained using Equation ( 28) and is illustrated in Figure 9.The error curves, such as those of the wind-load-alone cases, converged as the number of simulations increased.Inevitably, when considering the combined wind-wave conditions, the uncertainty of the fatigue damage caused by the stochastic process of the simulations increased further.Therefore, more simulations are required to reduce the simulation deviation.Thus, the number of required evaluations increases and the problem becomes intractable.

High-Dimensional Fatigue Analysis
Four variables were included in the four-dimensional fatigue analysis-wind speed, wind direction, wave height, and wave period.The stochastic processes of both wind and wave simulations could influence the CFD value, and the random effects of phase sets were also implicitly included in the 'dimension'.For the grid-based method, 2443 combinations (10-min short-term simulations) were conducted to evaluate the cumulative fatigue damage by summing the product of the fatigue damage derived from each environmental condition and the occurrence probability.Two sets of simulations for each condition were calculated (4886 simulations for the grid-based method), and the CFD results for each set were termed Grid 1 and Grid 2. The average of the two results, termed as the grid average, was considered as a representative for the grid-based method.We observed that the Grid 1 and Grid 2 results differed even under identical environmental conditions.Next, the MC approach was implemented by randomly sampling the environmental variables on the basis of the probability distribution of each variable and by calculating the 10-min short-term fatigue damage for each combination that was sampled.The CFD values of n samples were thus obtained.In total, 5000 samples were selected for the preliminary work.
Figure 10 presents the CFD results of the grid-based and MC methods with 95% confidence intervals for the hot spots around the circumference of the connections of the selected K and KK joints.The x-axis and y-axis represent the number of samples and the CFD result, respectively.The CFD values of the MC method are presented as solid lines.In Figure 5, the blue solid line is the actual CFD result derived using the 5000 samples, and the green and black lines are the upper and lower confidence bounds, respectively.The CFD results of the grid-based method are indicated by symbols (Grid 1 and Grid 2).The average of the grid-based method was termed the 'grid average'.Although more simulations are required to derive the long-term CFD of the grid-based method, the grid-based

High-Dimensional Fatigue Analysis
Four variables were included in the four-dimensional fatigue analysis-wind speed, wind direction, wave height, and wave period.The stochastic processes of both wind and wave simulations could influence the CFD value, and the random effects of phase sets were also implicitly included in the 'dimension'.For the grid-based method, 2443 combinations (10-min short-term simulations) were conducted to evaluate the cumulative fatigue damage by summing the product of the fatigue damage derived from each environmental condition and the occurrence probability.Two sets of simulations for each condition were calculated (4886 simulations for the grid-based method), and the CFD results for each set were termed Grid 1 and Grid 2. The average of the two results, termed as the grid average, was considered as a representative for the grid-based method.We observed that the Grid 1 and Grid 2 results differed even under identical environmental conditions.Next, the MC approach was implemented by randomly sampling the environmental variables on the basis of the probability distribution of each variable and by calculating the 10-min short-term fatigue damage for each combination that was sampled.The CFD values of n samples were thus obtained.In total, 5000 samples were selected for the preliminary work.
Figure 10 presents the CFD results of the grid-based and MC methods with 95% confidence intervals for the hot spots around the circumference of the connections of the selected K and KK joints.The x-axis and y-axis represent the number of samples and the CFD result, respectively.The CFD values of the MC method are presented as solid lines.In Figure 5, the blue solid line is the actual CFD result derived using the 5000 samples, and the green and black lines are the upper and lower confidence bounds, respectively.The CFD results of the grid-based method are indicated by symbols (Grid 1 and Grid 2).The average of the grid-based method was termed the 'grid average'.Although more simulations are required to derive the long-term CFD of the grid-based method, the grid-based method was used twice and the grid average was selected as the data for comparison because of the high computational cost (2443 simulations are required for the wind, wave, and structural analysis for one grid-based result).CFD values of the K joint (Figure 10a) derived from Grid 1 and Grid 2 were 4.108 × 10 and 4.395 × 10 , respectively, and the grid average was 4.251 × 10 .For the MC result, the value for the 5000 load cases was 4.240 × 10 .The fatigue lifetimes were 463, 432, 447, and 448 years for the Grid 1, Grid 2, the grid average, and the MC results, respectively.The variation between Grid 1 and Grid 2 indicated the uncertainty regarding the fatigue damage in four-dimensional fatigue analysis.Both the random effects of wind and wave simulations influenced the results of the case in which both wind and wave loads were considered.
One hot spot of the KK joint was also evaluated for validation (Figure 10b).The CFD ratio of Grid 2 to Grid 1 was 106.5%, thus revealing the random effect of the simulations.The CFD of the MC method (blue solid line) roughly converged after 500 samples and was near to the grid average.
The errors of Grid 1, Grid 2, and the MC values of both K and KK joints were estimated (Figure 11) by selecting the grid average for comparison.In Figure 11, the red solid line presents the grid average and is marked as the data line (0% error).For the K-joint case, the error estimations of Grid 1 and Grid 2 were approximately 3.38% .The grid-based method required 2443 simulations, whereas the CFD of the MC method converged after 1437 simulations, thus leading to an improvement in the computational efficiency.Moreover, the CFD of the MC method was within only 5% error after simulating 617 samples when compared with the grid average.A similar conclusion was drawn for the KK-joint case.CFD values of the K joint (Figure 10a) derived from Grid 1 and Grid 2 were 4.108 × 10 −8 and 4.395 × 10 −8 , respectively, and the grid average was 4.251 × 10 −8 .For the MC result, the value for the 5000 load cases was 4.240 × 10 −8 .The fatigue lifetimes were 463, 432, 447, and 448 years for the Grid 1, Grid 2, the grid average, and the MC results, respectively.The variation between Grid 1 and Grid 2 indicated the uncertainty regarding the fatigue damage in the four-dimensional fatigue analysis.Both the random effects of wind and wave simulations influenced the results of the case in which both wind and wave loads were considered.
One hot spot of the KK joint was also evaluated for validation (Figure 10b).The CFD ratio of Grid 2 to Grid 1 was 106.5%, thus revealing the random effect of the simulations.The CFD of the MC method (blue solid line) roughly converged after 500 samples and was near to the grid average.
The errors of Grid 1, Grid 2, and the MC values of both K and KK joints were estimated (Figure 11) by selecting the grid average for comparison.In Figure 11, the red solid line presents the grid average and is marked as the data line (0% error).For the K-joint case, the error estimations of Grid 1 and Grid 2 were approximately ±3.38%.The grid-based method required 2443 simulations, whereas the CFD of the MC method converged after 1437 simulations, thus leading to an improvement in the computational efficiency.Moreover, the CFD of the MC method was within only 5% error after simulating 617 samples when compared with the grid average.A similar conclusion was drawn for the KK-joint case.

Concluding Remarks
In this study, two sampling methods were presented for the fatigue estimation of the jacket substructure of an OWT-a conventional grid-based method and the MC method-to derive the fatigue damage of a K-type jacket substructure designed for use in the Fuhai Offshore Wind Farm in Taiwan.The long-term environmental statistics were based on a preliminary site survey conducted near the weather station to reproduce the environmental conditions.The uncertainty of the fatigue damage caused by the stochastic process of aerodynamic and hydrodynamic simulations were discussed, and the hot spots of both K and KK joints were selected for studying the uncertainties and the computational efficiencies of the cases in which only the wind and wave loads were considered.For the cases in which only wind loads were considered, the simulation error of the fatigue damage was in the range of 154.0 to 224.6% (maximum) or 42.1 to 59.6% (minimum) when compared with the long-term average.For the cases in which only wave loads were considered, the error was in the range of 225.2 to 325.0% (maximum) or 32.3 to 38.9% (minimum).The error was considerable when the number of simulations were limited, thus indicating that more simulations are required for obtaining the representative CFD.The simulation results of the grid-based method also indicated uncertainties for the combination of wind and wave loads.The ratio of Grid 2 to Grid 1 was 106.9% for the K joint and 106.5% for the KK joint.Moreover, the MC method was implemented by randomly sampling the environmental variables based on the probability distribution of each variable.The number of simulations can be potentially reduced using the MC method.The CFD results revealed that the MC approach achieved the same error level for both Grid 1 and Grid 2 (2443 iterations) after 1437 and 516 iterations for the K-and KK-joint cases, respectively.This indicates that the MC method has a high convergence rate.

Figure 1 .
Figure 1.Flowchart of the methodology; UBEM and S-N curve denote the unsteady blade element momentum method and the stress-cycle curve, respectively.

Figure 2 .
Figure 2. Scatter diagrams of wave height (H s ) and peak period (Tp) on the basis of a preliminary site survey that was conducted from December 1991 to August 1999 at the THL3 weather station, Taiwan.

Energies 2018, 11 , x 5 of 17 Figure 3 .
Figure 3. Offshore wind turbine (OWT) model with the K-type jacket substructure.(a) Wind tower, jacket substructure with the K joint (left) and KK joint (right), and four piles; (b) Top view of the substructure; (c) Front view of the substructure.

Figure 3 .
Figure 3. Offshore wind turbine (OWT) model with the K-type jacket substructure.(a) Wind tower, jacket substructure with the K joint (left) and KK joint (right), and four piles; (b) Top view of the substructure; (c) Front view of the substructure.

Figure 4 .
Figure 4. (a) Schematic of the HydroCRest simulation window; (b) Probability distributions in percentage of the directions of the total force and moment (OTM).

Figure 4 .
Figure 4. (a) Schematic of the HydroCRest simulation window; (b) Probability distributions in percentage of the directions of the total force and moment (OTM).

Figure 5 .
Figure 5. (a) Time-domain wind data reconstructed from power spectral densities (PSDs) based on site data; (b) Time-domain wind data with wind fluctuations added by randomizing the phase shifts.

Figure 5 .
Figure 5. (a) Time-domain wind data reconstructed from power spectral densities (PSDs) based on site data; (b) Time-domain wind data with wind fluctuations added by randomizing the phase shifts.

Figure 6 .
Figure 6.(a) Potential flow tower model with the colour bar indicating the dimensionless pressure coefficient (Cp).(b) Comparison of wind power curves.

Figure 6 .
Figure 6.(a) Potential flow tower model with the colour bar indicating the dimensionless pressure coefficient (C p ).(b) Comparison of wind power curves.

Figure 7 .
Figure 7. (a) Fatigue damage and characteristic fatigue damage (CFD), both dimensionless, of the 60 short-term simulations when the wind speed (U w ) is 7.5 m/s.The uncertainty among the 60 simulations was clear despite simulating under identical weather conditions; (b) Fatigue damage distribution over a specific wind speed.

Energies 2018 ,
11, x 12 of 17condition was smaller than those of the 7.5 m/s and 12.5 m/s wind load conditions, the former case still requires the highest amount of cases to converge within a 5% error.

Figure 8 .
Figure 8. Simulation deviation (error) for the wind-load-alone cases.The graph indicates that the error curve converges as the number of simulations increases.

Figure 8 .
Figure 8. Simulation deviation (error) for the wind-load-alone cases.The graph indicates that the error curve converges as the number of simulations increases.

Figure 9 .
Figure 9. Simulation deviations (error) for the three wave-load-alone cases.

Figure 9 .
Figure 9. Simulation deviations (error) for the three wave-load-alone cases.

Energies 2018, 11 , x 14 of 17 Figure 10 .
Figure 10.CFD results of the grid-based and Monte Carlo (MC) methods with 95% confidence interval for a hot spot of the (a) K and (b) KK joints.

Figure 10 .
Figure 10.CFD results of the grid-based and Monte Carlo (MC) methods with 95% confidence interval for a hot spot of the (a) K and (b) KK joints.

Figure 11 .
Figure 11.CFD error estimation of the grid-based and MC methods for (a) the K joint and (b) the KK joint when an average of two grid-based results (grid average) was selected as the data.

Table 1 .
Extreme and average values, relative standard deviation (CV), and the deviation of the fatigue damage for the wind-load-alone cases.

Table 2 .
[9]reme and average values, CV, and deviation of the fatigue damage for the wave-loadalone cases.Gaussian sea surface, the following simplification of ocean surface wave elevation is used extensively in offshore engineering[9]:

Table 2 .
Extreme and average values, CV, and deviation of the fatigue damage for the wave-load-alone cases.