Un- and reloading stiffness of monopile foundations in sand

Applied Ocean Research 84 (2019) 62–73

Contents lists available at ScienceDirect

Applied Ocean Research journal homepage: www.elsevier.com/locate/apor

Un- and reloading stiffness of monopile foundations in sand Martin Achmus a b

a,b,⁎

, Klaus Thieken

a,b

a

a

, Jann-Eike Saathoff , Mauricio Terceros , Johannes Albiker

T a

Institute for Geotechnical Engineering, Leibniz University Hannover, Germany ACP Grundbauplanung GmbH, Hannover, Germany

A R T I C LE I N FO

A B S T R A C T

Keywords: Monopile Stiffness Eigenfrequency Un- and reloading Offshore wind

The eigenfrequency of offshore wind turbine structures is a crucial design parameter, since it determines the dynamic behavior of the structure and with that the fatigue loads for the structural design. For offshore wind turbines founded on monopiles, the rotational stiffness of the monopile-soil system for un- and reloading states strongly affects the eigenfrequency. A numerical model for the calculation of the monopile’s behavior under unand reloading is established and validated by back-calculation of model and field tests. With this model, a parametric study is conducted in which pile geometry, soil parameters and load conditions are varied. It is shown that of course the rotational stiffness varies with mean load and magnitude of the considered un- and reloading span, but that for most relevant load situations the initial rotational stiffness of the monopile system, i.e. the initial slope of the moment-rotation curve for monotonic loading, gives a good estimate of the actual stiffness. Comparisons of different p–y approaches show that the ordinary API approach considerably underestimates the initial stiffness, whereas the recently developed ‘Thieken’ approach and also the ‘Kallehave’ approach give a much better prediction and thus might be used in the design of monopiles in sand.

1. Introduction In European sea regions, a great number of offshore wind farms have already been erected in recent years and, according to the plans of European governments, an even greater number shall be erected in the next years. In predominately sandy soils as present for instance in the German North Sea regions, the monopile foundation has proven to be an economic solution for the foundation of wind energy converters in water depths up to 30 m or even 40 m. A monopile is a single steel pipe pile of large diameter that is driven into the subsoil. For large water depths and the current generation of wind turbines with rated energy output greater than 8 MW, even in dense sands pile diameters of 7 m–8 m and embedded pile lengths of 25 m–40 m are necessary to fulfill the design requirements. A crucial requirement for offshore wind turbine foundations regards the stiffness under operational lateral loading due to wind and wave action, which strongly affects the eigenfrequency of the overall structure. For the usual “soft-stiff” design it must be ensured that the eigenfrequency lies sufficiently above the 1 P excitation frequency and below the 3 P excitation frequency resulting from the rotational frequency of the wind turbine in order to avoid resonance effects and associated high fatigue loads. During operation, a monopile is mainly loaded by a horizontal force

⁎

H, which acts at a certain distance (eccentricity h) above the embedment point and thus induces a combined H–M (horizontal load H and bending moment M) loading with respect to the embedment point. The effect of the vertical load V on the behavior under horizontal loading is small and can normally be neglected. The bearing behavior of the monopile for a certain eccentricity can be represented either by its loaddisplacement (H-y) or by its moment-rotation (M-θ) curve, wherein y is the pile deflection and θ is the pile rotation at the embedment point (Fig. 1). In Fig. 1 also the secant stiffnesses Ks,y and Ks,θ, which apply for a certain load (H, M), and the initial stiffnesses Ky,0 and Kθ,0, which describe the slope of the load-deformation curve at yHead = 0 or θHead = 0, respectively, are introduced. Due to the rather low frequency of wind and wave loadings, dynamic effects are usually neglected in the determination of stiffnesses. It is common practice to calculate the bearing behavior with a model in which the monopile is considered as a beam supported by springs with non-linear load-displacement relationships representing the soil around the pile. This approach is called p–y method and is recommended by the current offshore guidelines of API [1] and DNVGL [2]. However, p–y curves only describe the monopile behavior under monotonic loading, whereby actually wind and wave loads are cyclic loads and induce unand reloading of the monopile. Fig. 2 elucidates the difference between the secant stiffnesses derived from either “monotonic” or “cyclic” p–y

Corresponding author at: Institute for Geotechnical Engineering, Leibniz University Hannover, Germany. E-mail address: [email protected] (M. Achmus).

https://doi.org/10.1016/j.apor.2019.01.001 Received 22 August 2018; Received in revised form 23 November 2018; Accepted 2 January 2019 0141-1187/ © 2019 Elsevier Ltd. All rights reserved.

Applied Ocean Research 84 (2019) 62–73

M. Achmus et al.

Fig. 1. Monopile system, p–y method and resulting load-displacement curves.

Fig. 3. Measured and calculated eigenfrequencies of offshore wind turbines founded on monopiles (Damgaard et al. [4]). Fig. 2. Comparison of secant stiffnesses and un- and reloading stiffnesses.

foundation is considered rigid, and once the initial stiffness of the moment-rotation curve calculated with the p–y approach of the current offshore guidelines was considered. The measured frequencies all lied between the two predictions and thus indicated that even the initial stiffness calculated with the ordinary p–y approach underestimates the actual stiffness of the monopile foundation. An accurate prediction of monopile stiffnesses is important not only for new wind farm projects, but also for re-assessments of existing foundations in the scope of considerations regarding a possible lifetime extension of the structures. Such optimizations are based on the frequency response measured during turbine operation (as done by Damgaard et al. [4]) which are then used as validation basis for overall OWT models considering more realistic p–y relationships. The following numerical investigation regarding the un- and reloading stiffness of monopiles in sand soil was initiated and validated by the authors in the scope of such a project study for a wind farm in the North Sea. The numerical results presented hereinafter are used to provide recommendations regarding the choice of p–y approaches which can be used to realistically predict stiffnesses and with that eigenfrequencies of offshore wind turbines founded on monopiles in sand.

