PROBABILISTIC FINITE ELEMENT STIFFNESS OF A LATERALLY LOADED MONOPILE BASED ON AN IMPROVED ASYMPTOTIC SAMPLING METHOD

The mechanical responses of an offshore monopile foundation mounted in over-consolidated clay are calculated by employing a stochastic approach where a nonlinear p–y curve is incorporated with a finite element scheme. The random field theory is applied to represent a spatial variation for undrained shear strength of clay. Normal and Sobol sampling are employed to provide the asymptotic sampling method to generate the probability distribution of the foundation stiffnesses. Monte Carlo simulation is used as a benchmark. Asymptotic sampling accompanied with Sobol quasi random sampling demonstrates an efficient method for estimating the probability distribution of stiffnesses for the offshore monopile foundation.


Introduction
Soil is a heterogeneous material. Thus the seabed consists of various types of soil and the properties are subject to uncertainty and variation with position.
In heterogeneous materials, failure occurs through the weaker parts. This leads to different types of failure mechanisms in comparison with failure mechanisms in homogeneous soil. In the deterministic design practices, expensive geotechnical in-situ tests are carried out and only a quantile value of them is used as design parameter. In addition, the reliability of the structure may be unknown due to existing uncertainties related to the soil. Indeed, stochastic analysis provides a reliability-based design with fewer costs for tests and construction as well.
A monopile is one of the most common foundations for offshore wind turbines. This type of foundation is cost-effective where the water depth is less than 30 meters. In the new design procedure for offshore foundations, fatigue is considered as the main failure mode. Hence, the first natural frequency of the structure is an important parameter. In this regard, the stiffness of a foundation has a significant role. In the current design methods, a deterministic value of the foundation stiffness is considered by analysing a soil-structure system with deterministic properties for the materials. However, in reality the foundation stiffnesses are random with a probability distribution and cannot be fixed on a deterministic value due to various statistical and physical uncertainties related to material properties or modelling.
The primary focus of this paper is the uncertainties regarding soil properties and developing a method by which the probability distribution of the foundation stiffnesses can be assessed in an efficient manner. The uncertainties related to the model and the structural properties will not be considered in this study.
The p-y method was proposed about 55 years ago to model the soil as a system of uncoupled lateral springs (Reese, Matlock 1956;McClelland, Focht 1958). In this method, the lateral capacity of the pile, modelled as a beam, is determined in terms of displacement, bending moment and shear force by introducing the soil around the beam as uncoupled lateral springs known as a Winkler model. The p-y curves are the load-displacement relationships which represent the interaction between soil and pile via the springs (Reese et al. 1974;Bowles 1988;API 1993). A standard procedure for the design of laterally loaded piles proposed in the relevant guidelines for offshore engineering is the p-y method, cf., e.g. (API 2000;DNV 2004). Numerical methods such as finite-element or finite-deference analyses are commonly used for analysing the pile response. Kim and Jeong (2011) applied the three-dimensional (3-D) finite-element method (FEM) to determine the load distribution and deflection and the soil is assumed to consist of over-consolidated clay. A random field is employed to provide the spatial variation of the undrained shear strength. The foundation stiffnesses derived through the soil-structure interaction analysis are converted into horizontal and rotational springs at the monopile cap. An improved asymptotic sampling method associated with two types of random variable generation (normal and Sobol sampling) is performed to estimate the probability distribution of the stiffnesses. Standard Monte Carlo simulations are used to illustrate the efficiency of the proposed method and its accuracy.

