Macro- and micro-mechanics of granular soil in asymmetric cyclic loadings encountered by offshore wind turbine foundations

Abstract

Offshore wind turbine foundations are subject to 10⁷–10⁸ cycles of loadings in their designed service life. Recent research shows that under cyclic loading, most soils change their properties. Discrete Element Modelling of cyclic simple shear tests was performed using PFC2D to analyse the micromechanics underlying the cyclic behaviours of soils. The DEM simulation were first compared with previous experimental results. Then asymmetric one-way and two-way cyclic loading pattern attained from real offshore wind farms were considered in the detailed parametric study. The simulation results show that the shear modulus increases rapidly in the initial loading cycles and then the rate of increase diminishes; the rate of increase depends on the strain amplitude, initial relative density and vertical stress. It shows that the change of soil behaviour is strongly related to the variation of coordination number, rotation of principal stress direction and evolution of degree of fabric anisotropy. Loading asymmetry only affects soil behaviours in the first few hundred of cycles. In the long term, the magnitude of (γ_max − γ_min) rather than loading asymmetry dominates the soil responses. Cyclic loading history may change the stress–strain behaviour of a soil to an extent dependent on its initial relative density.

1 Introduction

1.1 Background and motivation

Offshore wind turbine foundations are subjected to a combination of cyclic and dynamic loading arising from wind, wave, 1P (rotor frequency) and 2P/3P (blade passing frequency) loads. Designing foundations for offshore wind turbines is challenging as: (1) these are dynamically sensitive structures in the sense that natural frequencies of these structures are very close to the forcing frequencies [1]; (2) cyclic loading can induce accumulated rotation of foundation but the tolerance for the total rotation at seabed is suggested as low as 0.5° in DNV [2]. The soil-foundation interactions under cyclic lateral loading have been investigated intensively by field tests [3], scaled model tests [4,5,6,7,8,9,10], FEM simulations [11] and DEM simulations [12, 13] with majority focused on mono-pile foundations as it is the main type of foundation for offshore wind turbines. Empirical relationships for the foundation tilt or stiffness were derived from the model tests but few attention was drawn on the soil disturbances due to the limitations of monitoring technique. With the aid of numerical modelling, stress and strain field within the soils could also be visualised in addition to the foundation responses. But due to the high computational costs, the number of loading cycles in the numerical simulations for the entire soil-pile system is usually limited.

The cyclic lateral load considered in the literature is mainly symmetric two-way loading or simply loading/unloading. Only a few considered the effects of asymmetric cyclic loadings, e.g. Zhu et al. [7, 8, 14]. Zhu et al. [7] found that the accumulated rotation of suction caisson foundation increases with decreasing ζ_c (the ratio of minimum moment to the maximum moment), i.e. two-way cyclic loading causes more accumulated rotations. It contrasts with the findings by Zhu et al. [8], where one-way cyclic loading causes more accumulated rotations. It demonstrated that it is necessary to perform a systematic study on the influence of asymmetric cyclic loading on soil-structure interactions. On the other hand, Jalbi et al. [15] reviewed the wind and wave loading data from 15 offshore wind farms and determined the maximum and minimum moment applied at the mudline of a wind turbine. It was found that ζ_c lie in the range between − 0.5 and 0.5 for these wind farms. This is a combined effect of wind load with high magnitude but low frequency and wave load with low magnitude but high frequency. With increasing water depth, the magnitude of wave load increases and may change the one-way cyclic loading caused by wind load into two-way loading. A schematic diagram illustrating the variations is given in Fig. 1a.

Fig. 1

a Combination of wind and wave load, b schematic diagram of soil stress conditions surrounding a monopile and c one-way and two-way cyclic loading patterns adopted in current DEM simulations

