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

Offshore wind turbine foundations are subject to 107–108 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.


Abbreviations a
Parameter defining the magnitude of fabric anisotropy a max Fabric anisotropy at maximum shear strain a min Fabric anisotropy at minimum shear strain d 50 Diameter at which 50% of the sample's mass is comprised of particles with a diameter less than this value e Void ratio G Parameter defined to quantify the degree of asymmetry of cyclic loading θ a Direction of the principal fabric σ Vertical stresses on the sample in a cyclic simple shear test τ max Shear stress at maximum shear strain τ min Shear stress at minimum shear strain 1 Introduction

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 1 3 73 Page 2 of 20 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.
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 macroscale 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.

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.

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 50 = 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.
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 micromechanism 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 50 = 0.18 mm, matching the d 50 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.
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).
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.

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 (oneway 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 . The shear modulus (G) of the sample can be calculated as 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 Simulation ID (γ min , γ max ) σ (kPa) Packing density 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]. With the same magnitude of (γ max − γ min ) in Series C, the two-way loading resulted in higher G than the oneway 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.

Number of Cycles
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.

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.

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.
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.

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.

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:  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 (2) E( ) = 1 2 1 + a cos 2 − a 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°. 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.
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 Fig. 9 Spatial distribution of contact normals (dense sample with γ max = 0.52% and σ = 100 k) 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.

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 .

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  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.

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. 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. 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.
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′ = (σ 1 ′ + σ 2 ′)/2, the deviator stress q = σ 1 ′ − σ 2 ′ and σ 1 ′, σ 2 ′ are the principal stresses. The CSL is obtained by leastsquare 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.
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.

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.