DEM simulations of transverse pipe–soil interaction on sand

Offshore oil and gas pipelines in deep water are generally laid directly on the seabed, since trenching is uneconomic. During operation these pipelines can be subject to repeated cycles of pressure- and temperature-induced axial expansion. This expansion can be accommodated through the formation of controlled lateral buckles at pre-determined locations, where the lateral displacements can be up to 10 or 20 pipe diameters (Bruton et al., 2007). Periods of production and shutdown give rise to cyclic lateral movements that are typically repeated hundreds of times during the life of the pipeline (Cheuk et al., 2008). This leads to the formation of soil mounds or ‘berms’ that progressively increase in volume as the pipe sweeps back and forth across the seabed, increasing its embedment while the seabed geometry continuously changes. Assumptions concerning the mobilised soil resistance are important in numerical modelling of controlled lateral buckling, which is required for a cost-effective pipeline design. Consequences of excessive lateral displacements include large bending strains in the pipe at the crown of a buckle, which may lead to significant fatigue damage. Recent incidents involving pipeline failure as a consequence of lateral buckling have occurred off the coast of Brazil (Almeida et al., 2001; Pasqualino et al., 2001), as well as in the North Sea and offshore West Africa (Bruton & Carr, 2011).

Soft normally consolidated clay is the most prevalent soil type for offshore developments in deep-water regions. Hence industrially relevant research on transverse (i.e. vertical and lateral) pipe–soil interaction has been dominated by modelling of pipeline response on clay, with most studies placing emphasis on undrained conditions. Physical modelling investigations of pipes on clay have included large-scale tests (Brennodden et al., 1986; Wagner et al., 1989; Verley & Lund, 1995; Cardoso & Silveira, 2010) as well as centrifuge model tests (Cheuk et al., 2007; Dingle et al., 2008; Cheuk & White, 2011). Numerical analyses, mainly using the continuum-based finite-element method, have been performed to simulate both small-amplitude (Aubeny et al., 2005; Merifield et al., 2008; Martin & White, 2012) and large-amplitude pipe displacements (Hesar, 2004; Yu & Konuk, 2007; Merifield et al., 2009; Wang et al., 2010; Chatterjee et al., 2012; Martin et al., 2013; Kong et al., 2018). These approaches can provide good predictions of pipeline response under a range of relevant conditions, and some studies have successfully replicated experimentally observed failure mechanisms in undrained clay soils.

Until now research on buckling of subsea pipelines (e.g. in the SAFEBUCK project) has tended to focus on deep-water clayey soils. However, pipe–soil interaction on sand is of primary concern for several facilities in the North Sea (siliceous sand), the North West Shelf of Australia, Timor Sea, South China Sea and Middle East (calcareous sand), as well as offshore Brazil (calcareous sand and clay). Owing to the fundamentally different mechanical behaviour exhibited by sand and clay soils under short-term loading, the above-mentioned physical and numerical modelling studies of pipelines on undrained clay contribute little to the understanding of transverse pipe–soil interaction on sandy seabeds.

Relevant studies involving physical modelling of pipelines on sand include in situ tests (Lambrakos, 1985), large-scale tests (Brennodden et al., 1986; Wagner et al., 1987; Palmer et al., 1988; Allen et al., 1989), centrifuge model tests (Zhang et al., 1999, 2001; Tian & Cassidy, 2010) and 1g model tests (Sandford, 2012; Li & Byrne, 2014). Originally, the motivation for performing lateral loading tests on pipe segments was to understand pipeline response during hydrodynamic loading. Indeed, the initial tests involved lateral displacements of less than one pipe diameter, D. In recent years, the research focus has shifted towards the investigation of pipe–soil interaction during lateral buckling. Hence pipe segments on sand have been tested to larger lateral displacements of up to 2D (Tian & Cassidy, 2010) and 4D (Sandford, 2012; Li & Byrne, 2014).

A number of semi-empirical load–displacement relationships have been developed to describe the soil resistance mobilised during lateral movement of a pipeline resting on sand. These relationships typically consist of separate terms that account for the frictional resistance to sliding, and the lateral resistance offered by a wedge of soil in the passive state (Wagner et al., 1987; Verley & Sotberg, 1994). The model proposed by Verley & Sotberg is widely used in industry and is still recommended by current design codes (DNV GL, 2017). It accounts for the initial lateral ‘breakout’ force and a constant residual force developed immediately after breakout, but it neglects the increasing soil resistance developed during large-amplitude cyclic lateral displacements of the pipe. More recent works have attempted to describe the transverse load–displacement response in the framework of macro-element plasticity modelling, which has previously been used for various types of shallow foundation on sand (Schotman, 1989; Tan, 1990; Nova & Montrasio, 1991; Gottardi et al., 1999). Following an early demonstration of the macro-element approach for pipes on clay (Schotman & Stork, 1987), macro-element models for pipes on sand have been developed by, among others, Zhang et al. (2001, 2002a, 2002b), Tian & Cassidy (2011) and Sandford (2012). Some of these models have also been incorporated into structural analysis programs (Tian & Cassidy, 2008; Sandford, 2012). Most of these models can replicate transverse pipe–soil interaction behaviour with reasonable accuracy for specific sands, but only for limited lateral pipe displacements. For instance the lateral displacement during the validation tests of Zhang et al. (2002b) reached only 0·2D. The model of Tian & Cassidy (2011) was used for pipe displacements up to 2D, but this value is still somewhat smaller than the typical displacements experienced by laterally buckled pipelines on the seabed; their model validation was also limited to monotonic pipe motions. Sandford (2012) validated his macro-element model against 1g model tests involving cyclic lateral displacements of up to 4D. Encouraging results were obtained, but only a few such analyses were performed.