curves and the un- and reloading stiffness. Evidently, the use of the monotonic and even more the use of the cyclic secant stiffness would considerably underestimate the stiffness. In practical applications, usually the monopile’s load-displacement characteristic for typical operational loads is calculated and from the obtained results a stiffness matrix for the embedment point is derived, which substitutes the pile-soil system in a dynamic calculation of the wind turbine structure. By using such a stiffness matrix, the interaction of horizontal and moment loading and its effect on the stiffness is taken into account. The entries in the stiffness matrix are usually determined either from the secant stiffnesses or from the initial stiffnesses of the load-deformation curves. It is questionable – and that is the subject of this paper – whether these stiffnesses represent the relevant un- and reloading stiffnesses. Kallehave et al. [3] reported for offshore wind turbines founded on monopiles in sandy soils at the Walney offshore wind farm that the prediction of the eigenfrequencies by use of the initial stiffness determined with the p–y approach given in the offshore guidelines revealed an under-prediction of around 5–7%. Therefore, this p–y approach obviously yields much too small monopile stiffnesses. This finding was confirmed by similar eigenfrequency measurements reported by Damgaard et al. [4]. Here the foundations consisted of monopiles located in layered soil with both cohesionless and cohesive soils. Fig. 3 compares predicted and measured eigenfrequencies. Regarding the calculated foundation stiffness, once rigid fixation of the tower at the sea bottom was assumed, which means that the monopile

2. State-of-the-art of the calculation of monopile behavior in sand soils As mentioned above, the behavior of monopiles is usually calculated by p–y methods (see Fig. 1). The first, nonlinear p–y formulation was proposed by McClelland & Focht [5]. Since then, several modified 63

Applied Ocean Research 84 (2019) 62–73

M. Achmus et al. m

formulations have been developed, the most established one by Reese et al. [6], who proposed sectional p–y curves calibrated on the “Mustang Island” field test (Cox et al. [7]). This formulation was incorporated in the regulations of the American Petroleum Institute until it was replaced by an approach utilizing a hyperbolic tangent function proposed by Murchison & O’Neill [8], which was calibrated on further field tests of different kinds of pile geometries and materials. According to the current offshore guidelines ([1,2,9]), the maximum mobilized soil reaction force per unit length of the pile pu depends on the regarded depth under seabed z, the submerged unit weight of the soil γ’, the pile diameter D and on the angle of internal friction ϕ’ of the sand:

pus = (c1 z + c2 D) γ ′ z

(1a)

pud = c3 D γ ′ z

(1b)

kKallehave =

(2)

with A = 3.0 − 0.8 z / D ≥ 0.9 for static (monotonic) loading. Here p is the soil resistance per unit length of the pile and y is the horizontal pile deflection. The parameter kAPI also given in [1] determines the initial slope of the p–y curve Epy (cf. Fig. 1):

Epy = kAPI z

⎜

⎟

⎜

0.5

⎟

(4)

Here, the recommended initial stiffness coefficient kKallehave depends on the pile diameter and on the initial bedding stiffness according to the current p–y method (kAPI∙z) in a reference depth z0 = 2.5 m. The reference diameter is equal to the diameter of the Reese et al. [6] tests D0 = 0.61 m. The dimensionless parameter m, which rules the course of the initial stiffness with the depth, is suggested to be set to m = 0.6. Applying this approach to the monopiles at the Walney offshore wind farm, Kallehave et al. [3] showed a better agreement of predicted and measured eigenfrequencies, although the measured values were still 2–3% higher. Thieken et al. [16] developed a new p–y method for piles in sand by means of parametric studies with a sophisticated numerical model. The basis of the approach is a p–y curve formulation valid for a pile of infinite length which exhibits a constant horizontal deflection. An iterative procedure is used to account for the effect of the actual pile deformation mode on the p–y curves. Herein, also the effect of the pile tip, i.e. the distance of the considered point to the pile tip, on the p–y curves is accounted for. The input parameters of this approach are the soil’s unit weight γ’, the friction angle ϕ’, the dynamic (small-strain) oedometric soil stiffness Esd or Eoed,0 (respectively the dynamic shear modulus G0 and the Poisson’s ratio ν) and the oedometric soil stiffness ES (or Eoed) obtained in standard oedometric compression tests. It has been thoroughly evaluated and is proven to give better predictions for the load-bearing behavior (in terms of deflection line, bending moment, local bedding resistance) of monopiles than other approaches (Achmus et al. [17]). Recently, another modified calculation approach for large-diameter piles in sand has been developed in the scope of the PISA project (Byrne et al. [18]). Herein, besides p–y springs also rotational springs and a pile tip spring are introduced to the beam-spring model. The spring characteristics are to be determined by calibration with a numerical model. This method is not further considered here, since details of this approach are not yet published.

The first mentioned equation applies to small depths (pus) and the second equation to greater depths (pud), the smaller of both values is to be considered. The influence of the relative density of the sand is described by the factors c1, c2 and c3, which are given in [1] dependent on the angle of internal friction. The p–y curve is described by the following equation:

k z p = A pu tanh ⎛⎜ API y ⎞⎟ ⎝ A pu ⎠

1 z D ⋅kAPI ⋅z 0⋅⎛ ⎞ ⋅⎛ ⎞ z ⎝ z 0 ⎠ ⎝ D0 ⎠

(3)

Therefore, the approach assumes a linear increase of the initial spring stiffnesses with depth. kAPI as given in [1] is dependent on the relative density ID and with that also on the angle of internal friction, which means that a direct correlation of stiffness and shear strength is assumed. For instance, for a relative density of ID = 60% (corresponding according to [1] to ϕ’ = 36°) kAPI = 25.5 MN/m3 applies, whereas for ID = 80% (corresponding to ϕ’ = 40°) kAPI = 41.9 MN/m3 can be read off the diagram given in [1]. Eqs. (1a), (1b) and (2) are mainly based on investigations on flexible pile-soil systems. For instance, Reese et al. [6] tested a 21 m-long steel tube pile having a diameter of 61 cm. Based on experience in the oil and gas industry, the p–y method proved to be sufficiently accurate for flexible piles with diameters up to two or even three meters. However, for larger pile diameters, several numerical investigations showed that the horizontal deflections are underestimated for extreme loads [10–12]. In contrast, as already mentioned in Section 1, experience from operating offshore wind farms with monopile foundations indicate that the deflections for (small) operational loads are overestimated and thus the foundation stiffness is underestimated [13,3,4]. In this regard, it has to be mentioned that it is still under discussion whether or not this larger stiffness partly results from dynamic effects such as inertia effects and damping. However, as mentioned above, it is current practice to neglect the dynamic effects. Likewise, also densification effects on the surrounding soil by small operational loads are discussed to be the reason for the larger foundation stiffness in turbine operation (LeBlanc et al. [14]). At least for the very dense and strongly overconsolidated sands in the North Sea, this effect should, however, play a subordinate role. Several researchers developed approaches to adapt the original p–y method in order to make it suitable also for large-diameter piles, see for details Thieken et al. [15]. The approach of Kallehave et al. [3] specifically addresses the observed stiffer behavior of monopile systems in sand under small operational loads. They proposed a modification of the initial stiffness parameter in the original p–y method according to Eq. (4):