To attain detailed information of soil-structure interactions under various conditions, it is better to consider a small zone of soils adjacent to the structure and replicate it in an element test with similar stress conditions and stress paths. Intensive investigations have been conducted using laboratory element tests to examine the accumulated deformation [16, 17] and variation of stiffness [18, 19] of soils under cyclic loading in different load conditions. However, these experimental studies have no access to the micro-mechanism underlying the soil responses. Discrete element modelling of cyclic tests were also widely adopted to link the macro-scale soil responses and underlying micro-mechanism during cyclic loading including [20,21,22]. However, these studies were limited to symetric two-way loadings, which are different from the real loading scenarios in offshore wind farms as found by Jalbi et al. [15]. This is a knowledge gap to fill in the current study.

1.2 Aim and scope of the paper

Majority of current offshore wind turbines are supported by monopile foundation, which is a large steel tube typically 30–40 m in length and 3–7 m in diameter. Unlike the slender piles for offshore structure, monopile tends to rotate rather than bend under lateral load or overturning moment. Therefore, the interactions between the monopile and the soil element in front of the monopile is analogic to a cyclic simple shear scenario, where the vertical stress remains almost constant, as shown in Fig. 1b. It is noted that the vertical stress applied on Block A may fluctuate slightly due to the compaction and rearrangement of soils above Block A, but this fluctuation is negligible comparing with the vertical stress levels considered in the current study. Therefore, cyclic simple shear tests were adopted to explore the soil responses under various cyclic loading patterns and its underlying micromechanism, in particular, to study the soil stiffness and deformation responses. The cyclic loading profiles considered are illustrated in Fig. 1c, which embrace the load scenarios summaries by Jalbi et al. [15] as well as previous studies.

Experimental cyclic simple shear tests on typical silica sand (RedHill 110) have been conducted in previous studies [23]. DEM simulations were first performed using PFC2D in the current study and compared with the experimental results. The micromechanical parameters were then analysed to find the relationships between the micromechanical parameters and macroscopic responses with particular focus on the influence of one-way and two-way loading patterns. Finally, original samples and cyclically sheared samples are simply sheared till the critical state. Variations of soil stress–strain response and volumetric response as a consequence of cyclic loading history were analysed under the critical state framework.

2 Description of DEM simulation programme

In the experimental study performed by Nikitas et al. [23], a cyclic simple shear apparatus was used for testing cylindrical soil samples with 50 mm in diameter and 20 mm in height as suggested in ASTM D6528 [24]. RedHill 110 Sand, a poorly graded fine grained silica sand with d₅₀ = 0.18 mm and particle size distribution (PSD) curve shown in Fig. 2, was tested as this soil has been used to carry out scaled model tests on different types of foundations [25, 26]. Strain controlled drained cyclic simple shear tests with symmetric two-way loading on sand with various relative densities, vertical stresses and shear strain amplitudes were performed.

Fig. 2

Particle size distribution (PSD) of test sand in experiments [23] and current DEM simulations

In the current study, a commercial DEM code PFC2D (Itasca, 2008) was used to perform the presented simulations. The experimental cyclic simple shear test is a three-dimensional problem; however the current study only models a thin slice of the sample in the middle. It is obvious that a two-dimensional simulation cannot accurately represent the three-dimensional granular soil. However, there is no intention in this paper to reproduce the physical test quantitatively, but analyse the similar underlying micro-mechanism because the major and minor principal stresses in a simple shear test lie in the loading plane considered in the 2D simulations.

The sample initially generated for testing is about 20 mm in height and 50 mm in width, similar to the sample dimensions in experiments [23]. It contains 8000 disks with size ranging from 0.1 mm to 0.3 mm and d₅₀ = 0.18 mm, matching the d₅₀ value in experiments. The PSD curve is also given in Fig. 2. Note that the PSD range is narrower in the DEM than that in experiments due to the unrealistic computational time required for simulating samples with the same particle grading as in experiments. Parameters used in DEM simulations are listed in Table 1. Two sets of samples were generated with radius expansion approach followed by 1D-consolidation (laterally confined) with particle friction coefficient, μ = 0 and 1.0, to generate relative dense and loose samples, respectively. Particle friction coefficient was then changed to 0.5 and specimens were brought to equilibrium again with the new frictional coefficient. Final void ratios (e) of the relative dense and loose samples at two different vertical stresses (σ) are also listed in Table 1. Note that due to different boundary velocity in consolidation stage, the loose sample reached slightly lower e at 50 kPa than that at 100 kPa.

