Comparative Study of Structural Reliability Assessment Methods for O ﬀ shore Wind Turbine Jacket Support Structures

: O ﬀ shore wind turbines (OWTs) are deployed in harsh environments often characterized by highly stochastic loads and resistance properties, thus necessitating the need for structural reliability assessment (SRA) to account for such uncertainties systematically. In this work, the SRA of an OWT jacket-type support structure is conducted, applying two stochastic methods to predict the safety level of the structure considering various design constraints. The ﬁrst method refers to a commercial ﬁnite element analysis (FEA) package (DesignXplorer © from ANSYS) which employs direct simulations and the six sigma analysis function applying Latin hypercube sampling (LHS) to predict the probability of failure. The second method develops a non-intrusive formulation which maps the response of the structure through a ﬁnite number of simulations to develop a response surface, and then employs ﬁrst-order reliability methods (FORM) to evaluate the reliability index and, subsequently, the probability of failure. In this analysis, ﬁve design constraints were considered: stress, fatigue, deformation, buckling, and vibration. The two methods were applied to a baseline 10-MW OWT jacket-type support structure to identify critical components. The results revealed that, for the inherent stochastic conditions, the structural components can safely withstand such conditions, as the reliability index values were found acceptable when compared with allowable values from design standards. The reliability assessment results revealed that the fatigue performance is the design-driving criterion for structural components of OWT support structures. While there was good agreement in the safety index values predicted by both methods, a limitation of the direct simulation method is in its requirement for a prohibitively large number of simulations to estimate the very low probabilities of failure in the deformation and buckling constraint cases. This limitation can be overcome through the non-intrusive formulation presented in this work.


Introduction
In line with the quest to battle against the energy crisis and global climate change, wind energy is found to provide a technically and economically feasible solution, making it the most promising renewable energy technology in Europe and beyond (i.e., China, United States of America (USA), etc.) [1]. Today, the vast majority of wind energy is generated from onshore wind farms. However, there is a limitation in the further growth of onshore wind farms as a result of visual impact and limited space available to deploy wind turbines. In comparison to land, there is more available space to deploy wind turbines at sea, and the wind shear is stronger at offshore locations, driving the wind industry to move offshore, currently aiming for offshore wind in Europe to reach 64.8 GW by 2030 [1][2][3][4].
With few exceptions, most of the existing offshore wind turbines (OWTs) are installed on monopile foundations in shallow waters (depths < 50 m). However, as turbines become larger (>5 MW) and as water depth increases beyond 35-50 m, monopiles become progressively uneconomical, primarily because of modal requirements that force structural dimensions to grow even beyond fabrication and installation capabilities currently available [5]. State-of-the-art research on monopile foundations can be found in References [6,7]. Space frame structures such as jackets, widely used in the petroleum industry, offer a lighter and yet stiff alternative to monopiles. Effective design of these structures from a turbine system dynamics perspective is a demanding task. Although jacket-type support structures could contribute to further developments in the offshore wind industry, research is still required to support their basic design and analysis [5,8], as well as to take into account the requirement for optimization due to both the serial production required and the nature of these unmanned structures. As a result of application-based research, several new methods and patents were developed aiming to resolve some of these issues [9][10][11][12][13][14][15][16][17][18].
In a conventional structural analysis, the environmental loads, geometric dimensions, material properties, and parameters are treated as deterministic quantities, which are often expressed in the design through characteristic values incorporating partial safety factors for consideration of uncertainties. This type of analysis is adequate in cases where randomness is relatively small. However, structures used for offshore deployment are generally complex and subjected to high uncertainties introduced by the environment, manufacturing, and operational processes. Randomness in design parameters for such systems should be systematically accounted for via stochastic modeling, allowing for a better informed and potentially optimized design that would take into consideration the life cycle of the structure and associated time-variant mechanisms of damage. Reliability analysis enables an efficient approach to better understand the response of the system to variation of an input parameter, enabling the systematic modeling of uncertainties occurring from highly complex models, loading conditions, etc. [19]. In particular, for offshore jacket foundations, soil-structure interactions and particularities of marine environments introduce considerable uncertainties (e.g., soil properties, wind, wave, and current loads, etc.), leading to the requirement of employing proper reliability analysis methods. A structural reliability model of OWT support structures requires two key components, i.e., (1) a structural model describing the structural behavior of OWT support structures, and (2) a suitable reliability assessment method which finds the reliability index β (or the probability of failure P f ) of critical members with respect to the constraints [20].
Structural models for OWT support structures can be roughly classified into two groups, i.e., one-dimensional (1D) beam models and three-dimensional (3D) finite element analysis (FEA) models.
In the 1D beam model, the support structure is discretized into a sequence of elastic Euler or Timoshenko beam elements. It is computationally efficient and capable of accurately modeling global structural behavior, such as deflections and modal frequencies [21]. However, it is incapable of accurately representing the structural responses at a local scale such as stress concentration effects [22]. In the 3D FEA model, the jacket foundations are usually defined by applying brick or shell elements. In comparison to the one-dimensional beam model, the three-dimensional finite element model has the capability of examining detailed stress distributions and accurately capturing structural responses within the structure. The three-dimensional finite element analysis model was extensively applied to Box-Behnken design DOE type, variational technology, and Monte Carlo simulations on the response surface [53,54]. More information on these probabilistic methods can be found in References [48,54]. The benefit of the tool is in its user-friendly interface and direct link with FEM modeling. However, it relies on direct simulations, which is a limitation for the order of magnitude probabilities of interest for offshore and marine applications, and which cannot be generalized and integrated with more specialized tools for more advanced analysis.
The second approach, which was developed by the authors, proposes a non-intrusive stochastic formulation, which links the structural model to a reliability analysis algorithm through approximation modeling. In this case, quadratic response surface modeling is adopted in order to map the response of the structure in the domain of stochastic variables. The benefit of this approach is that it allows for high-fidelity computational tools to be employed for the analysis in a modular way, thereby extending its applicability to various engineering problems requiring advanced modeling techniques, and further permitting coupling with analytical reliability methods such as FORM and SORM, allowing for low values of probability of failure to be calculated.
This work aims to carry out a comparison of the two methods described above for a typical OWT jacket-type support structure considering various design constraints and stochastic variables. This paper is the first of its kind to study the structural reliability of an OWT jacket-type support structure using two methods and then comparing the results, as there is no work which reported the use of the post-processing tool inbuilt in ANSYS applied to an OWT jacket-type support structure to determine the probability of failure, nor any work comparing results from the probabilistic FEA in ANSYS to a non-intrusive stochastic formulation, which is also documented in this paper. Furthermore, the proposed non-intrusive method is verified against the commercial package, allowing for further complex engineering problems to be investigated without the need for a fully integrated application-specific tool, rather than through a modular approach which can allow the use of simpler and more computational efficient tools.