3. Numerical model 3.1. General The numerical simulations were carried out with the finite element code Plaxis3D (Brinkgreve et al. [19]). Due to the symmetry conditions of the loading and the pile only one half of the system could be modeled. By means of preliminary analyses a sufficient mesh fineness with regard to solution accuracy was found and the model dimensions were chosen such that there is no impact of the boundaries. The final model is discretized with roughly 200,000 elements and has a width of 25 D in loading direction, a depth of 2 L and a length of 6 D perpendicular to the loading direction (Fig. 4). The pile was modelled as an open ended tubular pipe pile. For the steel material, linear-elastic behavior was assumed with a Young’s modulus E = 2.1E8 kN/m², a Poisson’s ratio ν = 0.27 and an effective steel unit weight γs = 68 kN/m³. Above the embedment point a lever arm was modelled in order to apply the horizontal load with a load eccentricity and thus introduce a moment load. An elasto-plastic contact interface was considered between the inside and the outside of the steel pile and the adjacent soil. The maximum shear stress in the contact surface τmax results from the product of the horizontal stress σH and the tangent of the contact friction angle δ = 2/3ϕ’, with ϕ’ being the internal friction angle of the soil. Thereby, slip between soil and steel, which might occur under greater loads, is enabled. In the numerical simulation, firstly the initial stress state in the soil was calculated with a k0-procedure, i.e. σH = k0 ∙ σV, with k0 set to k0 = 1 – sin ϕ’. Then the monopile was installed in a “wished-in-place” 64

Applied Ocean Research 84 (2019) 62–73

M. Achmus et al.

Fig. 5. Relation between shear modulus and shear strain.

parameters. Additionally, for un- and reloading a Poisson’s ratio νur has to be defined. A shortcoming of the material law as implemented in Plaxis3D is that the stress-dependency of the small strain modulus G0 and of the “normal strain” modules is ruled by the same exponent m, i.e. mG = m. However, for non-cohesive soils as considered here, this is a justifiable approach, since both exponents normally lie in a similar range. For sands and gravels, Benz [22] reports mG-values of 0.41 to 0.57. Typical m-values of sands lie between 0.5 for a very dense and 0.65 for a medium dense sand. If the shear strains are smaller than a threshold value, the soil stiffness is corrected in order to account for the small strain effect. This is governed by the dependence of the shear modulus G on the magnitude of shear strain γ. Here, the approach of Dos Santos and Correia [24] is used (Fig. 5):

Fig. 4. Finite element mesh with system denominations.

procedure, i.e. the shell elements representing the pile were activated. Doing that, effects of the installation procedure on the stress state of the soil are disregarded, which seems acceptable because such effects have only minor impact on the pile behavior under lateral loading (see also Li et al. [20]). As mentioned above, also the effect of vertical loads on the lateral bearing behavior can normally be neglected (see e.g. Achmus and Thieken [21]). However, to ensure realistic conditions, a vertical compressive load of 8 MN was then applied to the pile. In the last stage, the lateral and moment loading was applied with alternating virgin loading steps and un- and reloading cycles. 3.2. Constitutive law for the soil

G 1 = γ G0 1 + 0.385 γ

A special aspect of soil behavior is that the stiffness response to a certain load is dependent on the magnitude of shear strain induced. For the problem considered here this is a crucial feature, because the initial stiffness and the un- and reloading stiffness under small cyclic load amplitudes is significantly affected by the small strain stiffness. Here the Hardening Soil small strain stiffness (HSsmall) material law according to Benz [22] was applied. This model is an upgraded version of the sophisticated hardening soil model (HSM) according to Schanz [23], which is an elasto-plastic model with isotropic hardening and also enables the consideration of stress-dependent soil stiffness, for instance. The HSsmall model additionally takes the strain-dependency of soil stiffness into account. The stress-dependency of stiffness is ruled by the following equations, in which either the current maximum principal stress σ1 or the minimum principal stress σ3 is used as a stress measure ([19],[22]):

ref

The reference shear strain γref corresponds to the shear strain at which the soil stiffness is degraded to 72.2% with respect to the maximum value G0. In the HSsmall model, the reduction of G with shear strain is limited to the value Gmin = Eur/(2 (1+νur)) belonging to “normal” shear strain magnitudes. The shear strength of the soil is described by a Mohr-Coulomb failure criterion with the shear parameters ϕ’ and c’. Plastic strains are accounted for by isotropic hardening with a yield surface consisting of a cone and a cap in the principal stress space. For the cone, a non-associated flow rule is used by considering a dilatancy angle ψ. Further details on the HSsmall material law can be found in Benz [22]. 4. Validation of the numerical model

m

⎛σ ⎞ Eoed = Eoed, ref ⎜ 1 ⎟ p ⎝ ref ⎠

Benz [22] proved that the HSsmall material law yields realistic predictions even for problems with complex boundary conditions. For a horizontally loaded pile, Thieken et al. [16] showed very good agreement of numerical predictions with the results of the “Mustang Island” field tests of Cox et al. [7]. However, this model validation only regarded the pile behavior under monotonic loading. In the literature, only few tests of piles with small aspect ratios L/D – as typical for monopiles – with well-documented un- and reloading behavior can be found. Recently, Li et al. [20] reported suitable tests with piles of a diameter of D = 0.34 m, and in 2016 tests with monopiles D = 4.3 m were carried out at a test site close to Cuxhaven in Germany. The results of these tests were used for the validation of the aforementioned numerical model with respect to un- and reloading behavior. In both cases, first the model parameters decisive for the behavior under monotonic loading were adapted to the corresponding load-displacement curve. In particular, the stiffness parameters Eoed,ref, E50,ref and G0,ref were varied to optimally reflect the load-displacement

(5a)

m

⎛σ ⎞ E50 = E50, ref ⎜ 3 ⎟ p ⎝ ref ⎠