Table 1 DEM simulation parameters

In the modelling of the cyclic simple shear test, a drained condition is maintained, i.e. the vertical normal stress is kept at a prescribed value by stress-controlled top and bottom walls (moving in or out to maintain the required stress level). The left and right side walls were rotated about their midpoints to realise the simple shear condition. The DEM simulation program of the cyclic simple shear tests is provided in Table 2. Symmetric two-way loading were planned for comparison with experimental tests [23] and other experimental tests and DEM simulations in literature. Various soil packing densities and vertical stresses (Series A), various strain magnitudes (Series B) were considered. In addition, the effect of asymmetry of cyclic loading (i.e. one-way or two way loading) were investigated. Similar to the parameter, ζ_c, used by Zhu et al. [7], in this paper, strain ratio, η_c = γ_min/γ_max, was defined to quantify the degree of asymmetry of cyclic loading as strain-controlled cyclic tests were performed in the current study. Strain ratios, η_c, between 0.5 and − 1.0 were considered to embrace the load ratios encountered in real wind farms as summarised by Jalbi et al. [15] and previous work. For easy comparison, one series of simulations were planned with the same value of γ_max − γ_min but various η_c (Series C), and two other series were planned with the same γ_max but various η_c (Series D and E).

Table 2 DEM simulation program

Twelve measurement circles, as depicted in Fig. 3, were defined within the sample to measure the average stress, void ratio and coordination number. The results presented in the following sections are the average value from these twelve measurement circles.

Fig. 3

Schematic diagram of soil sample and locations of measurement circles

3 DEM simulation results

3.1 Macro-scale responses

The shear stress-shear strain curve forms hysteresis loops during cyclic loadings. Stress–strain relationships for the symmetric two-way loading (Series A, B, C) are similar. Thus only results from Series E were presented in Fig. 4. While η_c decreases from 0.5 to − 1.0 (one-way loading to symmetric two-way loading), stress–strain relationships also show different evolution trends. With η_c = 0.5, shear strain only reverses half way, shear stress reduces to a positive value slightly above zero in the first cycle. With the progress of cyclic loading, shear stress at both γ_max and γ_min deceases steadily. As a consequence, the negative τ_min at γ_min approaches the same magnitude of the positive τ_max at γ_max, even though both γ_max and γ_min are positive values. With decreasing η_c, the τ_min decreases more as γ_min decreases. However, the trends of τ_max are quite different, even though γ_max remains the same. When η_c > 0, τ_max decreases with cyclic loading; when η_c = 0, τ_max remains almost constant during cyclic loading; when η_c < 0, τ_max increases with cyclic loading. It demonstrates that the magnitude of γ_min affected the stress response at γ_max. Following 1000 number of loading cycles, no matter the symmetry of loading cycles (one-way or two-way), shear stress oscillates two-way symmetrically, i.e. negative τ_min at γ_min reaches the same magnitude as the positive τ_max at γ_max.

Fig. 4

Stress–strain relationships during cyclic loadings for simulation Series E

The shear modulus (G) of the sample can be calculated as

$$G = \frac{{\tau_{\mathrm{max} } - \tau_{\mathrm{min} } }}{{\gamma_{\mathrm{max} } - \gamma_{\mathrm{min} } }}$$

(1)

The variations of G in all simulations under cyclic loading are illustrated in Fig. 5. The magnitudes of shear moduli in the DEM simulations (e.g. 5–6.5 MPa for σ = 100 kPa, dense sample, γ = (− 0.52%, 0.52%)) are in the same range as the values from the experimental tests [e.g. 4–6 MPa for σ = 100 kPa, relative density = 50%, γ = (− 0.5%, 0.5%)] in Nikitas et al. [23]. Referring to Fig. 5a, for both loose samples, there is a clear increase in G under cyclic loading; while for both dense samples, G decreases obviously, which was not observed in the experiments. Reason for this could due to the irregular shapes of particles, wide particle grading and possible particle abrasion or crushing in experiments, which could lead to continuously densification of samples. However, these phenomena could not be easily replicated in DEM simulations. After 6000 cycles, G of the loose samples and the dense samples at a same σ approach a same constant. As shown in Fig. 5b, G increases dramatically with reducing γ_max as expected, and in all cases G reduces slightly during cyclic loading. The results correlated quite well with the observations from scaled model tests with different types of offshore wind turbine foundations [4,5,6, 27] and the field measuremens [28].