Further investigation of pipe–soil interaction on sand is necessary to overcome the shortcomings of current load–displacement models, and to obtain a more comprehensive understanding of pipelines undergoing large-amplitude cyclic lateral displacements on sand. Numerical analyses are potentially a valid and cheaper alternative to physical modelling, where test repeatability and high-quality data are often difficult to achieve. So the aim of the numerical modelling undertaken here is on one hand to predict quantitatively the development of lateral pipe–soil resistance during monotonic and cyclic loading, and on the other hand to capture qualitatively the kinematics at the grain scale, with a special focus on the development of failure mechanisms. Both aspects are challenging because of the highly irregular geometry of the soil surface adjacent to the pipe and its profound modifications over time, for example when several berms are progressively formed by the cyclic lateral movement of a buckling pipeline scraping the seabed, and due to the formation of extensive localised deformations (shear bands) in the soil.

Numerical analyses of transverse pipe–soil interaction on sand using conventional Lagrangian (or updated Lagrangian) finite-element methods are difficult because of the complex constitutive behaviour of sand at low stress levels, the need for periodic remeshing (Sandford, 2012) and the presence of localised deformations. The achievement of mesh-independent results requires the use and calibration of sophisticated constitutive models (e.g. non-local models), while remeshing is necessary to deal with large plastic soil deformations and the evolution of the soil surface profile. Although coupled Eulerian–Lagrangian finite-element analysis is growing in popularity and avoids the need for remeshing, it is computationally expensive even for pipes on undrained clay (Martin et al., 2013; Kong, 2015) and it still suffers from mesh-dependency if a simple sand model is used, such as a Mohr–Coulomb material with non-associated flow.

The distinct-element method (DEM) takes advantage of the particulate nature of granular soils (Cundall & Strack, 1979). It can readily simulate problems where large plastic deformations occur in such soils, and can provide insights into the statics and kinematics occurring at the grain level. Its use was initially restricted to simulations of element tests such as the triaxial test, with a focus on investigating fundamental aspects of soil response. In recent years, DEM simulations of field-scale geotechnical problems have become more common – for example, stability of jointed rock masses (Boon et al., 2014), earth pressure of retaining walls (Widuliński et al., 2011), tunnelling (Bym et al., 2013; Boon et al., 2015), offshore wind turbine monopiles (Cui & Bhattacharya, 2016), cyclic loading of railway ballast (Chen et al., 2015) and cone penetration testing (Arroyo et al., 2011; McDowell et al., 2012).

DEM simulations of a partially buried pipeline on sand were performed for the first time by Macaro et al. (2014). In that work, DEM was used to analyse a pipe segment interacting with a three-dimensional (3D) prismatic particle domain, and the results were compared with experimental test data from Sandford (2012). The pipe displacement history was the same in all of the simulations, consisting of an initial small vertical penetration followed by a lateral movement at constant embedment. The work was restricted to monotonic displacement-controlled analyses, with the pipe reaching a lateral displacement of only 1D.

In this paper a more extensive programme of DEM simulations is reported, considering the behaviour of a pipe segment on sand undergoing vertical penetration followed by various lateral displacements (small and large, monotonic and cyclic) with different types of control in the vertical direction. The main objectives are: (a) to study the pipe–soil interaction during monotonic and cyclic displacement-controlled analyses; (b) to validate the DEM approach for simulating cyclic lateral displacements under constant vertical load, V; and (c) to study the influence of the vertical loading history experienced by the pipe prior to the onset of cyclic lateral displacements.

The DEM simulations were performed using the open-source code YADE (Kozicki & Donzé, 2008; Šmilauer et al., 2015). As part of this study, new user-defined algorithms to generate a cylindrical object representing the pipe segment, and to treat its mechanical interaction with the soil particles, were implemented in the code. A full description of the DEM model, as well as the calibration and validation procedures that were employed, can be found in Macaro (2015). The main features of the soil and pipeline modelling are outlined below, including the interparticle contact law, numerical soil ‘sample’ preparation procedure, domain dimensions and simulation control algorithms. A summary of the key parameters adopted for the DEM simulations is given in Table 1.