(5b) m

⎛σ ⎞ Eur = Eur , ref ⎜ 3 ⎟ p ⎝ ref ⎠

(5c)

m

⎛σ ⎞ G0 = G0, ref ⎜ 3 ⎟ p ⎝ ref ⎠

(6)

(5d)

Here, Eoed is the stiffness under oedometric loading, E50 is a stiffness value which determines the behavior under triaxial compression states, Eur is the un- and reloading stiffness and G0 is the small strain shear modulus. The reference stress value is pref = 100 kN/m2. The reference values of the four stiffnesses and the exponent m are material 65

Applied Ocean Research 84 (2019) 62–73

M. Achmus et al.

curves measured, whereas other parameters were chosen as typical values for sand (m = 0.5, νur = 0.25, γ0.7 = 10−4). The angle of friction ϕ’ was chosen according to the relative density reported in the references, and the angle of dilatancy was estimated by the common expression ψ = ϕ’ – 30°. Afterwards, using these parameters, the un- and reloading cycles were back-calculated. Herein, the following ratio of stiffnesses for un- and reloading Eur and the stiffness E50 was used, which is typical for sand soils (Brinkgreve et al. [19]):

Eur =4 E50

(7)

It should be noted that not only Eur, but also the other material parameters and in particular the parameter G0,ref affect the behavior under un- and reloading. Fig. 7. Load-deflection curve from a field test with a monopile D = 4.3 m and from numerical back-calculation.

4.1. Pile tests of Li et al. [20] Li et al. [20] conducted field tests on steel pipe piles in almost homogeneous dense to very dense sand. The test pile considered here had an embedded length of L = 2.2 m and a diameter of D = 0.34 m, implying an aspect ratio of about L/D = 6.5, which is quite slender but not totally out of range for a monopile. The pile wall thickness was t = 14 mm. The horizontal load was applied with an eccentricity of h = 0.4 m. From different load levels reached, partial unloading with subsequent reloading was carried out, as depicted in Fig. 6. The soil parameters could be deducted from the results of a CPT test given. To get good agreement with the monotonic load-displacement curve, the parameters of the numerical model were adjusted within plausible ranges. Fig. 6 shows that a reasonable agreement of calculation and measurement could be achieved. From a load of around 20 kN, unloading to 9 kN was done, and from a load of 28 kN unloading to 17 kN was simulated. This means that the cyclic load spans were 55% and 40% of the previously reached maximum load, respectively. The comparison with the measured results shows that a much greater hysteresis was observed in the experiments, which partly might have been caused by some creep deflections occurring in the tests. However, the secant stiffnesses of the un- and reloading loops, which are decisive here, are quite similar in simulation and experiments.

the load-displacement curves can only be shown here with displacements in dimensionless form. The soil mainly consisted of fine to medium sands, which were in a medium dense to dense state in the upper half of the pile and in a dense to very dense state underneath. Several CPT tests at the test site were available from which the material parameters of the HSsmall material law could be derived. Using a reasonable parameter set, a good agreement between the experimentally determined load-deflection curve for virgin loading and the curve determined numerically could be obtained. The application point of the horizontal load was about 1 m above the embedment point. The maximum load applied was around 18 MN. The load-displacement curves from the field test and from the numerical simulation are depicted in Fig. 7. From two different load levels, full unloading and subsequent reloading of the pile was carried out. It is observed that up to an unloading stage of about 50% of the applied load, the curves of field test and numerical back-calculation match quite well. For further unloading, the numerical back-calculation underpredicts the reduction of pile head displacement and thus the pile relaxation. The deviations in this region might be caused by the fact that the chosen material law for the soil comprises only isotropic hardening. In reality, at greater unloading amplitudes the soil stiffness is degraded relatively strong due to plastic deformations, which necessitates a kinematic hardening law. Concluding, it can be stated that the numerical model is capable to represent the behavior of monopiles under un- and reloading with reasonable accuracy for cases of unloading up to at least about 50% of the formerly applied load. This is a typical loading situation of foundations of offshore wind turbines, where wind loads are oscillating around the wind load belonging to the mean wind speed. Therefore, it is concluded that the numerical model can be applied for further investigation of the un- and reloading stiffness of horizontally loaded pilesoil systems.

4.2. Pile tests in Cuxhaven/ Germany Gattermann et al. [25] reported about large-scale field tests on monopiles with diameters of D = 4.3 m and an embedded length of around 18.2 m (L/D ≈ 4.2), which were conducted in the scope of an industry research project. The results of these tests are confidential, but the authors had access to the data. To keep the required confidentiality,

5. Results for a reference system 5.1. System configuration and simulation procedure Investigations regarding the reloading stiffness of a pile-soil system were firstly carried out for a reference system. The monopile dimensions (D = 8 m, L = 32 m) and the soil conditions (very dense sand) and soil parameters were chosen as typical (based on the authors’ experience) for an offshore wind farm in the German North Sea with a water depth of around 35–40 m. The horizontal load was applied with an eccentricity of h = 48.0 m, which implies a dimensionless load eccentricity of h/L = 1.5 or h/D = 6. The wall thickness of the steel pipe pile amounts to t = 86.35 mm, which is the minimum value for a driving

Fig. 6. Load-deflection curve from a field test of Li et al [20] and from numerical back-calculation. 66

Applied Ocean Research 84 (2019) 62–73

M. Achmus et al.

be represented by the secant stiffnesses Ks,y and Ks,θ as defined in Fig. 8. Obviously, for the considered case the variation of the secant stiffness with maximum load is quite small. Fig. 9 further elucidates the results for the reference system. Here, the horizontal tangent stiffnesses, i.e. the current slopes of the loaddeflection curves, are depicted dependent on the current deflection. In Fig. 9 left, the results for 50% cyclic load span are depicted. The slope of the envelope curve for virgin loading monotonically decreases from the initial stiffness (Ky,0 ≈ 750 MN/m) to values less than half of the initial value at deflections greater than around 5 mm. The unloading curves start with a slightly higher and end with a slightly smaller stiffness than Ky,0. With increasing load (or displacement) level, the deviations from the initial stiffness increases. Also, for increasing un- and reloading span, the lower bound value decreases. Fig. 9 right shows the envelope curves of un- and reloading stiffnesses for different cyclic load spans. At the end of unloading, the smallest stiffnesses apply, whereas at the beginning of reloading the highest stiffnesses occur. For a small cyclic load span of 25%, the maximum and minimum tangent stiffnesses are only slightly greater or smaller, respectively, than the initial stiffness derived from the virgin loading curve. With increasing cyclic load span and with increasing head displacement, the minimum stiffness becomes smaller. However, for a cyclic load span of 50%, the minimum tangent stiffness is at maximum only 15% smaller than the initial stiffness, and even for a cyclic load span of 100%, the maximum deviation is 35%.