Fig. 5

Evolution of shear modulus during cyclic loading in DEM simulations

With the same magnitude of (γ_max − γ_min) in Series C, the two-way loading resulted in higher G than the one-way loading in the first few hundred cycles (Fig. 5c). The reason is that the true strain level for the one-way loading is larger than those for the two-way loading (i.e. 1.04% > 0.75% > 0.52%), and as expected G decreases with increasing strain level. However, following large number of loading cycles, soil losts memory of initial strain levels and the key parameter dominating the long term stiffness is the magnitude of (γ_max − γ_min), which are the same for the three simulations, so are the final G values. Similar phenomena also observed in the DEM simulations of pile-soil interaction under cyclic loading by Cui and Bhatacharya [12]. The main reason for this may be that asymmetry only appears in the first cycle and the remaining cycles appear to be pseudo-symmetric about a non-zero mean shear strain, (γ_max + γ_min)/2. The different soil responses caused by asymmetric loading in the first cycle are gradually eliminated by the pseudo-symmetric loading. Samples seem to reach a cyclically stable state with response center drifted to the non-zero mean shear strain and no obvious differences between samples can be observed.

With the decreasing values of η_c in Series D and E, G decreases obviously as the magnitude of (γ_max − γ_min) increases (Fig. 5d and e), which confirms the findings from Series C that the key parameter dominating the long term magnitude of stiffness is the magnitude of (γ_max − γ_min). It can also be observed that G decreases in dense samples and increases in loose samples in general. The exemptions are the dense samples with η_c > 0, where the magnitude of (γ_max − γ_min) may be not large enough to produce the same effect.

The evolutions of void ratio, e, of all simulations under cyclic loading are illustrated in Fig. 6. The increasing G of loose samples could be explained by the densifications of samples (reduction in e) as observed in Fig. 6a. The G values for the dense samples reduce obviously, but they only dilated slightly. There should be other causes for the remarkable decrease of G, which will be investigated further in the latter part of this paper. Moreover, e of the two samples at σ = 50 kPa coincide, which consists with the coincidence of their G values. However, for the two samples at σ = 100 kPa, e of the loose sample is still higher than that of the dense sample, which is also consistent with the comparison of their G values. As shown in Fig. 6b, with higher γ_max, dense sample dilates slightly more in the initial stage. However, for γ_max = 1.04%, sample dilates first and then contracts slightly. As seen from the monotonic simple shear tests in the latter section (Fig. 14), this is mainly because the specimen approaches failure from 2% of shear strain thus the sample is disturbed and re-arranged significantly during cyclic loadings. As observed from Series C in Fig. 6c, the three simulations with the same value of (γ_max − γ_min) but different values of η_c shows very subtle differences in void ratio. However, with increasing value of (γ_max − γ_min) in Series D/E (Fig. 6d, e), void ratio reduces/increases more obviously. This confirms again that the magnitude of (γ_max − γ_min) rather than η_c dominates the responses.

Fig. 6

Evolution of void ratio during cyclic loading in DEM simulations

3.2 Micro-scale mechanism

The observed macro-scale stress and strain responses should be underlain by the micro-scale (particle-scale) mechanism. In the following section, micro-scale parameters, including the coordination number, contact force network, fabric anisotropy, principal stress rotation, etc., were examined in details to bridge the micro–macro gaps.

3.2.1 Coordination number (N_c)