Sand grains were modelled as spherical particles with a particle size distribution (PSD) approximating that of Leighton Buzzard (LB) sand, grade 14/25 (median diameter, D₅₀ = 0·8 mm). The sand was modelled as dry, therefore fully drained conditions were assumed at all times, with no excess pore pressure development. This assumption is consistent with the laboratory tests that were used for validation. The PSD of LB sand grade 14/25 was used for an earlier DEM study (Macaro et al., 2014) to calibrate the particle contact parameters by way of simulations of triaxial tests performed by Schnaid (1990) on the same soil. For reasons of computational efficiency, the pipe loading simulations described in this paper were performed using scaling factors of 2 and 4 applied to the original PSD. Some of these DEM simulations were then compared against 1g model tests performed by Sandford (2012) using a 50 mm diameter pipe segment on LB sand grade DA30 (D₅₀ = 0·46 mm). Fig. 1 shows the experimental and numerical grading curves.

For the mechanical contact model in the DEM analyses, the Hertz–Mindlin no micro-slip solution (Mindlin, 1949; Cundall, 1988) was adopted for the normal and shear components of the contact forces. Moreover, to account for the non-sphericity of real sand grains, an ad hoc rolling resistance model was incorporated in the YADE code. This was based on the models of Iwashita & Oda (1998) and Plassiard et al. (2009), but modified to be consistent with the Hertz–Mindlin formulation (Macaro, 2015). The calibration of the parameters of the DEM contact model (see Table 1) entailed a series of 3D triaxial test simulations using a periodic cubic cell, comparing the DEM simulations with the results of experimental triaxial tests performed by Schnaid (1990) on LB sand grade 14/25. Further details of the calibration procedure can be found in Macaro et al. (2014), where it is also shown that the inclusion of a moment–relative rotation contact law is essential for replicating vertical and lateral pipe–soil interaction at small displacements to a good degree of accuracy.

The pipe segment was implemented in the YADE code as a circular cylinder of finite stiffness and diameter D = 50 mm. This diameter, which is 4 to 20 times smaller than typical pipeline diameters used offshore, was chosen in order to replicate the solid steel model pipe segment used in the laboratory tests of Sandford (2012). As shown in Table 1, the mechanical parameters for the pipe were based on the elastic properties of steel.

With regard to the interface friction between pipe and sand, relevant literature (Potyondy, 1961; Kishida & Uesugi, 1987; Randolph et al., 1994; Jardine et al., 2005; Lings & Dietz, 2005; Ho et al., 2011) shows that this friction depends on the level of dilation (for instance in the case of dense sand a peak angle different from the constant volume angle is exhibited), particle size to pipe roughness ratio (Kishida & Uesugi, 1987), level of confinement (Randolph et al., 1994; Ho et al., 2011) and any particle breakage (Ho et al., 2011). In these experimental works the measured interface friction is related to the internal friction angle or angle of shearing resistance exhibited by the sand in bulk, which is much easier to measure than the interparticle (grain to grain) friction angle. Here instead, the interface friction angle was assumed to be a fraction (2/3) of the interparticle friction angle to be consistent with the DEM approach. In fact the angle of shearing resistance is an average property of the bulk material which depends not only on the interparticle friction but also on (non-spherical) grain shapes, whose effect is (partly) captured by the rolling resistance contact model adopted. The bulk interface friction angle, like the bulk internal friction angle, depends on both interparticle friction and particle shape and varies with the amount of dilation undergone by the sand, so that in the literature at least two different angles, peak and ultimate state (constant volume) are measured by experiments such as the direct shear test, simple shear test, ring torsion and so on. The interface friction angle assigned here expresses the friction experienced by each sand grain in contact with the pipe, rather than a macroscopic soil–pipe friction angle averaged over the pipe. Since no experimental measure of pipe roughness from Sandford (2012) is available, 2/3 of interparticle friction falls in between the bounds experimentally measured for smooth and rough steel/sand interfaces. Also, from the numerical simulations of Sandford (2012) it emerges that for semi-rough to rough pipes, the pipe–soil behaviour during lateral displacements is not highly sensitive to the chosen interface friction angle.

A 3D prismatic particle domain was used, with boundary conditions to model a thin slice of soil perpendicular to the pipeline axis (the dimensions of this slice are discussed below). The sample preparation procedure consisted of: (a) random generation of non-contacting particles in the space enclosed by frictionless walls; (b) gradual application of gravity for particle settlement; and (c) generation of the pipe with its invert a small height above the settled particle domain. The homogeneity of the domain at the end of deposition was checked by assessing the average porosity within ten horizontal layers: no significant variation with depth was found (Macaro, 2015). The distribution of vertical stress on the floor of the container was also investigated. This was essentially constant, suggesting that there was no influence of the walls during deposition of the particles. It also confirmed the homogeneity of the sample in the horizontal direction.