Monte Carlo Simulation
As the name implies, Monte Carlo simulation techniques involve "sampling" at "random" to simulate artificially a large number of experiments and to observe the result. In the case of analysis for structural reliability, this means, in the simplest approach, sampling each random variable X i randomly to give a sample valuex i . The limit state function G(x) = 0 is then checked using the sample set of valuesx i . If the limit state function is violated (i.e., G(x i ) ≤ 0), the structure or structural element "fails". The experiment is repeated many times, each time with a randomly chosen vectorx ofx i values. If N trials are conducted, the probability of failure is given approximately as follows [55]: where n(G(x i ) ≤ 0) is the number of trials n for which (G(x i ) ≤ 0). Obviously, the number N of trials required is related to the desired accuracy for p f . It is clear that, in the Monte Carlo method, a game of chance is constructed from known probabilistic properties in order to solve the problem many times over, and from that to deduce the required result (i.e., the failure probability). The technique sketched above is the simplest Monte Carlo approach for reliability problems but not the most efficient. The basis for its application is outlined below. The probability of limit state violation may be expressed as follows [56]: where I[] is an "indicator function" which is equal to one if [ ] is "true" and zero if [ ] is "false". Thus, the indicator function identifies the integration domain. Equation (1) represents the expected value of I[]. If x j represents the j-th vector of random observations from f x (), it follows directly from sample statistics that Equation (3), is an unbiased estimator of J. Thus, Equation (3) provides a direct estimate of p f . In exploiting this procedure, three matters are of interest: how to extract most information from the simulation points, how many simulation points are needed for a given accuracy, or, conversely, how to improve the sampling technique to obtain greater accuracy for the same or fewer sample points [55]. The fundamentals of the Monte Carlo simulation method are well documented in literature [48,57]. The key functionality of Monte Carlo simulation techniques is the generation of random numbers with a uniform distribution from 0-1 [54]. The DesignXplorer© uses the advanced Mersenne Twister (MT) algorithm [58]. The MT algorithm is adopted for generating uniform pseudorandom numbers. In order to accommodate Mersenne-prime period, the MT is a modified variant of a previous generator, the twisted generalized feedback shift register (TGFSR). The MT generator was implemented in portable C-code and passed stringent statistical tests, including diehard. Its merits are due to the efficient algorithms that are unique to polynomial calculations over the two-element field [58]. ANSYS also uses the inverse probability method [48] to generate random numbers with arbitrary distributions. The implementation of the inverse cumulative probability function was verified with the results published by Kececioglu [59] with matching results even for probabilities as low as 10 −24 . It uses the Latin hypercube sampling technique [60,61] to make the sampling process more efficient, as well as to ensure that the tails of the distribution of the input variables are better represented.
The interpretation of the results of a Monte Carlo simulation analysis is based on statistical methods. The statistical procedures to calculate, for example, mean values and standard deviations, and to derive histogram plots are available in textbooks on statistics and probabilistic methods (e.g., Reference [62]). For the cumulative distribution function (CDF), the data are sorted in ascending order, and the CDF of the i-th data point, here denoted with F i , can be derived as follows [54]: In DesignXplorer©, this equation is solved numerically for F i . From the CDF curve, failure probabilities, as well as the inverse probability, can be derived.
In ANSYS, sensitivities are evaluated based on correlation coefficients [62]. The stronger a result parameter is correlated with a particular input variable, the more sensitive the result parameter is with respect to the input variable.

Accuracy and Validity
Monte Carlo simulation methods have two major sources of inaccuracy. Of course, the number of samples cannot be infinite, but must be limited. The error due to the limitation of the number of samples can be easily quantified. As described in Reference [62], confidence intervals from statistical results of a Monte Carlo simulation can be easily obtained. In the context of finite element analysis, another source of error, which is typically overlooked, is the error due to re-meshing the geometry if geometric uncertainties are included in the probabilistic model. Furthermore, minute changes in the geometry can cause the mesh density in critical areas of the geometry to change in a non-continuous fashion. This leads to a change in the finite element results much larger than would be expected for a small geometry change. This error due to re-meshing cannot be easily quantified, but it can be minimized, if mesh controls are applied such as specifying the number of mesh divisions in critical areas of the geometry. However, this is a manual and often time-consuming operation. A method avoiding the influence of the meshing error involves variational technology.

Advantages and Disadvantages
The Monte Carlo simulation method does not make any simplification or assumptions in the deterministic or probabilistic model. The only assumption it does in fact make is that the limited number of samples is representative to quantify the randomness of the result parameters. The error associated with this assumption is well quantifiable as outlined above. With an increasing number of samples, the Monte Carlo simulation method converges to the true and correct probabilistic result. The Monte Carlo simulation is, therefore, widely used as the benchmark to verify the accuracy of other probabilistic methods. Another advantage is the fact that the required number of simulations is not a function of the number of input variables.
The disadvantage of the Monte Carlo simulation method is its computational cost. As a rule of thumb, addressing low probabilities of failure requires about N sim ≈ 100 P f simulation loops, where P f is the targeted failure probability. Hence, if the targeted failure probability is low, then the required number of samples may be prohibitively large, making the Monte Carlo simulation method impractical for real engineering problems.

Variance Reduction
According to Reference [55], the variance σ 2 directly affects the variance of J 1 , and the number of samples N inversely affects the variance of J 1 . This means that the standard deviation ofĵ and, hence, of the Monte Carlo estimate decreases in proportion to N − 1 2 . By comparison, for one-dimensional integration, the error in standard deviation reduces as N −2 for the trapezoidal rule and N −4 for Simpson's rule [63]. Obviously, the slow convergence of the "crude" Monte Carlo method is a significant disadvantage for its practical application. One approach involves ways to reduce the variance σ 2 for the estimate of the probability of failure J. These are the so-called "variance reduction" techniques [64]. Apart from increasing the number of samples, variance reduction can be achieved only by using additional (a priori) information about the problem to be solved. The Latin hypercube sampling technique is an example of various variance reduction techniques. A good overview of the various strategies for variance reduction in general Monte Carlo work was given by References [64,65] in relation to some techniques for structural reliability calculations.

Sample Generation
For SSA, the sample generation is based on the Latin hypercube sampling (LHS) technique by default. In the properties view for a six sigma analysis cell, sampling type can be set to either LHS or WLHS (weighted Latin hypercube sampling). LHS is a more advanced and efficient form of Monte Carlo analysis methods. With LHS, the points are randomly generated in a square grid across the design space, but no two points share input parameters of the same value. This means that no point shares a row or a column of the grid with any other point. Generally, LHS requires 20% to 40% fewer simulations loops than the direct Monte Carlo technique to deliver the same results with the same accuracy. However, that number is largely problem-dependent [53].