The coordination number (N_c) is the average number of contacts surrounding each particle. It has a strong relation with the stress level within the sample [29]. The evolutions of N_c under cyclic loading for the four samples under various test conditions are shown in Fig. 7. It is clearly shown in Fig. 7a that the initial low N_c corresponds to initial low shear stress (thus low G in Fig. 5a). The increase in G for the two loose samples is related to the increase in N_c and the decrease in G for the dense samples agrees with the reduction in N_c. When values of G coincide, the coordination numbers also coincide. As observed in Fig. 7b, with higher γ_max, N_c reduces more and quicker with loading cycles, matching the reduction of G. As shown in Fig. 7c, with the same (γ_max − γ_min), the initial value of N_c was lowest with η_c = 0 (one-way loading) due to higher γ_max; however after 100 cycles, N_c decreases to similar value for all three η_c values, which agrees with the trend for the G in Fig. 5c. As seen from Fig. 7d, N_c does not change significantly for η_c > 0 but decrease more significantly with decreasing η_c for η_c ≤ 0, which agrees with the observations of G in Fig. 5d and e in Fig. 6d. N_c for Series E all increases remarkably with cyclic loading but the differences between various η_c are insignificant.

Fig. 7

Evolution of coordination number

In summary, considering all observations of G, e and N_c, we can see that lower e (denser packing) results in higher N_c, which in turn leads to higher stress level, thus higher G.

3.2.2 Contact force network

The evolution of the contact force network within the dense soil sample at σ = 100 kPa is illustrated in Fig. 8. Before cycling, the contact force network illustrates an approximately isotropic and uniform distribution (Fig. 8a). Once cycled to the maximum strain, e.g. + 0.52%, the stronger contact forces are aligned to the shorter diagonal direction and weaker contact forces can be observed at the other two diagonal corners (Fig. 8b). When the specimen is cycled to − 0.52%, the strong contact force direction rotated to the other diagonal direction, and the locations of the weaker contact force network also shifts accordingly (Fig. 8c). To statistically quantify the orientations of contact forces, the contact fabric is analysed Sect. 3.3.3; to statistically quantity both the magnitudes and orientations of the contact forces, the principal directions were examined in Sect. 3.3.4.

Fig. 8

Evolution of contact force network during the first simple shear loading cycle (dense sample with γ = (− 0.52%, 0.52%) and σ = 100 k)

3.2.3 Fabric

The spatial distribution of the contact force directions can be quantified by fabric. There are many evidences of the impact of fabric anisotropy on the characteristics of granular materials [30, 31]. Therefore, it is worth to analyse the evolution of soil fabric in the current cyclic loading conditions. Various definitions of fabric can be found in the literatures [32]. The one adopted in the current study is the Fourier approximation [33], which quantifies the distribution of contact direction per radian as:

$$E\left( \theta \right) = \frac{1}{2\pi }\left[ {1 + a\cos 2\left( {\theta - \theta_{a} } \right)} \right]$$

(2)

where a is a parameter defining the magnitude of fabric anisotropy and θ_a is the direction of the principal fabric. For an isotropic sample, a = 0 and E(θ) =1/2π, which is a circle with a uniform distribution of 1/2π per radian.

The histogram of spatial distribution of directions of particle contact normals for the dense sample at σ = 100 kPa is illustrated in Fig. 9. The Fourier approximation function is indicated by the red ellipse with the long axis indicating the major principal direction of fabric (θ_a). The major principal fabric direction of the initial sample is 92°. When sheared to the maximum strain (0.52% or 0.3°), the major principal fabric direction rotates to the diagonal direction (θ_a≈ 130.0°); when it is sheared to minimum strain (− 0.52% or − 0.3°), the major principal fabric direction rotates to the perpendicular diagonal direction (θ_a≈ 40.0°). Note that the sample boundaries are only rotated up to 0.6°; however, the major principal fabric direction rotates about 90°.

Fig. 9

Spatial distribution of contact normals (dense sample with γ_max = 0.52% and σ = 100 k)