Table 1 Parameters of very dense sand used in the calculations. Buoyant unit weight γ’ Reference oedometric stiffness Eoed,ref Reference stiffness E50,ref Reference un- and reloading stiffness Eur,ref Reference small strain shear modulus G0,ref Stress exponent m Poisson’s ratio for un- and reloading νur Shear strain γ0.7 Angle of internal friction ϕ’ Angle of dilatancy ψ

[kN/m3] [MN/m2] [MN/m2] [MN/m2] [MN/m2] [-] [-] [-] [°] [°]

11.0 120 120 480 268 0.5 0.25 0.0001 40 10

shoe with respect to pile buckling during pile driving recommended in API [9]. A change of the wall thickness over depth, which is usually executed in practice, was neglected. The pile was embedded in homogeneous very dense sand soil. The soil parameters are collected in Table 1. The small strain shear modulus G0 was chosen such that the ratio between the small strain oedometric stiffness Eoed,0 and the oedometric stiffness Eoed amounts to 4, which is a typical value for very dense sands. Assuming elastic behavior on the very small strain level, G0 can be calculated from Eoed,0 by the following equation stemming from elasticity theory:

G0 = Eoed,0

2 1 − νur − 2νur 2 2(1 − νur )

(8)

For νur = 0.25, this equation yields G0 = Eoed,0/3. In the numerical model, G0 and thus also Eoed,0 varies with σ3 and Eoed with σ1 (see Eqs. (5a) and (5d)), and hence the ratio Eoed,0/Eoed is not a constant. G0,ref was chosen here such that for the initial stress state with σ3 = k0 σ1 the desired ratio Eoed,0/Eoed = 4 applied. Firstly, for the purpose of comparison, a simulation with purely monotonic load increase (virgin loading) was carried out. Subsequently, simulations with load increase in steps followed by totally 10 intermediate un- and reloading cycles were carried out. Herein, the load was increased in constant steps after each cycle. The load steps were chosen such that the maximum load applied in the last cycle led to a pile head rotation of around 0.1°. Based on the authors’ personal experience with monopile designs, this value is believed to be a typical maximum value for operational conditions, for instance under maximum operational wind speed. For the reference system, five simulation series were carried out with un-/ reloading load spans (in the following also termed cyclic load span) of 25%, 50%, 75% and 100% (full unloading) of the previously applied maximum load. In the evaluation of the results, it should be kept in mind that the validity of the numerical model is at least questionable for large cyclic load spans of 75% and 100% (see also Section 4.2). The results are evaluated in terms of load-deflection and momentrotation curves. From these curves, the horizontal and rotational stiffnesses of the monopile-soil system can be derived.

5.3. Secant stiffnesses for un- and reloading The definition of the secant stiffnesses Ks,y and Ks,θ for the un- and reloading cycles is shown in Fig. 8. As already stated, from the graphs shown in Fig. 8 differences of these stiffnesses for different maximum load levels can hardly be seen. For a more profound analysis, the determined secant stiffnesses are plotted in Fig. 10 against the mean values of lateral load H and moment M for the considered un- and reloading loops, i.e. Hmean = Hmax – 0.5 ΔH and Mmean = Mmax – 0.5 ΔM. In these diagrams, furthermore, the secant stiffnesses derived from the virgin load-deformation curves are included. As to be expected for a system where partly plastic deformations occur, the numerically derived secant stiffnesses for un- and reloading are significantly larger than the secant stiffnesses for virgin loading, with the difference becoming greater with increasing load level. Interestingly, considering even cyclic load spans of up to 100%, the variation of un- and reloading stiffnesses with the magnitude of mean load or load span is quite small and the stiffness values lie in the range of the initial stiffnesses of the virgin load-displacement curves. This applies particularly for the most important rotational stiffness. This means that for the monopile-soil system considered the un- and reloading stiffness and thus the eigenfrequency of the overall system is almost independent of the current mean and cyclic load level as long as only one-way cyclic loading occurs. The results also indicate that the relevant system stiffnesses can be realistically approximated by the initial stiffnesses of the virgin load-deflection and moment-rotation curves.

5.2. Load-displacement curves Fig. 8 shows the load-deflection and the moment-rotation curves for the monotonically and the cyclically loaded reference system, the latter exemplarily for a cyclic load span of 50% of the previously applied load. First of all, it is observed that the cyclic envelope curve deviates only slightly from the virgin (monotonic) loading curve. Obviously, the effect of intermediate application of un- and reloading cycles on the system behavior under subsequent load increase is very small. The un- and reloading hysteresis, i.e. the area enclosed by the unand reloading load-displacement curves, is rather small, which indicates only small system damping for the considered loading states. However, it should be noted that the model used here was not validated with respect to realistic consideration of damping, but focuses on stiffness of the monopile-soil system. The un- and reloading stiffness can

6. Parametric study 6.1. General procedure and parameter variations In the scope of the parametric study, it shall be checked whether the finding for the reference system, that the un- and reloading stiffness can be approximated by the initial stiffness of the virgin loading curve, can be generalized. For that purpose, geometry, loading and soil conditions were varied:

• The geometry of the monopile is described by the pile diameter D,

the embedded pile length L and the wall thickness t. Firstly, the

67

Applied Ocean Research 84 (2019) 62–73

M. Achmus et al.

Fig. 8. Numerically obtained load-deflection and moment-rotation curves for monotonic and cyclic loading with 50% cyclic load span of the reference system.

Fig. 9. Tangent stiffnesses for the reference system for 50% cyclic load span (left) and bandwidth of tangent stiffnesses for different cyclic load spans (right).

• •

oedometric stiffness were considered by varying the parameter ratio Eoed,0/Eoed. Finally, once the internal angle of friction was varied between 38° and 42°, keeping the other parameters as in the reference case, and once the relative density of the sand was varied by considering combinations of friction angles and stiffness parameters typical for loose to medium dense, medium dense and dense sand. Regarding the effect of relative density on the parameters of the HSsmall material law, reference is made here to Brinkgreve et al. [26].