Latin Hypercube Sampling
The Latin hypercube method, originally presented in Reference [66], is a method that can represent multiple variables avoiding overlapping datasets. According to Reference [55], stratified and Latin hypercube sampling methods use the common idea that the sampling space be divided into strata (or hypercubes) and that only a few of the many possible samples may be selected in each strata (or hypercube) for inclusion in the estimation of the integral J. Application of the method initiates by dividing the distribution of each stochastic variable in n, non-overlapping intervals with equal probability. For each variable, one value should be randomly selected from each interval, and the analysis point derived from each dataset is then associated. The homogeneous allocation of intervals on the probability distribution function yields relatively small variance in the response, compared to the conventional Monte Carlo sampling. At the same time, the analysis is much less computationally demanding to generate [49].
The samples used from within each stratum (or hypercube) can be randomly selected from those that would otherwise be used; alternatively, for example, the mean or central point might be used as the one representative value. The evaluation of J then uses these selected values, weighted by their probability of occurrence. Latin hypercube sampling considers all the random variables [55,64,67,68].
For convenience, consider the special case in which all the random variables x i (i = 1, . . . , n) are each divided into the same number of strata, say N. It is convenient also to make the divisions such that the probability of x i falling in each stratum is 1/N. Now, select a representative value for each stratum and for each random variable. In principle, more values could be selected, but this can also be achieved by using a greater number of strata. To estimate J, the procedure is to use just one of the N representative values for each random variable x i (that is, selected from one of the strata) in combination with a selected representative value from each of the other random variables. This collection of representative values forms a sample set. The process is repeated, selecting from the remaining representative values, thereby forming another sample set, and so on. In the end, there will be N sample sets. These can be interpreted directly as the number N of sample vectors selected in Monte Carlo sampling and, thus, can be used directly in Equation (1), with evaluations of the limit state function using the present sample sets to estimate J (and, thus, p f ). Latin hypercube sampling can also be combined with other variance reduction techniques [68], such as importance sampling, although it is generally considered to be most suitable for small-scale simulation problems [55,56,69].

ANSYS DesignXplorer© for Structural Reliability Analysis
In order to perform the SRA, a stochastic parametric FEA model is firstly built in the ANSYS environment such that the different input parameters are assigned their corresponding distributions. Then, the developed FEA model is employed in performing a series of FEA simulations on the OWT support structure through the DoE module packaged in the DesignXplorer© facility in ANSYS, in order to map the response domain and then derive an appropriate response surface model. The purpose of the response surface is to interpolate the values in the multiple dimensions characterized by the DoE. There are several types of response surfaces as packaged in DesignXplorer© [53] which include the genetic aggregation, standard response surface full second-order polynomials, kriging algorithms, non-parametric regression, and the sparse grid. However, in this study, only the standard response surface full second-order polynomials with manual refinements is applied, along with the selection of the specification of three verification points. The response surface is interpolated from design points in the DoE. Also packaged within the response surface module is the goodness of fit metric, which is calculated for the DoE points and can be calculated for verification points to check how accurate the response surface can predict the design points. In all limit state cases, the predicted vs. observed chart should be checked to show the goodness of fit data for one or more outputs. Furthermore, the output data should be observed in order to ascertain whether most of the points fall on or near the line, which implies that the response surface predicts the values for most of the design points within its range including the verification points accurately. This indicates that the response surface is a good fit for the DoE points.
In this study, the six sigma analysis function of the DesignXplorer© module in ANSYS is used for the probabilistic assessment. Latin hypercube sampling (LHS) with number of samples equal to 1 × 10 7 is defined and applied to all the design constraint cases. The cumulative distribution function (CDF) is used to assess the probability of failure and, hence, the reliability of the component. The cumulative distribution function value at any given point expresses the probability that the respective parameter value remains below that point. Whether this probability represents the failure probability or the reliability of the component depends on how failure is defined [53]. The corresponding reliability index β is evaluated through appropriate statistical transformation [55].

A Modular Structural Reliability Assessment Algorithm
From the same simulations executed on the ANSYS DoE tool presented in Section 2.3, the results are imported into a MATLAB code that was developed for response surface modeling, which is presented in Section 3.2 [1,2]. Figure 1 depicts the flow chart of structural reliability assessment of the non-intrusive formulation [1,2]. In brief, the analysis starts from the definition of the system, associated limit states, and determination of which variables will be considered as stochastic in the analysis. Next, a number of FEA simulations are executed by varying inputs and recording outputs related to limit states in order to map the response of the structure in the domain of stochastic variables. An appropriate method is then selected in order to approximate the response of the system and link output variables with global inputs (such as wind and wave loads). Finally, the reliability index is calculated using MCS or FORM methods. In this instance, the iterative FORM algorithm is adopted to allow for the calculation of low probabilities of failure.

A Modular Structural Reliability Assessment Algorithm
From the same simulations executed on the ANSYS DoE tool presented in Section 2.3, the results are imported into a MATLAB code that was developed for response surface modeling, which is presented in Section 3.2 [1,2]. Figure 1 depicts the flow chart of structural reliability assessment of the non-intrusive formulation [1,2]. In brief, the analysis starts from the definition of the system, associated limit states, and determination of which variables will be considered as stochastic in the analysis. Next, a number of FEA simulations are executed by varying inputs and recording outputs related to limit states in order to map the response of the structure in the domain of stochastic variables. An appropriate method is then selected in order to approximate the response of the system and link output variables with global inputs (such as wind and wave loads). Finally, the reliability index is calculated using MCS or FORM methods. In this instance, the iterative FORM algorithm is adopted to allow for the calculation of low probabilities of failure.