The evolutions of θ_a at γ_max and γ_min for Series A and D are plotted in Fig. 10; the other Series show similar trends to Series A thus not included here. Note that the initial θ_a at γ = 0 for the four samples are 83° (loose, σ = 50 kPa), 165° (loose, σ = 100 kPa), 80° (dense, σ = 50 kPa) and 92° (dense, σ = 100 kPa). It can be observed from Fig. 10a, that the rotation of θ_a occurs slowly in loose samples. For the loose sample with σ = 50 kPa, θ_a at γ_max starts from about 110° and increases slowly to 130° at large cycle number; while for the loose sample with σ = 100 kPa, θ_a at γ_max starts from about 140° and decreases slowly to 130°. However, for both dense samples, θ_a at γ_max reaches 130° in the first loading cycle and remains at 130°. Rotation of θ_a at γ_min shows similar trend: rotation of θ_a occurs slowly in loose samples until it reaches 40°; θ_a in dense samples arrives at 40° in the first cycle. In the other simulation Series, rotation of θ_a does not show obvious difference with various γ_min or η_c except for the Simulation D-6 with η_c = 0.5, where θ_a at γ_min does not reduce to 40° as other simulations, but stays at 130° for γ_min for about 10 cycles, then reduces slowly to 65° following 6000 loading cycles. The stabilised fabric following cyclic loading is not a unique feature for drained tests, but also observed in undrained cyclic tests [34, 35]. It is also interesting to observe that θ_a at γ_max and θ_a at γ_min do not symmetric about 90°, but 85° for all simulations, despite their loading asymmetry. One hypothesis for this phenomenon is the rotation direction in the first quarter of the first cycle. Further investigation will be carried out in the future to verify this hypothesis.

Fig. 10

Evolution of major fabric direction at maximum strain and minimum strain

The magnitude of fabric anisotropy is quantified by a in Eq. (2). It has been demonstrated previously [30, 31] that larger magnitude of anisotropy can result in higher shear stress. As illustrated in Fig. 9, a relates to the aspect ratio of the red ellipse. At γ_max, the long axis is along 130° and short axis is along 40°; at γ_min, the long axis and short axis swop locations. Thus, the variation of fabric anisotropy from γ_max to γ_min can be quantified by the change of axis length along one principal direction, or a_max+ a_min. On the other hand, fabric anisotropy is an indicator of shear stress magnitude, while G is determined by (τ_max − τ_min)/(γ_max − γ_min); therefore, a_max+ a_min could be an indicator of G for the same level of (γ_max − γ_min). The evolutions of a_max+ a_min for the five simulation Series are shown in Fig. 11. It is clear in Fig. 11a that a_max+ a_min for the dense samples drops significantly, which agrees with the fact that G for these samples also drops obviously. And a_max+ a_min for the loose samples increases clearly, which matches the increase in G (Fig. 5a). With increasing γ_max, degree of anisotropy increases dramatically (Fig. 11b), which reflects the increasing shear stress level. However, as (γ_max − γ_min) also increases significantly, G reduces. In Series C, (γ_max − γ_min) remains the same, thus a_max+ a_min for η_c = 0 is lower than the other two cases in the first 100 cycles and then approaches similar values with other two cases (Fig. 11c), which agrees with the trend for G for this Series (Fig. 5c). In Series D and E, a_max+ a_min increases with decreasing η_c, which results in increasing (τ_max − τ_min) (Fig. 4); however, as (γ_max − γ_min) also increases with decreasing η_c, the resultant G is actually decreasing. Thus, change of G is a combined result of variation of degree of anisotropy and strain amplitude, and in Series D&E, strain amplitude has a greater impact than the degree of anisotropy.

Fig. 11

Difference of magnitude of fabric anisotropy between maximum strain and minimum strain

3.2.4 Rotation of principal directions