Numerical specimens with different relative densities were obtained by assigning a fictitious value of the interparticle friction angle during deposition, with high and low porosities generated using friction angles of 70° and 0·1°, respectively. The particles with the 70° friction angle led to the loosest specimen, with a final porosity of 0·497 (taken as the minimum relative density, RD = 0%). When the 0·1° friction angle was used, the particles settled in their densest configuration, with a final porosity of 0·390 (taken as the maximum relative density, RD = 100%). Then, after the particles had settled in static equilibrium, the value of interparticle friction calibrated from the triaxial tests (26°) was assigned and the pipe was loaded. The final porosity at the end of the deposition phase employed for the present analyses was 0·480, corresponding to a relative density (for the numerical specimen) of 18·6%. This is close to the relative density of the sand specimens employed in the experiments of Sandford (2012), where a mean value of 16% (based on standard laboratory tests for minimum and maximum density) is reported at a depth of 20 mm (0·4D) from the soil surface (see Figure 5.6 in Sandford (2012)). The computational times required for sample generation varied between 10 and 35 h depending on the number of cores employed and the particle mechanical parameters.

Preliminary DEM simulations involving monotonic pipe displacements in the vertical and lateral directions were performed with various particle domain thicknesses and particle sizes to identify appropriate values for these parameters (Macaro, 2015). The use of periodic boundaries in the out-of-plane direction – that is, along the pipe axis – was attempted since it would be best to minimise the thickness of the domain. However, the use of periodic boundaries in the YADE code is incompatible with the pipe object and therefore smooth rigid walls in the out-of-plane direction had to be employed instead.

As shown in Fig. 2, the domain thickness adopted in the simulations presented here is 5D₅₀. This value was chosen after a parametric study using thicknesses of 3, 5 and 20 times D₅₀, and corresponds to the minimum required for the results to be independent of the domain thickness – that is, ensuring a negligible influence of the frictionless walls in the out-of-plane direction. Table 2 lists the dimensions, particle counts and other properties of the DEM samples in Fig. 2. In this and subsequent figures the particles are shaded in layers according to their vertical positions at the start of the simulation. The initial height of each layer was taken as D/10 (5 mm) for the cyclic sideswipe analysis (Fig. 2(b)) and D/4 (12·5 mm) for the other analyses.

The scaling factor that was applied to the PSD of the reference material, LB sand grade 14/25 (see Fig. 1), depended on the dimensions of the particle domain in question. Preliminary simulations (Macaro, 2015) assessed the influence of PSD scaling factors of 1, 2 and 4 on the pipe–soil response. A scaling factor of 2 (giving D/D₅₀ = 31·2) was used for the cyclic sideswipe analysis (Fig. 2(b)), and no significant difference with respect to the unscaled PSD was observed. Because of the small particle size, the computational time for this simulation was 7 weeks using eight cores of a 16-core workstation with 32 GB of RAM. Such a small PSD would have been impractical for the deeper domains, giving even longer run times. A PSD scaling factor of 4 was therefore used for the monotonic sideswipe and cyclic constant V analyses. With this larger PSD, computational times were reduced to 4 weeks for the monotonic sideswipe test and 3 weeks for the one-cycle constant V test that was used for validation. The disadvantages of a larger PSD were (a) a limited but tolerable difference in pipe–soil response, and (b) less detailed observations of the kinematics in small areas of interest – for example, within the soil berm. To the authors’ knowledge, the D/D₅₀ ratios adopted here (31·2 and 15·6) are still the largest used so far in DEM simulations of buried or unburied pipelines in sand; for comparison see Calvetti et al. (2004) and Yimsiri & Soga (2006).

Displacement control was used to perform DEM simulations of both monotonic and cyclic sideswipe tests in which the pipe was subjected to lateral displacement, u, at constant embedment, w (Fig. 3). These analyses were carried out by applying constant lateral displacement rates, $u˙$⁠, of 5 mm/s (monotonic test) and ± 10 mm/s (cyclic test) to the pipe, while maintaining its vertical position fixed.

Load control was required for the DEM simulations involving lateral pipe displacements at constant vertical load (Fig. 3). As shown in the experimental work of Sandford (2012), maintaining a constant vertical load on the pipe throughout large lateral displacement cycles can be challenging owing to the continuous changes of geometry of the surface of the sand bed. For the DEM simulations at constant V that are presented here, a servo control algorithm was employed.

Contact viscous damping was used in all of the DEM simulations. A viscous contact damping ratio of 0·2 was used for both the normal and shear contact directions. During preliminary studies (Macaro, 2015) it was observed that damping mainly influences the preparation of the numerical sample, by reducing the rearrangement of particles during gravitational settling – that is, the higher the damping, the looser the sample obtained at the end of the generation procedure. During the simulations involving movement of the pipe, for the same sample, the only noticeable effect of increased damping is to attenuate the fluctuations of the resultant vertical and lateral forces on the pipe.

This section describes a displacement-controlled analysis in which the pipe was subjected to a monotonic lateral displacement at constant embedment. Of particular interest is the pipe response just after the onset of lateral movement. The DEM predictions are compared with experimental results obtained by Sandford (2012).