embedded length was varied between 24 m and 44 m, thereby keeping the other parameters constant as used in the reference system. This means that the aspect ratio L/D and therewith the effect of the flexibility of the pile-soil system is investigated. Secondly, the wall thickness t is varied, which also changes the system flexibility, but for a given L/D-ratio. Finally, the pile diameter was varied between 4 m and 10 m. Herein, L/D = 4 was applied as in the reference system, i.e. the pile length was adapted to the chosen diameter. In the reference system, the eccentricity h of the horizontal load with respect to the embedment point was set to h = 48 m, i.e. h/ L = 1.5 and h/D = 6. To investigate the effect of the eccentricity, h was varied between 32 and 112 m, thereby keeping L and D constant as in the reference system. Regarding soil conditions, stiffness parameters and the internal angle of friction of the sand were varied. Firstly, reference oedometric stiffnesses Eoed,ref between 80 MPa and 160 MPa were applied, and secondly, different ratios of small strain stiffness to the

Table 2 shows all system configurations which were considered in the parametric study. For each of these configurations, the same calculations as for the reference system were carried out. Both horizontal and rotational secant stiffnesses were determined dependent on maximum load level and cyclic spans of loads and/or deformations. For the evaluation of the results, the rotational secant stiffness Ks,θ is used in the following, since this value affects the eigenfrequency of the wind turbine system more than the horizontal secant stiffness. Ks,θ-

68

Applied Ocean Research 84 (2019) 62–73

M. Achmus et al.

Fig. 10. Secant stiffnesses against the cyclic load spans for different relative levels of unloading.

values for different mean rotations θmean in the un- and reloading loops were determined. Fig. 11 shows the results exemplarily for the reference system. From these graphs, secant stiffnesses for θmean-values of 0.005°, 0.01°, 0.015°, 0.02° and 0.05° and cyclic load spans of 25%, 50%, 75% and 100% were read off. With these values, the ratios of the secant stiffness and the initial stiffness of the virgin loading curve Ks,θ/ Kθ,0 were determined. These values indicate the effect of varying loading conditions on the rotational stiffness and can be compared for systems with different pile geometries and soil conditions.

6.2. Results for varied pile geometries and eccentricities of the load Fig. 12 shows the results in the aforementioned form (Ks,θ/Kθ,0) for the systems with variable pile geometries and load eccentricities. The graphs show quite similar results for all considered systems. The ratio of un- and reloading stiffnesses and initial stiffness is very close to 1 for small cyclic load spans. With increasing mean rotation, i.e. with increasing maximum load, the ratio is even slightly greater than 1. The greater the relative cyclic load span is, the smaller is the stiffness ratio. For cyclic load spans greater than 50%, the stiffness ratio is slightly

Table 2 System configurations considered in the parametric study. Designation

Reference Pile L/D

Pile t

Pile D

Load h

Soil Eoed

Soil Eoed,0/Eoed

Soil ϕ’ Soil rel. density

Pile

Load

Soil

D (m)

L (m)

L/D –

t (mm)

h (m)

Eoed,ref (MN/m2)

Eoed,0/Eoed (MN/m2)

ϕ’ °

m –

8 8 8 8 8 8 8 8 8 10 6 4 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8

32 24 28 36 40 44 32 32 32 40 24 16 32 32 32 32 32 32 32 32 32 32 32 32 32 32 32 32 32 32

4 3 3.5 4.5 5 5.5 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4

86.35 86.35 86.35 86.35 86.35 86.35 102.35 70.35 54.35 86.35 86.35 86.35 86.35 86.35 86.35 86.35 86.35 86.35 86.35 86.35 86.35 86.35 86.35 86.35 86.35 86.35 86.35 86.35 86.35 86.35

48 48 48 48 48 48 48 48 48 48 48 48 112 96 80 64 32 48 48 48 48 48 48 48 48 48 48 48 48 48

120 120 120 120 120 120 120 120 120 120 120 120 120 120 120 120 120 80 100 140 160 120 120 120 120 120 120 60 80 100

4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 3 3.5 4.5 5 4 4 5 4.5 4.2

40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 38 42 33 35 38

0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.65 0.6 0.55

69

Applied Ocean Research 84 (2019) 62–73

M. Achmus et al.

absolute pile diameter, a more pronounced tendency of smaller stiffness ratios with decreasing diameter can be seen. However, in the range of typical monopile diameters from 6 m to 8 m, the differences are also marginal. Thus, the general findings for the reference system, that the un- and reloading stiffness can be well approximated by the initial stiffness, hold for monopiles in dense sands for all geometry and loading conditions of practical relevance. 6.3. Results for variation of soil characteristics Fig. 13 depicts the results for the reference system in conjunction with systems of varied soil stiffness characteristics. A change of the oedometric reference modulus Eoed,ref shows only for the small cyclic load span of 25% a noticeable effect on the considered stiffness ratio with a slight tendency of increasing ratio with increasing soil stiffness (Fig. 13 left). A change of the small strain stiffness value affects both the initial stiffness of the system and the un- and reloading stiffnesses. As a result, the ratio of these stiffnesses is only marginally affected (Fig. 13 right). Tendencially, the bandwidth of stiffness ratios increases with increasing small strain stiffness. However, the differences are small. Fig. 14 shows the effects of varying shear strengths of the soil. An increase of the internal friction angle leads to an increase of the stiffness ratio (Fig. 14 left). The obvious reason is that less soil plastification due to increased strength leads to higher un- and reloading stiffness, whereas the initial stiffness of the system is unaffected by the shear strength. However, Fig. 14 right (referring to the reference and the last three data sets given in Table 2; for ease of presentation relative densities of 50%–80% were assigned to the data sets) shows that there is almost no effect on the stiffness ratio if a change in shear strength is accompanied by a corresponding change in stiffness, as to be expected if the relative density of a material varies. Therefore, the general findings

Fig. 11. Evaluation of rotational secant stiffnesses for the reference system.