Stiffness of the foundation
The main horizontal loads acting on the wind turbine are induced by the wind and waves. These loads lead to a horizontal force as well as an overturning moment at the base of the turbine tower which is fixed at the monopile cap. In order to obtain the stiffness of the foundation, an equivalent coupled-spring (ECS) model with horizontal and rotational springs has been considered.
A 5 MW offshore wind turbine in 20 m water depth and with a tower height of 80 m has been considered in this analysis. It is assumed that the wind force is providing a quasi-static load due to the mean wind speed so that wave and wind produce an oscillating force around the deformed shape due to mean wind force.
In this study, the pile has been subjected to a quasistatic wind load of 2 MN and a variation of the wind force of ±20% (= ±0.40 MN). Further, the quasi-static value of wave forces is 0 MN with a variation of ±1 MN. A load combination is considered such that horizontal forces at the pile cap are and Thus, variations in the shear force and bending moment are ΔQ = 2.80 MN and ΔM = 120 MNm, assuming that the wave and wind forces act 20 m and 100 m above the seabed, respectively.
The secant stiffnesses of the horizontal and rotational springs are derived from the deflection and rotation at the monopile cap as: (1) Here, Δ represents the changes of the parameter. Moreover, y cap and θ cap are the lateral displacement and rocking rotation of the monopile cap obtained from FEM analysis, whereas Q and M are the external shear force and bending moment at the monopile cap, respectiveliy.

Material properties and geometry
The Euler-Bernoulli beam theory is applied to model a single offshore wind turbine monopile foundation. Table 1 presents the tubular monopile geometry and material properties which are required in the model. The material properties of the soil (moderately over-consolidated clay) are listed in Table 2. The undrained shear strength, c u , is considered as a random field with lognormal distribution. of large diameter piles by lateral load transfer method (p-y curve). The FEM can be applied to a diversity of subjects in geotechnics such as stability analysis of reinforced soil structures (Asaoka et al. 1994), prediction of a safety factor (Ugai et al. 1995), and seismic evaluation Zhang, Kimura 2002). Wakai et al. (1999) applied the 3-D finite-element method in the elasto-plastic regime to investigate the lateral loaddisplacement relationship and the bending-strain distributions along piles. Likewise, Zhang et al. (2000) applied 3-D elaso-plastic FEM to investigate the nonlinear mechanical behaviour of a pile foundation. K. Georgiadis and M. Georgiadis (2012) used the 3-D FEM to analyse the response of laterally loaded piles with the developed p-y curves in the case of lateral loading of piles near the crest of undrained clay slopes.
Reliability based design of the structures subjected to random structural properties or random loads is one of the most interesting approaches, which usually involves estimation of probability distributions of events, especially at the low probabilities. A large number of methods have been established to perform reliability analysis for estimating the probability distributions of rare events. Simulation methods developed for the purpose include standard Monte Carlo (SMC) and its derivatives such as advanced Monte Carlo methods (Melchers 1999;Fishman 1996;Guan, Melchers 2001;Au, Beck 2001;Schuëller et al. 2004;Koutsourelakis et al. 2004;Au 2004;Pradlwarter et al. 2007). Asymptotic sampling (AS) is a type of advanced Monte Carlo method which has been proposed for high-dimensional reliability analysis (Bucher 2009). Sichani et al. (2011a) applied this method for high-dimensional dynamics problem such as wind turbines. In a similar study, asymptotic sampling has been proposed as an efficient method for estimating low, first passage probabilities of high-dimensional nonlinear systems (Sichani et al. 2011b). Several researches have been conducted on stochastic analysis of laterally loaded piles. As an example, Chan and Low (2009) presented a stochastic analysis on a laterally loaded pile with nonlinear soil and pile behaviour. Barakat et al. (1999) studied a reliability-based optimization of laterally loaded piles. Low et al. (2001) employed a stochastic nonlinear p-y method to analyse the responses of laterally loaded piles. Haldar and Babu (2008) performed a study on the effect of soil spatial variability on the response of a laterally loaded pile in undrained clay, where the soil shear strength has been considered as a random field. Andersen et al. (2011) proposed a probabilistic approach to estimate the stiffness and natural frequency of a simple model of a wind turbine founded on a monopile. In the same study by Andersen et al. (2012), the asymptotic sampling method was proposed to estimate rare events of the first natural frequency of an offshore wind turbine on a monopile foundation.
In this study an offshore wind turbine with a monopile foundation is investigated with focus on its mechanical response. A nonlinear p-y curve is integrated into a FEM, In Table 2, μ Ln , σ Ln and Y are the lognormal mean value, lognormal standard deviation and the correlated variables with standard Gaussian distribution, respectively. These lognormal parameters are determined by using a transformation from normal statistical parameters: where: μ is the normal mean value and COV is the coefficient of variation of c u , which are given as Andersen et al. (2012): Thus, µ increases with depth (x). The correlated normal variables Y are calculated based on Cholesky decomposition: Here, is a sequence of uncorrelated standard Gaussian (normal) or quasi (Sobol) random variables, U i , i = 1, 2, …, n, for the nodes of the FEM model. T is the lower triangular matrix calculated by Cholesky decomposition from the covariance matrix (ρ): (5) with the elements of matrix ρ given as: (6) Here ρ ab is an exponential, normalized auto-covariance function as suggested in the geotechnical literature (JCSS 2006). Further, τ ab = |η-ξ| is the absolute distance between points η and ξ in the soil medium measured in the depth direction (i.e. τ = x a -x b ), and d is a correlation parameter considered 5 m, cf. Haldar and Babu (2008). Figure 1 illustrates an example of a one-dimensional random field of undrained shear strengths to be applied into the probabilistic model.