The full analysis consisted of three stages. In the first stage, the pipe was pushed vertically into the soil at a rate of 5 mm/s until an invert embedment w = 0·24D was reached. Then, over a further penetration of 0·01D, the vertical displacement rate was gradually decreased to 0·5 mm/s to reduce fluctuations in the resistance. In the second stage, the vertical displacement was kept fixed (Δw = 0), while the pipe was moved laterally at a rate of 0·5 mm/s until a horizontal displacement u = 0·1D was reached. This small displacement rate allowed observation of the rapid changes that occur in the vertical and horizontal loads immediately after lateral breakout begins (Zhang et al., 2002a; Sandford, 2012). In this stage the data acquisition frequency was initially 10⁶ Hz. This was gradually decreased to 10³ Hz when fewer data points were required. In the third stage, the horizontal displacement rate was increased to the usual value of 5 mm/s until a final horizontal displacement u = 0·9D was reached. The run time using an eight-core workstation with 32 GB of RAM was about 4 weeks.

The purpose of this ‘sideswipe test’ simulation was to define the yield surface, or bearing capacity envelope, in the vertical and horizontal (V, H) load space. Experiments on shallow foundations under combined (V, H) loading have shown that the load path during a sideswipe test tracks through a series of yield surfaces (Tan, 1990). This is because as V reduces, downward plastic displacements develop to balance upward elastic displacements (to maintain Δw = 0), and the downward plastic displacements cause hardening or expansion of the yield surface (Tan, 1990; Martin & Houlsby, 2000, 2001). In typical pipe–soil interaction problems, the virgin penetration (plastic) stiffness is much smaller than the vertical unload–reload (elastic) stiffness. Hence the amount of yield surface hardening during a sideswipe is negligible and in practice it may be assumed that the load path followed in a sideswipe test closely approximates the current yield surface. Gottardi et al. (1999) discuss the correction to be applied to the sideswipe (V, H) loads to account for finite soil elasticity. For the present analysis, the difference between the corrected and uncorrected yield surfaces is negligible.

Figure 4 shows the load–displacement response during the vertical penetration stage. The embedment of the pipe invert relative to the undisturbed soil surface level, w, is normalised by the pipe diameter, D. The vertical force on the pipe per unit length, V/L, is normalised by γ′D², where γ′ is the effective unit weight of the soil. For the DEM sample γ′ = G_s(1 − n)γ_w = 13·5 kN/m³ where G_s is the specific gravity of the particles and n is the porosity of the sample after deposition (see Tables 1 and 2). The predicted load–displacement curve in Fig. 4 is typical for a pipe on sand, with the bearing capacity increasing approximately linearly with embedment (Zhang, 2001). It emerges that the DEM prediction and the experimental curve are quite close, especially for embedments up to 0·15D. At larger embedments the predicted resistance starts to fall below the measured resistance.

The pipe response at the onset of lateral movement is shown in Fig. 5(a) while the same plot with the x-axis magnified 10 times is shown in Fig. 5(b). As the pipe is first moved laterally, the horizontal load rises to a peak and then drops, while the vertical load decreases and remains approximately constant when the remaining horizontal load is mobilised. The peak horizontal loads in the DEM simulation and the experimental test, being 0·11V₀ and 0·12V₀, respectively, exhibit close agreement. The peak horizontal load is reached almost immediately (u = 0·00015D) in the DEM simulation, while a larger displacement is necessary to mobilise the maximum horizontal resistance in the experiment (u = 0·0038D).

In Fig. 5(c) the DEM prediction and the experimental result for the sideswipe load path are plotted in normalised (V, H) space. A reasonably good overall agreement between the shape of the curves can be observed in the initial part of the test, where the load path tracks along the (V, H) yield surface. However, the spacing of the data markers confirms that, as already noted, the lateral displacement required to traverse the yield surface is much larger in the experiment than in the DEM simulation. In Fig. 5(d) the loads V and H are normalised by the bearing capacity under pure vertical load, V₀. Note that, due to high fluctuations exhibited in the DEM simulations, the V₀ adopted for the DEM results was determined from a straight trend-line; see Fig. 4. The initial difference observed at V/V₀ ≈ 1 is because in the experiment there was some relaxation of vertical load while the pipe was held in position prior to the start of lateral movement.

The results of the experimental sideswipe test and its DEM simulation show that, after an initial peak in H while V is decreasing, both load components decrease and approach steady residual values (Fig. 5(c)). This corresponds to a point on the yield surface referred to as the parallel point (Tan, 1990; Martin, 1994; Zhang, 2001). Martin & Houlsby (2000) discuss the analogy between macro-element models for shallow foundations and strain-hardening plasticity models of the type used in critical state soil mechanics. The parallel point plays the same role as the critical state, and coincides with the peak of the plastic potential function in (V, H) space. In the context of transverse pipe–soil interaction, the parallel point marks the transition between pipe heave occurring during lateral displacement at low vertical loads (zone of strain softening) and pipe settlement occurring during lateral displacement at high vertical loads (zone of strain hardening). In this study it is clear from Fig. 5(c) that the parallel point is not located at the peak of the yield surface, indicating that the flow rule is not associated. This is in agreement with the findings from previous laboratory investigations of shallow foundations on sand (Tan, 1990; Gottardi et al., 1999) and partially embedded pipes on sand (Zhang, 2001; Sandford, 2012).