Multivariate Regression
Regression analysis is a statistical process for establishing the relationship between a dependent variable and one or more independent variables. Taking the ultimate limit state (ULS) stress constraint as an example, the dependent variable (i.e., maximum von Mises stress, , ) and independent variables (i.e., thrust load , tilting moment , torsional moment , total hydrodynamic load (wave and current) , soil 1 Young's modulus , soil 2 Young's modulus , soil 3 Young's modulus , and weight of rotor-nacelle assembly (RNA) ) are assumed to have the following functional relationship:

Multivariate Regression
Regression analysis is a statistical process for establishing the relationship between a dependent variable and one or more independent variables. Taking the ultimate limit state (ULS) stress constraint as an example, the dependent variable (i.e., maximum von Mises stress, σ VM, max ) and independent variables (i.e., thrust load x 1 , tilting moment x 2 , torsional moment x 3 , total hydrodynamic load (wave and current) x 4 , soil 1 Young's modulus x 5 , soil 2 Young's modulus x 6 , soil 3 Young's modulus x 7 , and weight of rotor-nacelle assembly (RNA) x 8 ) are assumed to have the following functional relationship: where (a 0 , a 1 , . . . , a 16 ) are 17 regression coefficients for a quadratic regression. For other types of constraints, i.e., deformation, buckling, vibration, and fatigue, expressions similar to Equation (5) can be derived. Multivariate regression can be used to obtain the regression coefficients for different limit states [49].
Taking the stress constraint case as an example, 300 stochastic FEA simulations were performed, obtaining 300 samples. The R-square value obtained in the regression analysis performed in this case was consistently higher than 0.95, which is relatively high and indicates a satisfactory multivariate regression to be used in this study.

FORM (First-Order Reliability Method)
Having obtained the performance function from regression, FORM is used to calculate the reliability index β. The FORM Hasofer and Lind algorithm is summarized below [70].

1.
Define the performance function for different limit states. For example, in the ULS (stress constraint) case, Equation (5) is substituted into Equation (21), which yields the following performance function: 2.
Set the mean value point as an initial design point, i.e., x i,k = µ x i i = 1, 2, . . . , n, and calculate the gradients ∇g(X k ) of the limit-state function at this design point. Here, x i,k refers to the i-th element in the vector X k of the k-th iteration, and µ x i is the mean value of the i-th element.

3.
Calculate the initial reliability index β using the mean-value method, i.e., β = µ g σ g and its direction cosine α.

4.
Compute a new design point X k and U k , function value, and gradients at this new design point.
Repeat steps 4-5 until the convergence of the reliability index β is achieved.

Sources of Loads
OWTs are exposed to several load sources usually imposed by their environment. In design codes, such as DNV-OS-J101 [71] and IEC 61400-3 [72], loads which are relevant to the design of OWT jacket support structures are outlined. Calculations and formulations of each load produced by the environment in this study are based on DNV-RP-C205 [73]. The loads relevant to OWT support structures can be roughly classified into six groups, i.e., (1) aerodynamic loads transferred from the rotor, (2) wind loads on the tower, (3) inertia loads, (4) current loads, (5) wave loads, and (6) hydrostatic loads [24]. These loads are illustrated in Figure 2.
Gravitational and inertia loads are static and dynamic loads which result from seismic activity, the entire structure's weight and imposed equipment on the structure under operation, and the exerted forces on the structure due to operations [74]. The inertia loads, as a result of the support structure and rotor-nacelle assembly (RNA) weight at the tower top, can have a significant impact on buckling and considerably influence the eigen-frequencies of the OWT support structure [24].
Aerodynamic loads are dynamic and static loads resulting from the airflow interaction with the moving and stationary wind turbine components. The load magnitude varies, and this depends on the air density, the mean wind speed and turbulence across the rotor plane, and the aerodynamic shapes of the components of the wind turbine and their interactive effects. These loads can be computed using aeroelastic models [74]. The aerodynamic loads imposed by the rotor are transmitted to the top of the tower and are usually decomposed into a load matrix defined in the referential axis of the wind turbine. The fatigue loads were determined by employing the damage equivalent load (DEL) method, which was described in Reference [75]. The support structure was subjected to thrust load directed along the diagonal of the jacket base and applied at the assumed RNA center of mass (CM). This loading set-up resembled a yaw error situation and produced an additional torsional moment at the tower top. The presence of a horizontal offset for the RNA CM led to a gravity moment at the tower top [8]. Wave load and rotor thrust vectors were directed along the jacket base diagonal, causing maximum compression and tension on the opposite legs [5].
Environmental loads are those caused by environmental phenomena, i.e., the set produced by wave force, current force, and wind force. Wave loads are dynamic loads produced by waves and their interaction with the OWT support structure. Le Mehaute's graph indicates the validity range of each wave theory and can be employed to determine the most applicable theory for given characteristics of waves, such as wave period (T), wave height (H), and water depth (d) [74]. The regular wave theories are ideal for designs where a single wave method is adopted, and this is a widely used method in offshore structure design. In this case, a regular wave with the appropriate period and height is used to represent an extreme wave. This method provides an easy analysis in estimating the extreme response of offshore structures. On the other hand, the random ocean wave is defined by an energy density spectrum. The ocean wave energy content and its distribution over a random wave frequency range are described by the wave energy spectrum. Therefore, the random wave design method may be ideal for designing floating structures [76]. In this study, the wave profile data presented in Table 1 fall within the area of validity of the Stokes third-order wave theory. Stokes wave theory is one of the most commonly used methods in the analysis of offshore structures because of its accuracy in predicting the kinematic properties of the wave.  A typical bottom-fixed jacket offshore structure is considered to be dominated by drag owing to the wavelength, water depth, structural shape, and size; thus, Morison's equation is adopted. An assumption that the total force due to waves acting on a structure can be computed by the linear superimposition of the drag and inertia forces is made in Morison's equation. For Morison's equation to be applied to determine the wave loads acting on the structure, the member diameter, , has to be less than one-fifth of the wavelength, , i.e., For slender members such as jacket members submerged in water, both current and wave loads  A typical bottom-fixed jacket offshore structure is considered to be dominated by drag owing to the wavelength, water depth, structural shape, and size; thus, Morison's equation is adopted. An assumption that the total force due to waves acting on a structure can be computed by the linear superimposition of the drag and inertia forces is made in Morison's equation. For Morison's equation to be applied to determine the wave loads acting on the structure, the member diameter, D, has to be less than one-fifth of the wavelength, λ, i.e., For slender members such as jacket members submerged in water, both current and wave loads can be determined using Morison's equation [71].
where the drag force is the first term while the inertia force is the second term; ρ w is the density of water; C d and C m are the drag and inertia coefficients of the pile, respectively, and their corresponding values are 1.25 and 2.0, respectively [5]; D is the cylinder diameter; a x and u x are the horizontal wave-induced or current-induced acceleration and velocity of water, respectively, which were derived from the Stokes third-order theory. Current loads are dynamic loads which are produced by the flow of water from external sources, apart from surface waves. Current can cause drag loads on the support structure. The current velocity u c (z) is described by a sub-surface current exponential profile from the mean sea level (MSL) to the seabed, d. The total current velocity of a submerged, stationary, tubular member, in the absence of vortex shedding, can be calculated as follows: where u c,MSL is the current velocity at MSL, and d is the depth of water from MSL to the seabed. The wave particle velocity and the current velocity are summed up in Morison's drag term in Equation (14), assuming that the wave and current are aligned. The wind loads on the tower structure are a result of the drag and a function of the average wind velocity V(z). The wind shear is generally represented by a power law profile, defined by the following equation: where V r denotes the reference wind speed at the elevation of nacelle altitude Z r , and α is the roughness coefficient with a value of 0.115 for the offshore site [24]. Wind loads acting along the tower can be calculated as follows: where C D,T is the coefficient of drag of the tower, with a typical value of 1.0 [24], and D(z) represents the external diameter of tapered tower at height z. Having determined the wind loads, the moments due to the action of the wind loads can then be calculated. Hydrostatic pressure, which acts on the submerged sections of the jacket support structure, is a constant normal load and varies with the depth of water. The hydrostatic force can be expressed as where F h is the hydrostatic force, g is the acceleration due to gravity, and h is the depth of the water.

Site Specifications and Design Load Cases
The design load cases should take into account a set of design situations considering the most critical conditions that an OWT support structure can be subjected to, combining extreme or standard external conditions with operational states of the wind turbine or other operational modes (such as installation, transportation, fault, or maintenance) [74]. For the structural design of OWTs, IEC61400-3 [72] outlines 32 DLCs (design load cases) consisting of all operational modes of an OWT, which include start-up, standard operation, shut down, and 50-year extreme conditions. These DLCs can be roughly categorized into two main groups, namely, ultimate and fatigue. Basically, the typical load cases adopted in the structural design of OWTs are mainly the fatigue load considering normal sea conditions and the ultimate load considering 50-year extreme conditions. Both the fatigue and the ultimate load cases are considered in this study.

Fatigue Load Case
Due to the nature of rotor operation, as well as hydrodynamic loading, a significant source of periodic loadings is imposed during the lifetime of an OWT. OWT support structures are, thus, prone to fatigue failure [24,77]. A widely used fatigue DLC correlates to an operating state within the normal sea state (NSS) and normal turbulence model (NTM), where the site is assumed to have no current and both the wave height and the cross zero period are obtained via a probability density function of the site. DLC 1.2, as prescribed in the IEC standard [72], are generally regarded as the governing fatigue DLCs for OWT support structures; therefore, they are considered in this study as fatigue load cases [24,78]. According to the IEC standard [79,80], the safety factor for fatigue load is equivalent to 1.0.

Ultimate Load Case
For the extreme hydrodynamic loads experienced by the OWT, that produced by the 50-year return period is widely considered as a critical ultimate load case. It was demonstrated from a previous study that wind loading generally governs the design of OWTs rather than hydrodynamic loading [81]. Therefore, the critical load case for ULS is mostly assumed to be defined by the parked turbine, under the 50-year EWM (extreme wind model) with 50-year RWH (reduced wave height) and ECM (extreme current model). The loading characteristics, as described above, correspond to the IEC 61400-3 [72] DLC 6.1b and 2.1 GL regulation [82], respectively. Load safety factors for gravitational and other loads (such as environmental loads) are given as 1.1 and 1.35, respectively [79,80].

Site Specifications
The information available from a previous study for a deep-water site located off the coast of Netherlands designated as the K-13 deep-water site [80] was adopted for this study, as listed in Table 1.

Baseline 10 MW OWT
There is continued interest in increasing the size of offshore turbines given the high cost of the support structure and associated operation and maintenance (O&M) activities. The baseline 10-MW OWT [5,8,83] is an envisioned simple reference 10-MW turbine scaled from a 5-MW NREL (national renewable energy laboratory) reference turbine [84]. The main properties were calculated through basic mechanics considerations and scaling laws. The RNA (rotor-nacelle assembly) is supported by an 88.4-m tapered tubular tower. Its basic characteristics are summarized in Table 2, and further details can be found in References [5,8].

Parametric FEA
A parametric finite element analysis model of the OWT support structure and the soil system was developed using ANSYS software package, a widely used multi-purpose finite element analysis software. The parametric FEA model enables modeling OWT support structures with stochastic variables, such as environmental loads and material properties. The environmental loads, i.e., wind, wave, and current loads, are applied to the structure in an uncoupled way, ignoring the fluid-structure interactions, as illustrated in Figure 2. The parametric FEA model is applied to the jacket support structure of the baseline 10-MW OWT, of which the main parameters are summarized in Table 3. The geometry, material, mesh, and boundary conditions used in the FEA model are presented in the upcoming subsections. Table 3. Baseline 10-MW support structure parameters [5,8,80,83]. MSL-mean sea level; OD-outer diameter; TP-transition piece; RNA-rotor-nacelle assembly; UL-ultimate limit; FL-fatigue limit. The configuration of the jacket is usually a structure with four legs, inter-connected with bracings welded with diameters up to 2 m. A transition piece is connected to the tower bottom and transfers the forces from the tower into the jacket structure, which is anchored via piles or suction caissons to the seabed at each leg. The jacket structure has a vertical height of 66 m, and it is joined with the tower through a 7-m-long transition piece. The transition piece is envisioned as a deck reinforced by stringers, which hosts a central cylindrical shell supported by four tubular struts. The piles are embedded 42 m in the soil and are intended to be driven through the legs and, thus, feature a batter angle [5,8].

Mesh
A crucial step in the coupling of the model to be applied for the FEA is the mesh generation. The ANSYS software package provides a powerful and reliable structure mesh generator capable of developing a mesh, which is consistent in all of the structure with minimal computational requirement.
In order to ensure accuracy of results, a mesh sensitivity study was carried out [1,24]. This is depicted in Table 4 and Figure 3. From the mesh sensitivity study, the equivalent (von Mises) stress converged at the mesh size of 2 m for the soil and 0.5 m for the support structure, corresponding to a total number of elements of 190,222. Therefore, 2 m and 0.5 m are deemed as appropriate mesh sizes for the soil and the support structure, respectively. As can be observed, the maximum displacement and buckling load multiplier are not sensitive to change in mesh size. The generated mesh for the FEA model is depicted in Figure 4. for the soil and the support structure, respectively. As can be observed, the maximum displacement and buckling load multiplier are not sensitive to change in mesh size. The generated mesh for the FEA model is depicted in Figure 4.

Material
The jacket structure is made of S355 structural steel with a density of 8500 kg/m 3 (increased by 8% to cover for secondary steel appurtenances, coatings, and welds that are not covered for in the thickness data of the support structure), Young modulus of 210 GPa, a Poisson's ratio of 0.3, and a

Material
The jacket structure is made of S355 structural steel with a density of 8500 kg/m 3 (increased by 8% to cover for secondary steel appurtenances, coatings, and welds that are not covered for in the thickness data of the support structure), Young modulus of 210 GPa, a Poisson's ratio of 0.3, and a nominal yield strength of 355 MPa. The tower and the transition piece are made of structural steel of the same properties [5,24].

Soil Profile
Other than the OWT support structure, a crucial section of the detailed parametric model comprises the soil-structure interaction, which is a factor of design of significant importance for the structural response of the piles-jacket-tower assembly. The soil modeled in the FEA model is cylindrical in shape, with a diameter of 140 m and height of 50 m and consisting of six different strata (see Table 5). The soil material is modeled using the Drucker-Prager strength linear model [85]. This model is used to calculate the yield strength of each soil strata. The yield strength of the soil, σ ys , according to the Drucker-Prager model, is given by where φ represents the internal friction angle, and c represents the cohesion value. The coefficient of friction C f between the soil and the pile is determined as follows [28]:

Boundary Conditions
The wind and wave loads are applied to the surface of the tower and the surface of support structure submerged into the water, respectively. Additionally, the following three boundary conditions are also defined: (1) the bottom of the soil model is fixed against translation in all directions; (2) lateral boundaries of the soil model are fixed against lateral translation; (3) an augmented Lagrangian formulation-based frictional contact [86,87] is established between the soil and pile with suitable friction coefficients in order to enable soil-structure interaction. See Figure 5. The wind and wave loads are applied to the surface of the tower and the surface of support structure submerged into the water, respectively. Additionally, the following three boundary conditions are also defined: (1) the bottom of the soil model is fixed against translation in all directions; (2) lateral boundaries of the soil model are fixed against lateral translation; (3) an augmented Lagrangian formulation-based frictional contact [86,87] is established between the soil and pile with suitable friction coefficients in order to enable soil-structure interaction. See Figure 5.

Modal Analysis
The modal frequencies of the OWT support structure calculated from the present FEA model are compared against those reported in the literature. The comparison results are presented in Table 6. As can be observed in Table 6, the first two frequencies of both fore-aft and side-to-side bending modes calculated from the present FEA model show good agreement with those reported in References [8,83], with the maximum percentage difference (11.11%) observed for the second side-to-side bending mode. This confirms the validity of the present FEA model. See Figure 6.

Deflections in Static Analysis
This section of the validation aims to examine the deformation behavior of the support structure in a static analysis. In the literature [5,8], the RNA weight and a 3.4-MN thrust load were applied on the tower top. The displacements at the RNA elevation and the base of the tower are compared with reference values, and results are presented in Table 7. Under the loaded condition, the deflections at the RNA and at the tower base are measured with respect to the location of the RNA and the center of the tower base, respectively.
As can be observed from Table 7, the results of the present FEA model match well with those reported in Reference [8], for deflections at both the RNA and the tower base, with a maximum difference of 3.44% observed for the deflection at the tower base. This further confirms the validity of the present FEA model. See Figure 7.

Deflections in Static Analysis
This section of the validation aims to examine the deformation behavior of the support structure in a static analysis. In the literature [5,8], the RNA weight and a 3.4-MN thrust load were applied on the tower top. The displacements at the RNA elevation and the base of the tower are compared with reference values, and results are presented in Table 7. Under the loaded condition, the deflections at the RNA and at the tower base are measured with respect to the location of the RNA and the center of the tower base, respectively. As can be observed from Table 7, the results of the present FEA model match well with those reported in Reference [8], for deflections at both the RNA and the tower base, with a maximum difference of 3.44% observed for the deflection at the tower base. This further confirms the validity of the present FEA model. See Figure 7.

Implementation of Structural Reliability Assessment on the OWT Support Structure
In this section, the structural reliability of the OWT support structure is assessed considering three limit states according to DNV-OS-J101, i.e., (1) ultimate limit state (ULS), (2) fatigue limit state (FLS), and (3) serviceability limit state (SLS). The fully parametric FEA model presented in Section

Implementation of Structural Reliability Assessment on the OWT Support Structure
In this section, the structural reliability of the OWT support structure is assessed considering three limit states according to DNV-OS-J101, i.e., (1) ultimate limit state (ULS), (2) fatigue limit state (FLS), and (3) serviceability limit state (SLS). The fully parametric FEA model presented in Section 5.1 is used to perform FEA modeling taking account of stochastic variables. The results are then post-processed using two methods, i.e., (i) using the probabilistic FEA tool packaged in ANSYS referred to as the DesignXplorer©, and (ii) applying the non-intrusive formulation as described in Section 2 and References [2,49].

Constraints and Design Criteria
In the design of OWT support structures, three limit states should be considered according to DNV-OS-J101, i.e., (1) ULS (ultimate limit state), which accounts for the maximum load-carrying capacity (i.e., buckling and yielding stress), (2) FLS (fatigue limit state), which accounts for failure as a result of cyclic loads, and (3) SLS (serviceability limit state), which accounts for tolerance criteria (i.e., vibrations and deflections) applicable to standard use. Thus, the structural reliability of OWT support structures in this study takes into consideration five design criteria, i.e., stress (ULS), vibration (SLS), buckling (ULS), deformation (SLS), and fatigue (FLS).

Stress Constraint
The ULS (ultimate limit state) defines the ability of the structure to resist yielding. In terms of ultimate limit state, the maximum stress in the support structure σ VM, max should not exceed the allowable stress limits σ VM, allow . The limit state function for the von Mises criterion can be expressed as follows: The allowable stress σ VM, allow can be expressed as follows: where σ y is the yield strength, with a value of 355 MPa for steel S355, and γ m is the material safety factor, with a value of 1.1 suggested by DNV-OS-J101 standard [71]. Thus, the allowable stress σ VM, allow is 323 MPa.

Fatigue Constraint
The fatigue limit state is a crucial phenomenon in OWT support structures as they are subjected to significant cyclic loads. OWT support structures normally have a long service period that may exceed 20 years. This, in conjunction with the inspection intervals, affects the reliability requirement of the jacket structural design. According to the S-N curve method, the number of loading cycles to failure N, can be obtained by where A and m are the intercept and the slope of the S-N curve on the log-log plot, respectively; ∆S is the stress range. The two parameters in Equation (23), i.e., intercept A and slope m, are generally given by design standards, e.g., DNVGL-ST-0126 [29]. The performance function of fatigue reliability assessment based on the S-N curve method can be expressed as follows: where subscripts f and SN denote the fatigue limit state and S-N curve method, respectively; N is the number of loading cycles to failure given by Equation (23), and N t is the number of loading cycles expected during the given design life. The number of cycles during the design life N t may be determined as a function of rated rotor speed n rated and availability η a (98.5%) on the location chosen [24,88]. Therefore, considering a design life of 20 years, the number of cycles can be expressed as follows: The minimum fatigue safety ratio f sr, min must be greater than the allowable fatigue safety ratio f sr,allow , which is equivalent to the product of one and the material PSF (partial safety factor) γ m, f for fatigue.
f sr,min ≥ f sr,allow .
Since the material's PSF for fatigue limit state is 1.15 [29], f sr,allow is equivalent to 1.15. Hence, the limit state function based on the fatigue safety ratio can be expressed as follows:

Deformation Constraint
Excessive deflections influence the serviceability of OWT support structures and, therefore, should be avoided. The allowable deflection d allow must exceed the maximum deflection d max to ensure overall structural stability. This can be expressed as The limit state function for deflection criteria can be expressed as The allowable deflection d allow is given by the following empirical equation as suggested by DNV-OS-J101 [71]: where L is the height of the jacket.

Buckling Constraint
The slenderness of the support structure coupled with the large RNA mass at the tower top forces the investigation of the risk of instability as a result of buckling. The ULS static analysis results are employed as pre-stress loads. The load multiplier L m , which is defined as the critical load divided by the present load applied, should exceed the permissible load multiplier, L m,allow in order to avoid failure. This design criterion is given as According to the DNV standard [29], 1.4 is adopted as the L m,allow value in this study. The limit state function for buckling criteria can, thus, be expressed as

Vibration Constraint
One of the critical concerns for OWTs on jacket support structures is the resonance phenomenon. In order to guard against such a phenomenon, the first eigen-frequency f 1st of the OWT support structure must be separated sufficiently from rotor induced frequencies f 1P and blade-passing frequency f 3P .
A soft-stiff structural design, whereby the modal frequency lies between the rotor upper bound f 1P ( f 1PH ) and lower bound f 3P ( f 3PL ) frequencies, is presently the most usual and cost-effective design for a jacket. Therefore, the resonance constraint should be expressed as follows [5,24,89]: 0.133 ≤ f 1st ≤ 0.233.
The limit state function for the resonance criterion at the rotor upper bound can, thus, be expressed as On the other hand, it is expressed at the rotor upper bound as follows: The lower g-function value is chosen to be applied in the reliability assessment.

Stochastic Variables and FEA
The stochastic variables considered in this study are presented in Table 8. Eight stochastic variables were considered for all cases (i.e., wind thrust, tilting moment, torsional moment, the weight of RNA, hydrodynamic load, and three Young's moduli were considered at three soil strata). The total hydrodynamic load is considered differently, based on its dependence on the wave height, wave period, and current speed. To find the peak hydrodynamic force acting on the system, a Weibull distribution is applied to the significant wave height and the current speed, while the lognormal distribution is applied to the wave period [2,90,91]. The distribution is then fitted into a Weibull-equivalent distribution. In the presence of observed data, distribution fitting algorithms (such as Akaike information criterion, Bayesian information criterion, and Kolmogorov-Smirnov) are often employed to determine the shape coefficients of the most appropriate statistical distributions [19]. The coefficient of variation (COV) of all stochastic variables was assumed to be 0.1 except for the stochastic variables with Weibull distribution. The mean values of these stochastic variables are presented in Table 2, Table 4, and Table 5. This is summarized in Table 9, which depicts the means and the corresponding standard deviations. Having defined the stochastic variables, the FEA model presented in Section 5.1 is then used to perform FEA simulations of the OWT support structure with the help of the design of experiments (DoE) module in ANSYS. It enables the input parameters to be designated as stochastic parameters, having different types of distributions. DoE indicates the variation of design points with the output parameter.

Results and Discussion
Figure 8a-e depict the probability density function (PDF) histogram and cumulative density function (CDF) of the output parameters for each design constraint from DesignXplorer©.
For the stress constraint case (ULS), the probability that maximum equivalent stress remains less than σ VM, allow = 323 MPa, P f (σ VM ≤ 323MPa), is about 0.99999867 which means that there is 1 − 0.99999867 = 1.327655 × 10 −6 probability that the von Mises stress exceed the limit value of 323 MPa, i.e., P f (σ VM > 323MPa) = 1.327655 × 10 −6 . Therefore, the probability of failure of the critical component of the structure is 1.327655 × 10 −6 considering the ULS (stress constraint) design criteria. The corresponding reliability index value is β = 4.7.
For the fatigue constraint case (FLS), the probability that the minimum safety factor is less than f sr,allow = 1.15, P f ( f sr ≤ 1.15), is about 5.6461 × 10 −7 . Therefore, the probability of failure of the critical component of the structure is 5.6461 × 10 −7 considering FLS design criteria (using time-invariant analysis). The corresponding reliability index value is β = 4.87.
For the deformation constraint (SLS), the DesignXplorer© simulation could not produce results within the region of P f (d ≤ 0.54), i.e., d allow = 0.54. The value P f (d > 0.54), which is the value of the probability of failure, is expected to be very small, and this is identified as a limitation of the direct MCS requiring a prohibitively large number of simulations to estimate results of this magnitude. Similarly, for the buckling constraint case (ULS), the simulation could not produce the results within the region of P f (L m ≤ 1.4), as this value is expected also to be very small.    The results obtained from the non-intrusive formulation for each design constraint are depicted in Figure 9 and Table 10. The results obtained from a time-invariant fatigue reliability assessment The results obtained from the non-intrusive formulation for each design constraint are depicted in Figure 9 and Table 10. The results obtained from a time-invariant fatigue reliability assessment employing the FORM are the reliability index, β = 5.077, and probability of failure, P f = 0.192 × 10 −6 .   As evident from the results above, the reliability assessment performed on the structure shows that the present model, as designed and for the modeling of the stochastic variables considered, satisfies the recommended reliability assessment criteria, as the reliability indices for all the design constraints considered are within the design thresholds. The target reliability index for OWT support structures is typically 3.71, corresponding to a probability of failure of 10 −4 according to the DNV [91].

Results and Discussion
In this work, the time-variant fatigue reliability assessment was carried out using the fatigue analysis guidelines in the DNV standard [71,92], and this is presented in Figure 10. In this study, the thickness-corrected cathodic-protected D curve given by DNV-OS-J101 [71] is chosen for fatigue analysis. The intercept (A) and slope (m) of the S-N curve used for studying the fatigue life of steel structure in seawater for > 10 6 are given as 15.606 and 5, respectively. From the time-variant fatigue reliability index curve, it can be observed that the structure maintains a reliability index of 3.92, exceeding the defined threshold of target reliability as specified by DNV for the nominal 20 years of operation. It can be observed that the fatigue limit state dominates the design of OWT support structures, as the reliability index for the fatigue limit state is lower than other reliability index values (see the time-variant fatigue reliability assessment at time = 20 years).  As evident from the results above, the reliability assessment performed on the structure shows that the present model, as designed and for the modeling of the stochastic variables considered, satisfies the recommended reliability assessment criteria, as the reliability indices for all the design constraints considered are within the design thresholds. The target reliability index for OWT support structures is typically 3.71, corresponding to a probability of failure of 10 −4 according to the DNV [91].
In this work, the time-variant fatigue reliability assessment was carried out using the fatigue analysis guidelines in the DNV standard [71,92], and this is presented in Figure 10. In this study, the thickness-corrected cathodic-protected D curve given by DNV-OS-J101 [71] is chosen for fatigue analysis. The intercept (A) and slope (m) of the S-N curve used for studying the fatigue life of steel structure in seawater for N > 10 6 are given as 15.606 and 5, respectively. From the time-variant fatigue reliability index curve, it can be observed that the structure maintains a reliability index of 3.92, exceeding the defined threshold of target reliability as specified by DNV for the nominal 20 years of operation. It can be observed that the fatigue limit state dominates the design of OWT support structures, as the reliability index for the fatigue limit state is lower than other reliability index values (see the time-variant fatigue reliability assessment at time = 20 years). As observed in Table 11, the values of safety index, , predicted by the non-intrusive reliability analysis formulation were marginally higher than those predicted by the ANSYS probabilistic analysis for the different design constraints. However, the difference in the predicted reliability index values from both methods remains consistently below 5%; hence, the agreement of the two methods can be considered sufficient. Furthermore, considering the observed time required, the FORM simulations run instantly compared to the LHS (as packaged in the ANSYS DesignXplorer©), which requires about 15 min to run on a high-performance computer (i.e., for the 1 × 10 7 samples).
Additionally, a case study is performed to investigate the effects of COV of the stochastic variables on the fatigue reliability index. In this case, two additional values of COV are considered, i.e., a 5% increase and 5% decrease. The calculated fatigue reliability index over 20-year service life with different values of COV is depicted in Figure 11. As can be seen from Figure 11, the reliability index is sensitive to the value of COV. The higher value of COV means higher uncertainties in stochastic variables, resulting in lower reliability.  As observed in Table 11, the values of safety index, β, predicted by the non-intrusive reliability analysis formulation were marginally higher than those predicted by the ANSYS probabilistic analysis for the different design constraints. However, the difference in the predicted reliability index values from both methods remains consistently below 5%; hence, the agreement of the two methods can be considered sufficient. Furthermore, considering the observed time required, the FORM simulations run instantly compared to the LHS (as packaged in the ANSYS DesignXplorer©), which requires about 15 min to run on a high-performance computer (i.e., for the 1 × 10 7 samples).
Additionally, a case study is performed to investigate the effects of COV of the stochastic variables on the fatigue reliability index. In this case, two additional values of COV are considered, i.e., a 5% increase and 5% decrease. The calculated fatigue reliability index over 20-year service life with different values of COV is depicted in Figure 11. As can be seen from Figure 11, the reliability index is sensitive to the value of COV. The higher value of COV means higher uncertainties in stochastic variables, resulting in lower reliability.  A global sensitivity analysis is also presented from the ANSYS probabilistic analysis for each of the design constraints through post-processing the results. These depict the most critical and important parameters that can radically change the response and reliability performance of the structure. Figure 12a-e show how sensitive the input design variables are to the structure based on the different design criteria. The bar heights indicate sensitivity in these charts.
In Figure 12a, for the probability distribution of the maximum equivalent von Mises stress, it can be observed that the thrust magnitude and the torsional moment are positively correlated with direct impact. The RNA weight and the hydrodynamic load are negatively correlated and have little influence but inversely. The Young's modulus of the soil system and the tilting moment are neutral and have limited or no impact.
In Figure 12b, for the probability distribution of the deformation at the jacket top, it can be observed that four variables are positively correlated with direct impact. These are the thrust load, hydrodynamic load, and tilting and torsional moments, which have significant influence. The Young's modulus of the soil system has little or no impact. The RNA weight is negatively correlated and has a significant impact.
In Figure 12c, for the probability distribution of the minimum output parameter safety factor, it can be observed that Young's modulus of the soil system is neutral and has little or no impact. The RNA weight and hydrodynamic load are positively correlated with direct impact; the RNA weight has significant influence while the hydrodynamic load has quite some impact. The thrust and the torsional moments are negatively correlated and have significant influence but inversely. The tilting moment has very little or no impact. A global sensitivity analysis is also presented from the ANSYS probabilistic analysis for each of the design constraints through post-processing the results. These depict the most critical and important parameters that can radically change the response and reliability performance of the structure. Figure 12a-e show how sensitive the input design variables are to the structure based on the different design criteria. The bar heights indicate sensitivity in these charts.
In Figure 12a, for the probability distribution of the maximum equivalent von Mises stress, it can be observed that the thrust magnitude and the torsional moment are positively correlated with direct impact. The RNA weight and the hydrodynamic load are negatively correlated and have little influence but inversely. The Young's modulus of the soil system and the tilting moment are neutral and have limited or no impact.
In Figure 12b, for the probability distribution of the deformation at the jacket top, it can be observed that four variables are positively correlated with direct impact. These are the thrust load, hydrodynamic load, and tilting and torsional moments, which have significant influence. The Young's modulus of the soil system has little or no impact. The RNA weight is negatively correlated and has a significant impact.
In Figure 12c, for the probability distribution of the minimum output parameter safety factor, it can be observed that Young's modulus of the soil system is neutral and has little or no impact. The RNA weight and hydrodynamic load are positively correlated with direct impact; the RNA weight has significant influence while the hydrodynamic load has quite some impact. The thrust and the torsional moments are negatively correlated and have significant influence but inversely. The tilting moment has very little or no impact. Appl. Sci. 2020, 10 As can be seen from Figure 12d, the RNA weight and the thrust load have the highest impact and are positively and negatively correlated, respectively, with the other variables having little or no impact on the load multiplier probability distribution. It can also be observed that, in Figure 12e, the structure self-weight has the highest direct impact, with the thrust load and the second soil strata having very little inverse impact on the distribution of the output parameter, the first eigenfrequency, and, hence, the reliability of the structure.

Conclusions
In this study, a comparative analysis on the structural reliability of an OWT support structure between a commercial package and a purpose-developed, non-intrusive stochastic formulation was presented. Both methods start by defining the limit states taking account of the various design constraints. The design constraints considered in this study were as follows: stress (ULS), fatigue (FLS), deformation (SLS), buckling (ULS), and vibration (SLS). A 3D (three-dimensional) parametric FEA (finite element analysis) model of a typical OWT jacket-type support structure was developed, taking account of soil-structure interaction and stochastic variables (i.e., soil properties, wind, wave, and current loads). The series of FEA simulation results were post-processed using two methods. In the first method, the design of experiment and response surface module in the DesignXplorer© were used in order to map the response domain and then derive an appropriate response surface model. The probabilistic assessment was performed using the six sigma analysis function employing the Latin hypercube sampling method in computing the probability of failure, which was then translated to a reliability index. In the second method, multivariate quadratic regression was employed to derive the performance function expressed in terms of global stochastic variables using results from the FEA simulations. Next, FORM was employed to calculate the reliability index and then the probability of failure.
The following conclusions can be drawn from the present study:


The two methods, i.e., the probabilistic tool package in ANSYS and the developed non-intrusive formulation, match well on the calculated reliability indices calculated for all limit states, crossverifying applicability of both methods for the structural reliability assessment of OWT support structures.  The proposed non-intrusive formulation for SRA allows faster calculations, assessment of timevariant fatigue reliability assessment, and calculation of low values of probability of failure in complex engineering systems such as offshore and marine structures.  This study also allowed for verification of the non-intrusive formulation, which can allow multiple problems of complex nature to be solved probabilistically, alleviating the need for an As can be seen from Figure 12d, the RNA weight and the thrust load have the highest impact and are positively and negatively correlated, respectively, with the other variables having little or no impact on the load multiplier probability distribution. It can also be observed that, in Figure 12e, the structure self-weight has the highest direct impact, with the thrust load and the second soil strata having very little inverse impact on the distribution of the output parameter, the first eigen-frequency, and, hence, the reliability of the structure.

Conclusions
In this study, a comparative analysis on the structural reliability of an OWT support structure between a commercial package and a purpose-developed, non-intrusive stochastic formulation was presented. Both methods start by defining the limit states taking account of the various design constraints. The design constraints considered in this study were as follows: stress (ULS), fatigue (FLS), deformation (SLS), buckling (ULS), and vibration (SLS). A 3D (three-dimensional) parametric FEA (finite element analysis) model of a typical OWT jacket-type support structure was developed, taking account of soil-structure interaction and stochastic variables (i.e., soil properties, wind, wave, and current loads). The series of FEA simulation results were post-processed using two methods. In the first method, the design of experiment and response surface module in the DesignXplorer© were used in order to map the response domain and then derive an appropriate response surface model. The probabilistic assessment was performed using the six sigma analysis function employing the Latin hypercube sampling method in computing the probability of failure, which was then translated to a reliability index. In the second method, multivariate quadratic regression was employed to derive the performance function expressed in terms of global stochastic variables using results from the FEA simulations. Next, FORM was employed to calculate the reliability index and then the probability of failure.
The following conclusions can be drawn from the present study: • The two methods, i.e., the probabilistic tool package in ANSYS and the developed non-intrusive formulation, match well on the calculated reliability indices calculated for all limit states, cross-verifying applicability of both methods for the structural reliability assessment of OWT support structures.

•
The proposed non-intrusive formulation for SRA allows faster calculations, assessment of time-variant fatigue reliability assessment, and calculation of low values of probability of failure in complex engineering systems such as offshore and marine structures.
• This study also allowed for verification of the non-intrusive formulation, which can allow multiple problems of complex nature to be solved probabilistically, alleviating the need for an integrated, application-specific tool to be developed, instead of with the proposed modular formulation.
Author Contributions: Writing-original draft preparation, methodology, software, validation, and formal analysis, A.A.S.; methodology, software, validation, formal analysis, and writing-review and editing, L.W.; methodology and writing-review and editing, A.M. and K.S.; supervision, methodology, software, validation, formal analysis, funding acquisition, conceptualization, and writing-review and editing, A.K. All authors have read and agreed to the published version of the manuscript.