Finite element model of the soil and pile response
The basis of simulation has been carried out by using the FEM. A nonlinear p-y curve as presented by Low et al. (2001) is assumed to present the static response of the elastic pile, and undrained clay is considered. The following explanation describes in detail the steps involved in using the p-y curve and the FEM for analysis of a tubular offshore wind turbine monopile foundation.

Nonlinear p-y curve
The Matlock model (Reese, Matlock 1956) is employed to present a nonlinear p-y curve for clay. According to the model, which is based on the absolute value of the pile deflection, y, the soil resistance can be determined as: where y 50 is the deflection at one-half the ultimate soil resistance. It is determined as: where ε 50 is the strain corresponding to one-half the maximum principal stress difference in the soil. Further, the residual and ultimate soil resistances per unit length of the cylindrical pile, p rest and p ult , are calculated as: where: γ΄ is the submerged unit weight of soil; c u is the undrained shear strength of the soil; z is the depth; D is the pile diameter; x r is the depth at which failure is transmitted from wedge-type to flow-around-pile failure; and J is an empirical dimensionless parameter. More details can be found in the work by Andersen et al. (2012) as well as Reese and Matlock (1956).
It is worth mentioning that the sign of p in Eqn (7) needs to be changed if the pile deflection is negative. Moreover, the minimum value of y is considered 10 -4 (m), which means that if y < 10 -4 m then y =10 -4 m is considered. The secant stiffness k of the soil per unit length along the pile is calculated from the resistance divided by the deflection, i.e.: (10)

Finite-element model of the pile
The simulation has been carried out by utilization of the FEM. The monopile is modelled as an Euler-Bernoulli beam and discretized into finite segments as illustrated in Figure 2. Soil pressure along the monopile is considered as horizontal springs which act at the nodes.
A plane beam element with two nodes and cubic interpolation of the lateral displacement field has been adopted. Each node has two degrees of freedom, i.e. the horizontal translation and the in-plane rotation. The element properties are given in terms of the bending stiffness EI and the element length l elem . The reaction forces from the soil act via horizontal springs at each node of the pile FEM model, and the external loads (shear force and overturning moment) are exerted on the top of pile as shown in Figure 2. The problem is solved by iteration until balance between internal and external forces has been obtained, or until collapse of the soil has been identified, i.e. if convergence fails. Figure 3 shows the responses of the monopile with material properties as listed in Tables 1 and 2 and a deterministic value of c u = 150 kN/m 2 . A shear force of Q = 10 3 kN and a bending moment of c u = 10 5 kN.m 2 has been considered at the pile cap.