Figure 6 shows the two horizontal load–displacement curves (numerical and experimental) as the sideswipe is continued to larger lateral displacements, up to u = 1D. Here the comparison between the DEM prediction and the experimental data from Sandford (2012) is less satisfactory, with the analysis overpredicting the measured lateral resistance by a factor of about 3 on average, and showing considerable fluctuations. The difference could be due to a larger berm forming in the DEM simulation than in the experiment, although direct comparison is not possible because Sandford (2012) did not record the evolution of the surface profile for this test. It is also possible that more accurate DEM predictions could be obtained by employing a smaller scaling factor for the numerical PSD.

An important outcome from this validation analysis is that the DEM model is inherently able to capture the non-associativity of the sand behaviour and its manifestation as a non-associated flow rule for the macroscopic pipe–soil system. Although the evaluation of a plastic potential function goes beyond the scope of this work, further DEM analyses could be used to perform load- or displacement-controlled probes aimed at examining the orientation of the incremental plastic displacement vectors at different points on the yield surface. This could inform the development of an analytical formulation of the plastic potential function for a macro-element plasticity model.

As an initial test of the ability of the DEM approach to handle cyclic lateral loading, an analysis involving several lateral sweeps at constant embedment, w, was performed up to a final lateral displacement exceeding u = 3D. For this simulation, the PSD of the reference material (LB sand grade 14/25) was scaled up by a factor of only 2. This is sufficiently small to be representative of the actual grain size (Macaro, 2015).

The simulation was performed under displacement control, subjecting the pipe to the prescribed (w, u) path illustrated in Fig. 7(a). The usual lateral displacement rate of 10 mm/s was used. Each of the three sweeps at constant w was preceded by oblique penetration at an angle of 20° to the horizontal, imposing an incremental penetration of 0·1D. The aim of this simulation was to make a preliminary evaluation of the performance of the DEM model at large lateral pipe displacements, and to test features such as cyclic reversals of direction, active berm growth and dormant berm collection. The chosen pipe trajectory was designed for this purpose.

Figure 7(b) shows the normalised horizontal load–displacement response from the DEM analysis. Front views of the particle domain near the start and at the end of each lateral sweep are depicted in Fig. 8. The results illustrate the typical mechanical behaviour of a pipe cyclically swept over a sand surface. Initially, the horizontal load increases as the pipe is pushed obliquely into the sand (points A→B). As the test proceeds at constant w = 0·1D (B→C), the pipe scrapes across the surface and a berm of soil builds up in front of the pipe, gradually increasing the lateral resistance. After the first reversal, the previously ‘active’ berm is left ‘dormant’ (Cheuk et al., 2007). The lateral resistance returns to zero and starts gradually rising with the opposite sign (C→D). During the return sweep at constant w = 0·2D, another berm forms on what is now the leading side of the pipe, continually increasing the lateral resistance (D→E). Following the second reversal, the horizontal load changes sign and builds up rapidly during the oblique leg (E→F). When the pipe is moved for a further 2D at constant w = 0·3D (F→G), the lateral resistance is higher than in the first outward sweep (B→C), particularly after the active berm has merged with the dormant berm that was created by the reversal at C. Towards the end of the simulation, the large oscillations observed in the horizontal load are due to the proximity of the domain boundary on the right.

Figure 9 shows velocity vectors for particles lying within a central section of the domain 1·1D₅₀ thick (the entire domain is illustrated in Fig. 2(b)). The vectors are determined by considering the change of displacement occurring over the time the pipe moves 0·1D. A large region of the soil domain is affected by the initial oblique movement of the pipe (Fig. 9(a)). This indicates that the chosen depth of the DEM sample was not quite sufficient, even for the small amount of vertical displacement to which the pipe was subjected. The very large vectors in Figs 9(d), 9(e), 9(h) and 9(l) are believed not to be the result of fluctuations or numerical instabilities, but the expected high velocities of particles falling from the top of the berm. During berm collection (Fig. 9(j)), the particles of the active berm initially slide over those of the dormant berm. Then they all move together, sliding over the soil below, which remains practically stationary (Fig. 9(l)). At larger displacements, such as in Fig. 9(m), the particles in the central part of the berm move together, with the particles at the top and bottom of the berm moving more slowly. This is because the middle particles move with the pipe, while the lower particles are constrained by the stationary soil beneath.

The particle position and velocity plots during large pipe displacements demonstrate that the proposed DEM modelling approach is able to capture the large plastic deformations and complex free surface evolutions taking place in the sand. It is worth noting that these phenomena cannot readily be captured by continuum-based finite-element or finite-difference analyses, as typically used for pipeline modelling, without regular remeshing and careful attention to constitutive modelling. In conclusion, some key aspects of the pipe–soil interaction successfully modelled in this initial cyclic DEM simulation were: (a) deposition of a dormant berm following the first cyclic reversal (point D in Figs 7–9); (b) collection of the dormant berm during the second outward sweep (points F→G); and (c) large lateral displacements, up to several pipe diameters (points C and G).