The average stresses within the specimen during cyclic loading were monitored and the principal stresses and principal directions were determined. It is found that the major principal directions and the major fabric directions incline similar angles to the horizontal: 130º for both the major principal directions and the major fabric directions at γ_max; 40º for the major fabric directions and 50º for the major principal directions at γ_min. The rotations of the major principal directions in the first three cycles were illustrated in Fig. 12. Note that the initial major principal direction of the loose sample with σ = 100 kPa has an approximate horizontal orientation, while the remaining three samples have approximate vertical orientations. The rotation of the major principal direction of the loose sample with σ = 50 kPa is slower than the other three samples in the first three cycles (Fig. 12a). The rotation of the major principal direction for various γ_max yields same values at maximum and minimum strain, but forms larger loops for larger γ_max (Fig. 12b). Apart from the loose sample with σ = 50 kPa, for all the remaining samples with symmetric two-way loading as shown in Fig. 12a–c, majority of the principal direction rotation occur at about half the γ_max (γ_min) in the loading (unloading) stages. Outside the strain range of (γ_min/2, γ_max/2), principal directions remain almost constant. Due to this feature, the rotation of principal direction with η_c > 0 in Series D and E shows some interesting characteristics. In particular, in Simulation D-6, where η_c = 0.5, as γ_min is only half of γ_max, rotation of principal direction has not start yet before the loading changes direction in initial cycles. As a result, major principal direction remains at about 130° in the first three cycles (black curves in Fig. 12d). Rotation of principal direction for this sample gradually started after ten cycles and approaches 80° at the end of 6000 loading cycles. This delayed rotation of principal stress direction agrees with the delayed rotation of principal direction of fabric (Fig. 10b). In Simulation E-6 (loose sample with η_c = 0.5), as sample is relative loose and there are more spaces for particle movements and rearrangements, even the γ_min is only half of γ_max, rotation of principal direction starts from the first cycle to a smaller extent (Fig. 12e). But the shapes of the direction-strain loop are still different from majority simulations: the delay of rotation when loading direction changes at γ_min is not observed in Simulation D-6 and E-6. This is because the principal direction has yet stabilised at γ_min as other simulations with lower γ_min.

Fig. 12

Rotation of major principal direction in the first three loading cycles

3.2.5 Particle movements

Particle movements during 6000 loading cycles for a few representative simulations are displayed in Fig. 13. It can be observed clearly from Fig. 13a that for the dense sample with η_c = 0.5 (Simulation D-6), the particles flow anticlockwise slightly within the box due to the low value of (γ_max − γ_min), i.e. low intensity of vibration. As the η_c decreased to − 1.0 (Fig. 13b), i.e. (γ_max − γ_min) increases, particle movements are more remarkable and there are local vortices formed. For both simulations, no clear volumetric change can be observed, which is consistent with the void ratio observation. However, for the two simulations of loose samples (Fig. 13c and d), in addition to the local vortices, clearly inward particle movements can be observed along the top and bottom boundaries, which is again consistent with the observations of void ratios. Moreover, the particle displacements are more intensive with increasing (γ_max − γ_min), because there are more voids within loose samples allowing particle movement and rearrangements.

Fig. 13

Particle convections in 6000 loading cycles (amplified by 2.0)

3.3 Effect of cyclic loading on p-y curve and critical state analysis

In the design of monopile for the offshore wind turbine, soils surrounding monopile may be replaced by a series of independent springs with their reaction forces applied to monopile specified by p–y curves. The p–y curve could be attained from the conversion of a τ–γ curve [36]. Therefore, it is worth to check the variations of τ–γ curve (thus p–y curve) following 6000 cycles of loading.

Figure 14 shows the normalised τ-γ curves for the original dense and loose samples (without cyclic loading) with σ = 50 kPa and 100 kPa as well as the evolutions of e during monotonic simple shear tests until shear strain γ = 52%. Dense samples show higher initial stiffness and higher peak stress ratios as expected. Specimen with similar initial e shows slightly higher peak stress ratio and initial stiffness under lower σ. The e responses of dense samples showed clearly dilation and finally approached the critical e, while loose samples also showed slightly dilation.

Fig. 14

Monotonic simple shear tests of original samples

Similar monotonic simple shear tests were also performed on those four samples followed 6000 cycles of symmetric loadings with shear strain amplitude γ = ± 0.52%. The stress–strain and volumetric responses during monotonic simple shear test were illustrated in Fig. 15. Comparing Fig. 15b with Fig. 6a, it is interesting to observe that loose specimen densifies during cyclic loading, with e approaches constants similar to the values for dense samples; however, these values are not the critical e. The cyclic loading equalises the sample density which leads to very similar stress–strain response and volumetric responses of these four samples as observed in Fig. 15.