Validation of finite element model
For verification of the results presented in this study, the modified finite difference method (FDM) reported by Andersen et al. (2012), based on an FDM model proposed by Low et al. (2001), is used to estimate the probabilistic foundation stiffnesses. A nonlinear p-y curve proposed in Section 3.1 with the load combination in The probability density functions for the horizontal and rotational stiffnesses at the pile cap obtained from 7000 Monte Carlo simulations using the FDM are given in Figures 4 and 5 for a pile in over-consolidated clayey soil with the material properties and geometry as defined in Section 2. It can be seen that the results for histograms and lognormal distributions of this study are in close agreement with those in Andersen et al. (2012). For more details, results of Table 6 in Andersen et al. (2012) can be compared with the results in Table 3. It can be seen that mean values of stiffnesses obtained in

Concept of asymptotic sampling
The original work on the asymptotic sampling (AS) method was proposed by Bucher (2009) in order to analyse high-dimensional reliability. The basic idea of AS is to enforce more samples in the target domain by increasing excitation power so that more samples in the low-probability domain are obtained. Based on the asymptotic estimation of small probabilities and the scaling of them, the target probability is estimated (Bucher 2009). In this regard, the standard deviations of the random variables are increased artificially by the factor f -1 to scale the samples into the low-probability domain. This leads to an estimation of a reliability index β( f ) related to scaled samples, known as the scaled reliability index. Asymptotically, there is a linear relation between scaled and un-scaled reliability indices as β( f ) = f . β(1) (Bucher 2009). Therefore, this relationship enables the calculation of β(1) by estimating β( f ) using extrapolation techniques and curve fitting. The advantages and limitation of this method can be listed as: -Moderate accuracy; -Simple implementation; -Very low memory requirement; -Provides the CDF for all threshold levels; -Lack of accuracy in estimation of very low probabilities (possibility of bias). Applying the AS method, several values of f ∈]0;1] are chosen as f i , i{1, ..., w}, where w is the number of f factors, cf. (Sichani et al. 2011a). The estimation process is performed for each selected threshold level , where l{1, ..., w} and k is the number of considered threshold levels to estimate the cumulative distribution function (CDF) of samples (here the stiffness of the monopile).
The utilized methodology in this study has been presented in Figure 6. The process is performed for each  Andersen et al. (2012) factor of f i and a set of (f i , β l (f i )) is generated. Curve fitting is conveyed on the total sets of (f, β l ( f ) ) to extrapolate β l (1) which is respective to f =1 for each threshold level .

Curve fitting
The curve for extrapolation and estimation of β l (1) is proposed by Bucher (2009): Here, coefficients A and B are determined by regression analysis. This curve can satisfy asymptotically linear behaviour as mentioned before. In a general case, Eqn (11) can be rewritten as (Sichani et al. 2011b): where the power C can be estimated through an optimization process (Sichani et al. 2011b). This curve is proposed for more complicated cases, e.g. a wind turbine where the linear curve in Eqn (11) seems to have some errors in estimating the probability (Sichani et al. 2011b).
In this study, based on distribution-curve-behaviour condition (DCBC), Eqn (11) or Eqn (12) is proposed as an improvement of the AS technique.
As expected, the probability of event should be reduced by decreasing the threshold levels at the left tail of the CDFs. If it comes out vice versa, there is an error in the extrapolation process. This fitting is performed basically by Eqn (12) and there is a constraint for switching to Eqn (11) when the probability is increased by decreasing the thresholds. This procedure improves the estimation with more accuracy at high probabilities by Eqn (12), which is strengthened by three parameters, and at low probabilities by Eqn (11), where there are fewer points to fit and a rough fitting with only two parameters provides better accuracy.

Random variables
The random field presented in Section 2 has been used in the simulations. Both standard Gaussian (normal) and quasi random (Sobol) variables (Bucher 2009) have been considered within the process of the random field generation and applied to simulate the samples (i.e. the foundation stiffness, cf. Section 3). This procedure for AS is presented in Figure 6. A comparison between the two mothods of sampling is discussed in the following.