smaller than 1. For spans of 50%, the stiffness ratio is in all cases greater than 98%. Even for spans of 100%, the stiffness ratios are almost always greater than 95%. Only for monopile diameters of 6 m and 4 m and mean rotations of 0.05°, stiffness ratios of 94% and 93%, respectively, apply. The aspect ratio L/D, the pile wall thickness t and the eccentricity of the load affect the stiffness ratios only marginal. With regard to the

Fig. 12. Stiffness ratios for monopiles in dense sand with varying geometries and load eccentricities. 70

Applied Ocean Research 84 (2019) 62–73

M. Achmus et al.

Fig. 13. Stiffness ratios for a monopile D = 8 m, L/D = 4 in sand with varying stiffness characteristics.

validation of the results by field measurements is highly desirable. The results of this study indicate that the un- and reloading stiffness of monopiles in sand, which is decisive regarding the eigenfrequency of the overall system, is in cases of practical relevance almost identical with the initial stiffness of the virgin load-displacement curve or, depending on the loading condition, only slightly smaller. This means that in the design of a monopile a p–y method must be applied which accurately captures the initial stiffness Epy of the p–y curves. The stiffnesses or the stiffness matrix derived from that method can then be applied in the determination of the structure’s eigenfrequency and in dynamic calculations of the overall structure. Fig. 15 shows for the reference system a comparison of the Epy-values determined numerically (FEM), with the original (API) p–y method and with the methods of Kallehave et al. [3] and Thieken et al. [16] presented in Section 2. The Kallehave approach gives diameter-dependent values and yields for the considered large-diameter pile (D = 8 m) much higher Epy-values than the API approach. In the Thieken approach, Epy depends on the small strain stiffness Eoed,0 (or G0) and not on the pile diameter. Moreover, Epy is affected by the pile deformation mode. The stiffness is increased close to the soil surface and in the vicinity of the pile tip, whereas it is decreased in the vicinity of the rotation point. Therefore, as shown in Fig. 15, close to the surface Epy is greater than the Kallehave values, whereas at greater depths significantly smaller values apply. In the vicinity of the pile tip, again greater Epy-values arise from the Thieken approach. The comparison with the FEM results shows that in general the distribution of initial stiffness agrees with the numerical simulation results. The very high values close to the rotation point have no practical effect. Theoretically, Epy = p/y becomes infinite at y = 0. However, since both y and p are close to zero in this region, the actual value of Epy is of minor importance. Close to the pile tip, also an increase of Epy is found. The increase assumed in the Thieken approach is much greater, since here – in contrast to the FEM solution – also the effect of a toe shear force is

for the reference system also hold for monopiles in medium dense sands. 7. Discussion Two main effects govern the ratio of the stiffness for un- and reloading and the initial stiffness of a monopile system. The initial stiffness depends on the small strain stiffness of the soil at the initial stress state. In contrast, the tangent stiffness at the beginning of unloading depends on the small strain stiffness for the stress state reached under the preceding load step. Since the soil stiffness increases with increasing stresses, the system stiffness for very small un- and reloading spans should be greater than the initial system stiffness. That is the reason that for the small cyclic load span of 25% in almost all considered cases stiffness ratios Ks,0/Kθ,0 greater than 1 were determined. For greater load spans, the effect of increased stress level is counteracted by the fact that the secant stiffness of an un- and reloading loop is smaller than the tangent stiffness at the beginning of unloading. The difference is the greater, the greater the un- and reloading span is. Thus, for the greatest considered load span of 100%, the smallest stiffness ratios were obtained. The numerical simulation results confirmed the expectations stated above, but it also showed that the differences between the initial stiffness and the un- and reloading stiffnesses of a monopile system are quite small. For cyclic load spans less than 50%, which are normally to be expected in turbine operation, the un- and reloading stiffness is at maximum 2% greater and 2.5% smaller than the initial stiffness for all systems considered. Even for spans of 100%, at maximum around 5% smaller stiffnesses were determined. However, it must be admitted that the validity of the numerical model for such high cyclic load spans is questionable. Based on the results of the back-calculations of field tests presented in Section 4, it is believed that the above mentioned finding holds for un- and reloading spans up to at least 50%. Of course, further

Fig. 14. Stiffness ratios for a monopile D = 8 m, L/D = 4 in sand with varying internal angles of friction and relative densities. 71

Applied Ocean Research 84 (2019) 62–73

M. Achmus et al.

Fig. 15. Comparison of the initial p–y stiffnesses for monopiles in dense sand.

Fig. 16. Comparison of the initial rotational stiffnesses Kθ,0.

soil cannot be taken into account in this approach.

covered. Fig. 16 compares the initial rotational stiffnesses Kθ,0 calculated with the numerical model and calculated with the API, the Kallehave and the Thieken p–y approach. The soil parameters of the reference system have been applied. In order to get a complete picture, monopile diameters of 6 m and 8 m and load eccentricities of 48 m and 112 m are considered, i.e. with respect to Table 2, additional system configurations were investigated. Evidently, the API approach considerably underestimates the initial system stiffness for all cases. On the contrary, both the Kallehave and the Thieken approach comply much better with the numerical results. Actually, the Thieken approach predicts a stiffer system behavior than the Kallehave approach and therewith yields even better agreement with the numerical results than the Kallehave approach. The obvious reason for this is that the Thieken approach directly accounts for the actual small strain soil stiffness also applied in the numerical simulation, whereas the Kallehave approach uses only the internal friction angle as an input, thereby assuming a certain correlation between strength and stiffness of the soil. Therefore, the actual ratio of small strain stiffness to oedometric stiffness for a certain

8. Conclusions The eigenfrequency of an offshore wind turbine tower founded on a monopile depends strongly on the rotational stiffness of the monopilesoil system under un- and reloading. It was shown that the behavior under un- and reloading can be captured sufficiently accurate with a numerical model, which uses the HSsmall material law for the soil. This material law accounts for strain- and stress-dependent stiffness and considers, in conjunction with a Mohr-Coulomb failure criterion, plastic deformations by an isotropic strain-hardening approach. From parametric studies with this model, the following main conclusions can be drawn for monopiles in sand soils:

• The secant stiffness of the monopile for un- and reloading, which is the relevant stiffness in a dynamic calculation of the wind turbine structure, varies only slightly with the current mean load and the magnitude of the cyclic load span. Compared to the initial stiffness

72

Applied Ocean Research 84 (2019) 62–73

M. Achmus et al.

• •