A DEM simulation involving one full cycle of lateral displacement at constant vertical load was conducted to replicate more realistically the loading conditions experienced by a laterally buckling pipeline in the field. Experimental results obtained by Sandford (2012) were used for validating the numerical predictions. In the DEM analysis the pipe was first pushed vertically into the soil a rate of 10 mm/s until an invert embedment w = 0·17D (8·5 mm) was reached. The bearing capacity at this point was recorded as V₀. The pipe was then retracted at a rate of 1 mm/s until the vertical force, V, reduced to give an overloading ratio R = V₀/V = 10. Using the servo control algorithm described previously, the embedment w was adjusted to keep the vertical load constant (ΔV = 0) while the pipe was moved laterally at $u˙$  = 10 mm/s until a horizontal displacement u = 2D (100 mm) was reached. The pipe was then moved back to its initial lateral position, u = 0, again with w being varied as required to maintain ΔV = 0.

The loading sequence in the DEM simulation was identical to that in the experiment of Sandford (2012), although the pipe displacement rates were different. In the experiment the rates were 0·1 mm/s for vertical penetration, 0·001 mm/s for vertical unloading and 0·05 mm/s for the lateral sweeps. In addition, the experiment used a smaller displacement rate of 0·001 mm/s during the early part of each lateral sweep, to assist the stable operation of the vertical load control system.

Figure 10(a) shows the trajectory of the pipe invert in the DEM simulation and in the experimental test, while Fig. 10(b) compares the two normalised horizontal load–displacement curves. The DEM prediction shows that during the outward lateral sweep, the pipe embedment w undergoes small undulations but remains approximately constant overall. The horizontal force H increases gradually as the pipe builds up a berm of soil while scraping the surface. When the direction of movement is reversed, the pipe experiences a rapid increase of penetration, and the horizontal force starts increasing in the opposite direction. In the latter part of the return sweep, the pipe moves slightly upwards to keep the vertical load constant.

Comparison between the numerical predictions and experimental results is provided in Fig. 10. It is important to emphasise that in the DEM analysis, no additional tuning of the particle and pipe micromechanical properties was performed after the initial calibration based on triaxial test simulations since an authentic prediction is attempted here. The pipe embedment (Fig. 10(a)) at the end of the outward sweep is 8·5 mm in the simulation and 7·5 mm in the experiment. At the end of the return sweep, the embedment is 12·7 mm in the simulation and 12·1 mm in the experiment. However, the predicted lateral resistance during the outward sweep is generally larger than the measured value by a factor of 2 or more (Fig. 10(b)). The agreement between the predicted and measured forces is somewhat better on the return sweep. This suggests that a possible improvement to the current DEM model would be to increase the width of the sample (at its extreme lateral position the pipe is only 1D away from the domain boundary).

Figures 11 and 12 show further comparative results from the DEM analysis and the experimental test, in terms of the soil profile and the particle velocities. These figures reveal the failure mechanism in the soil at various stages of the pipe movement, with a good match between simulation and experiment. As the pipe starts its lateral motion, the sand particles initially move towards the free surface (snapshot A). The particles displaced by the pipe then follow a predominantly horizontal trajectory and form a berm (B and C). When the pipe reverses (D), the berm is left dormant and a new berm is built up on the other side of the pipe. Some particles fall from the inside edge of the dormant berm, giving rise to large velocity vectors in that area. In the remaining snapshots (E and F), the active berm merges with the small dormant berm that was created during the initial vertical penetration stage.

Three cyclic lateral displacement analyses at constant vertical load were performed for the same initial embedment but various overloading ratios, R. The purpose of this set of simulations was to investigate the effect of the vertical loading history on the subsequent trajectory and horizontal load–displacement response of the pipe. Because the DEM contact stiffness values have little influence on the predicted pipe–soil response (see Figure 7.28 in Macaro (2015)), these properties were adjusted to reduce the run times for this set of analyses. The Young's moduli for the sand particles and the pipe segment were set to 70 MPa and 200 MPa, respectively, three orders of magnitude smaller than the values in Table 1 that were used for the previous analyses in this paper. No other properties were altered.

The simulations consisted of the following loading stages.

• (a)

  Vertical penetration at a rate of 10 mm/s until an invert embedment w = 0·17D was reached; the bearing capacity at this point was recorded as V₀.

• (b)

  Vertical unloading at a rate of 1 mm/s until a specified overloading ratio R = V₀/V was reached. Analyses were performed for R = 5, 10 and 20 to study the subsequent lateral breakout behaviour of lightly, moderately and heavily overloaded pipes.

The vertical load was then kept constant (ΔV = 0) while the pipe was moved laterally at a rate of 10 mm/s through the following sequence of movements.

• (c)

  Out to u = 2D.

• (d)

  Back to u = 0.

• (e)

  Out to u = 2D.