Results and discussion
Simulations have been carried out on the constructed finite element model and based on the flowchart shown in Figure 6 in order to estimate the probability distribu- Fig. 6. Asymptotic sampling process tion of the monopile stiffness. For comparison, a standard Monte Carlo (SMC) simulation has been performed with 10 6 simulations. This amount of simulations can estimate a probability of 10 -5 having a COV = 3.2. This corresponds to a reliability level for serious safety class at failure type II (NKB 1978). The AS methods based on normal and Sobol sampling are compared and verified with SMC simulation as a benchmark. Figure 7 shows the CDF of the horizontal stiffness of the foundation with 10 4 simulations and a curve fitting based on Eqn (12) for the AS method. It seems that the curves estimated by AS match the results of the SMC approach for probabilities higher than 7 × 10 -4 , whereas the curves have some errors in estimating lower probabilities. These errors might be due to the lack of samples in low probabilities so that curve fitting is performed by less support points and three parameters A, B and C. Figure 8 illustrates a curve fitting based on Eqn (11) for AS in the same manner as Figure 7. As shown in this figure, there are some errors in high probabilities (more than 7 × 10 -3 ) compared with Figure 7. However, the errors in Figure 8 are less than the ones presented in Figure 7 for lower probabilities. This might be due to curve fitting is carried out by only two coefficients: A and B. This is cruder than Eqn (12) with three coefficients, which leads to a poor estimation in high probabilities where there are more samples. But stronger estimates are observed in lower probabilities where few samples occur. Figure 9 presents an improved curve fitting by Eqns (11) or (12) plus a constraint (DCBC). It seems that the combined form covers more accurate extrapolation compared with Eqn (11) or Eqn (12). Figure 10 shows the errors in the thresholds at the same probabilities for the Sobol sampling compared with SMC and for curve fitting performed in three ways. As shown in this figure, the DCBC-based curve has smaller errors compared with results from Eqns (11) and (12). The DCBC-based curve Fig. 7. Probability distribution of ky using Eqn (12) by AS and SMC Fig. 8. Probability distribution of ky using Eqn (11) by AS and SMC Fig. 9. Probability distribution of ky using Eqn (11) or (12) by AS and SMC Fig. 10. Error in thresholds value at the same probability using Eqns (11), (12) and DCBC-based curve coincides with Eqns (12) and (11) for high and low probabilities, respectively. It is worth to mention that the non-concurrence part of the dashed curve in Figure 10 belongs to the transition region. The error in this region should be smaller than or equal to the maximum error. Figures 11 and 12 illustrate the CDFs of the rotational and horizontal stiffness of the foundation, respectively, with ten times 10 4 simulations utilizing DCBC-based curve fitting. The curves show the average of the ten simulation sequences. For a comparison between normal and Sobol sampling, Figure 13 shows the coefficient of variation (COV) of estimated probabilities for rotational stiffnesses. Likewise, these COVs are calculated by using ten sequences of 10 4 simulations. As shown in the figures, the accuracy of the results obtained by Sobol sampling is much better than the accuracy achieved with normal sampling, and Sobol sampling is more efficient due to a lower coefficient of variation.

Conclusions
The probability distribution of the stiffnesses for an offshore wind turbine monopile foundation has been determined. Nonlinear p-y curves have been used in combination with the finite element method to estimate the mechanical response of a laterally loaded pile. Asymptotic sampling (AS) and standard Monte Carlo (SMC) methods have been applied to estimate the cumulative distribution functions for the resulting horizontal and rotational foundation stiffnesses measured at the pile cap. The results show the efficiency as well as the accuracy of the AS method versus SMC simulation. The number of simulations is crucial, especially for highly nonlinear problems where the simulation process can be very time consuming. It has been found that an improved AS method combining curve fitting and quasi random (Sobol) variables is more efficient than the SMC technique for analysis of such problems. Fig. 11. Probability distribution of ky using Eqns (11) and (12) and 10 times simulations by AS and SMC Fig. 12. Probability distribution of kr using Eqns (11) and (12) and 10 times simulations by AS and SMC Fig. 13. COVs of estimated probabilities for ky. A comparison between normal and Sobol sampling, R is referred to R-squared term in curve fitting