Fig. 15

Monotonic simple shear tests of cyclically sheared samples

To have a better understanding of the relationship of the stress ratios and e at the peak and critical states, the stress paths for the eight simulations as well as the critical state line (CSL) were plotted in Fig. 16, where the mean stress s′ = (σ₁′ + σ₂′)/2, the deviator stress q = σ₁′ − σ₂′ and σ₁′, σ₂′ are the principal stresses. The CSL is obtained by least-square fitting of the e-log(s′) values or the q–s′ values of the eight simulations at the critical state. It can be observed that the e for one relative loose sample (σ = 100 kPa) is just below the critical state e, while the other loose sample (σ = 50 kPa) is below further (Fig. 16a). The e values for the two dense samples are well below the critical state e. In the simple shear tests, their mean stresses first increase to the peak state value with small increases of e, then decrease significantly with further increases of e until reaching the critical state line. It can be seen that the loose sample with σ = 100 kPa has similar response to real loose sand; the loose sample with σ = 50 kPa has similar response to medium sand; and the two dense samples behave similar to very dense sand. As shown in Fig. 16b, these four samples reach different peak failure envelops due to various initial e values.

Fig. 16

Stress path for original samples and cyclically sheared samples

When the two loose samples subjected to 6000 cycles of symmetric shearing loadings, their e values reduce significantly with mean stresses remain unchanged; while the two dense samples remain similar e and mean stresses during cyclic loading (Fig. 16c). When the cyclically sheared samples are subjected to simple shear tests, all four samples have similar responses: reaching a same peak failure envelop due to similar e and then dropping onto the critical state line (Fig. 16d). Therefore, cyclic loading equalises the distances between the sample state points and the critical state line. Based on the critical state soil mechanics theory, the distance between a soil state and the critical state line dominates soil responses, thus the four samples post cyclic loading show very similar responses.

4 Conclusions

Series of DEM simulations of cyclic simple shear tests with various strain profiles (one-way or two-way) were performed to explore the variations of soil characteristics under cyclic loading and its underlying micro-mechanism. It has been found that:

• Shear modulus for loose soil increases rapidly in the initial loading cycles as a result of soil densification, and then the rate of increase diminishes when void ratio approaches a constant; shear modulus of the very dense samples in DEM decreases but dilation is not the main reason.

• Under the same vertical stress and strain amplitude, shear moduli and void ratios of dense sample and loose sample approach the same values at larger number of cycles.

• Shear modulus increases with increasing vertical stress and relative density, but decreases with increasing strain amplitude as expected.

• Higher shear stress level and shear modulus are correlated to higher coordination number and higher magnitude of fabric anisotropy.

• Accumulated rotation of principal direction of fabric and major principal direction occurs slowly in loose sample, but reaches the final value in the first cycle in dense sample. In each symmetric loading cycle, majority of the principal direction rotation occurs between half of γ_min and half of γ_max in spite of the magnitude of γ_max. There are delays in the reverse of principal direction rotation immediately after the boundary rotation is reversed apart from the case with η_c = 0.5 where the vibration intensity is low and sample has not stabilised at γ_min.

• Loading asymmetry only affects soil behaviours in the first hundreds of cycles. In long term, the magnitude of (γ_max − γ_min) rather than loading asymmetry dominates the soil responses.

• Changes in the stress–strain behaviour (thus p-y curve) of soil due to cyclic loading are significantly dependent on the initial relative density and the distance between the state point and the critical state line.

This study verifies the capability of DEM in analysing the micro-mechanism underling the soil stiffness variation during cyclic loading. More studies will be carried out in the future work, in particular, non-circular disks will be considered to explore the particle shape effect. An empirical model for evolution of shear modulus will be established to include the impact of all macroscopic and microscopic parameters. These findings provide basis for improving the design codes for offshore wind turbine foundations.