determined as the initial slope of the load-deformation curve for monotonic loading, the secant stiffness slightly increases with increased mean load, but decreases with higher cyclic load spans. For most relevant cyclic load spans of less than 50% of the maximum load, the actual secant stiffness is at maximum 2% greater and 2.5% smaller than the initial stiffness. This shows that the initial stiffness of the monopile-soil system gives a good and in usual cases sufficiently accurate estimate of the relevant stiffness. Of course, in order to get a good prediction of eigenfrequencies of the wind turbine structure, a p–y method must be used which is suitable for the calculation of the initial stiffness of the monopile. It was shown that the ordinary API p–y approach yields a poor prediction of rotational stiffness, whereas both the Thieken approach and the Kallehave approach seem suitable for monopiles in sand soils. Compared to the numerical results, both approaches still underpredict the initial rotational stiffness. However, this might account for the fact that for greater cyclic load spans the secant stiffness is slightly smaller than the initial stiffness of the monopile system.

[10]

[11]

[12]

[13]

[14] [15]

[16] [17]

[18]

References [1] American Petroleum Institute - API, Recommended Practice 2GEO - Geotechnical and Foundation Design Considerations, (2014). [2] DNVGL, Support Structures for Wind Turbines Offshore, Standard DNVGL-ST-0126, DNV GL, Oslo, Norway, 2016. [3] D. Kallehave, C. Le Blanc Thilsted, M.A. Liingaard, Modification of the API p-y formulation of initial stiffness in sand, Proceedings of the 7th International Conference for Offshore Site Investigation and Geotechnics (2012). [4] M. Damgaard, M. Bayat, L.V. Andersen, L.B. Ibsen, Assessment of the dynamic behavior of saturated soil subjected to cyclic loading from offshore monopile wind turbine foundations, Comput. Geotech. 61 (2014) 116–126. [5] B. McClelland, J.A. Focht, Soil Modulus for Laterally Loaded Piles, Trans. ASCE, ASCE, Reston, VA, 1958, pp. 1049–1063. [6] L.C. Reese, W.R. Cox, F.D. Koop, Analysis of laterally loaded piles in sand, Proc. of the Offshore Technology Conference (OTC) (1974). [7] W.R. Cox, L.C. Reese, B.R. Grubbs, Field testing of laterally loaded piles in Sand, Proc. of the Offshore Technology Conference (OTC) (1974). [8] J.R. Murchison, M.W. O’Neill, Evaluation of p-y-relationship in cohesionless soils, in: J.R. Meyer (Ed.), Analysis and Design of Pile Foundations, ASCE, New York, 1984, pp. 174–191. [9] American Petroleum Institute - API, Recommended Practice for Planning, Designing

[19] [20] [21] [22] [23] [24]

[25]

[26]

73

and Constructing Fixed Offshore Platforms – Working Stress Design, 21st edition, (2007) and Errata and Supplements 1, 2 and 3. K. Abdel-Rahman, M. Achmus, Finite element modeling of horizontally loaded monopile foundations for offshore wind energy converters in Germany, Proceedings of the 1st International Symposium on Frontiers in Offshore Geotechnics (ISFOG) (2005) 391–396. J. Wiemann, K. Lesny, W. Richwien, Evaluation of the pile diameter effects on soilpile stiffness, Proceedings of the 7th German Wind Energy Conference (DEWEK), (2004). S.P.H. Sørensen, Soil-structure Interaction for Non-slender, Large-diameter Offshore Monopiles. PhD Thesis, Aalborg University, Denmark, Department of Civil Engineering, 2012. T. Hald, C. Mørch, L. Jensen, C. LeBlanc Bakmar, K. Ahle, Revisiting monopile design using p-y curves results from full scale measurements on Horns Rev, DONG Energy A/S, European Offshore Wind Energy Conference and Exhibition, (2009). C. LeBlanc, G.T. Houlsby, B.W. Byrne, Response of stiff piles in sand to long-term cyclic lateral loading, Gotechnique 60 (2) (2010) 79–90. K. Thieken, M. Achmus, K. Lemke, M. Terceros, Evaluation of p-y approaches for large diameter monopiles in sand, Int. J. Offshore Polar Eng. 25 (2) (2015) 134–144. K. Thieken, M. Achmus, K. Lemke, A new static p-y approach for piles with arbitrary dimensions in sand, Geotechnik 38 (4) (2015) 267–288. M. Achmus, K. Thieken, K. Lemke, Evaluation of a new p-y approach for piles in sand with arbitrary dimensions, 3rd International Symposium on Frontiers in Geotechnics (ISFOG) (2015). B.W. Byrne, R.A. McAdam, H.J. Burd, C.T. Houlsby, C.M. Martin, et al., PISA: New design methods for offshore wind turbine monopiles, Proceedings of the 8th International Conference Offshore Site Investigations and Geotechnics (OSIG) (2017). R.B.J. Brinkgreve, E. Engin, W.M. Swolfs, PLAXIS 3D Manual, (2014). W. Li, D. Igoe, K. Gavin, Field tests to investigate the cyclic response of monopoles in sand, Proc. Inst. Civil Eng. – Geotechn. Eng. 168 (GE5) (2015) 407–421. M. Achmus, K. Thieken, On the behavior of piles in non-cohesive soil under combined horizontal and vertical loading, Acta Geotech. 5 (3) (2010) 199–210. T. Benz, Small Strain Stiffness of Soils and Its Numerical Consequences. PhD Thesis, University of Stuttgart, 2007. T. Schanz, Zur Modellierung des mechanischen Verhaltens von Reibungsmaterialien, University of Stuttgart, Habilitation, 1998. J.A. Dos Santos, A.G. Correia, Reference threshold shear strain of soil and its application to obtain a unique strain-dependent shear modulus curve for soil, Proceedings of the 15th International Conference on Soil Mechanics and Geotechnical Engineering, (2001), pp. 267–270. J. Gattermann, V. Herwig, C. Moormann, VIBRO Project – Vergleich des lateralen Tragverhaltens von vibrierten und geschlagenen Stahlrohrpfählen in sandigen Böden, Proceedings of HTG Kongress 2015 in Bremen (2015) 173–184. R.B.J. Brinkgreve, E. Engin, H.K. Engin, Validation of empirical formulas to derive model parameters for sand, Proceedings of the 7th European Conference on Numerical Methods in Geotechnical Engineering (2010) 137–142.