• (f)

  Back to u = 0·5D.

• (g)

  Out to u = 2·5D.

Figure 13(a) shows the normalised vertical load plotted against the normalised embedment. This confirms that the vertical load remains nearly constant throughout the cyclic lateral displacement phase of each analysis.

Figures 13(b) and 13(c) show the pipe invert trajectories and the normalised horizontal load–displacement curves for the three simulations. Corresponding particle position plots are depicted in Fig. 14. These DEM simulations show that during breakout – the early part of the first lateral sweep – the lightly overloaded pipe (R = 5) moves downward, the heavily overloaded pipe (R = 20) moves upward and the moderately overloaded pipe (R = 10) maintains nearly constant embedment. Because of the limited domain size used for these analyses, detailed quantitative comparisons with experimental data cannot be carried out. There is, however, a broad similarity in behaviour with respect to the experimental tests performed by Zhang et al. (2001) and Sandford (2012). These researchers found that when a vertically overloaded pipe on sand was moved laterally, the transition between initially downward and initially upward pipe motion occurred at a critical overloading ratio around R = 10. The numerical predictions presented here are able to replicate such behaviour.

At large lateral displacements, and during cyclic reversals, the trajectory of the pipe and the horizontal load–displacement response are governed not by the initial vertical overloading ratio, but by the weight of the pipe (i.e. the vertical load V = V₀/R that is plotted in normalised form in Fig. 13(a)). Figs 13(b), 13(c) and 14 show that the heaviest pipe (the one with R = 5) experiences a large increase in embedment during the lateral cycles, with upward movement only occurring at the ends of the second and third cycles. As it moves downwards, this pipe builds up a larger berm than the lighter pipes, although the berm formation and the resulting upward movement of the pipe, as well as the large resistance developed in the final outward sweep, are clearly influenced by the proximity of the domain boundary on the right. The medium-weight pipe (R = 10) remains at approximately the same vertical level during the initial part of the first outward sweep, before beginning to rise. In each subsequent sweep it undergoes penetration followed by upward movement. The lightest pipe (R = 20) follows an upward trajectory for most of the first outward sweep. Thereafter it experiences limited penetration in the first half of each sweep, and scrapes only a thin layer of sand to create a shallow trench and small berms. As expected, the soil resistance mobilised during the cyclic lateral loading of this pipe is the smallest of the three cases considered.

Numerical analyses using the 3D DEM have been used to study pipe–soil interaction during large lateral displacements (both monotonic and cyclic) of a pipe segment partially embedded in sand. The DEM contact law parameters were taken from a previous calibration based on triaxial test simulations, and were not varied during the analyses presented here. Therefore the simulations can be considered predictions of the pipe–soil response obtained by the DEM. Where possible, these predictions were compared with results from experimental model tests to assess the extent to which the predictions capture (or fail to capture) the experimentally observed behaviour.

The early part of a monotonic sideswipe analysis at constant w gave a predicted (V, H) load path that matched well with experimental data. The location of the parallel point was identified, again showing good agreement between prediction and experiment. For larger sideswipe displacements to around 1D the comparison was not as good, with the predicted lateral resistance being higher than the measured resistance.

A cyclic sideswipe analysis was used to model several horizontal sweeps at progressively increasing w, reaching lateral displacements of more than 3D. A smaller particle size allowed understanding of the mechanisms occurring at a particle level. Specifically, the formation and growth of a soil berm as the pipe moved laterally, and its deposition and subsequent collection after one cycle, were able to be captured directly by the DEM simulations. Large displacements of the soil and complex changes in the surface profile were also captured directly.

A cyclic lateral displacement analysis at constant V was used to perform further validation of the DEM approach, replicating more realistically the conditions experienced by a laterally buckling pipeline in the field. The soil failure mechanisms at various stages of the cyclic loading history were revealed from particle position and particle velocity plots, and these again showed good agreement with experimental observations.

Three further constant V analyses investigated the influence of the initial vertical overloading ratio on the pipe trajectory during lateral breakout. Although detailed quantitative comparisons with experimental data were not possible, the predictions of the DEM model were broadly consistent with observations from available experimental studies: lightly overloaded pipes follow a downwards trajectory during the initial lateral displacement, whereas heavily overloaded pipes experience an upwards trajectory. During large cyclic lateral displacements, the heaviest pipe formed a deep and steep-sided trench with large berms, whereas the lightest pipe created a shallow trench with small berms.

The use of the DEM for this problem is attractive because it allows investigation of the onset of localisation and the formation of geometrically complex failure mechanisms occurring in the soil (including cyclic pipe reversals leading to berm formation) without requiring any modification of the basic numerical procedure. The main limitation of the method, however, is its currently high computational cost, which is likely to confine its use to the research community for some time ahead.

Many of the simulations reported in this article were run on the high-performance computer managed by the Advanced Research Computing facility at the University of Oxford. The first author would also like to acknowledge the Engineering and Physical Sciences Research Council (EPSRC), EU H2020 RISE Hercules grant (778360) and BP for the financial support received.