Decoupled Analysis of a Multi-Layer Flexible Pipeline Buried in Clay Subjected to Large Lateral Soil Displacement

Abstract

Multilayered flexible subsea pipelines may experience significant lateral movements due to manmade and environmental geohazards. These pipelines incorporate several structural and protective layers to resist different loads, and may require additional protection such as trenching, rock placement, or burial. In practice, simplifications are considered due to the complexities and uncertainties involved in the multi-layer pipe structure and the surrounding soil, compromising the pipe structure or the soil behavior. These simplifications are applied either on the pipe by assuming a rigid section or on the soil by representing it as elastic springs, which may result in inaccuracies. This study proposes a decoupled methodology combining the Coupled Eulerian–Lagrangian (CEL) model for soil displacement with a small-strain finite element analysis of the flexible pipe. This approach aims to accurately capture cross-sectional deformations and local stresses due to soil movement while maintaining reasonable computational effort. A parametric analysis was conducted to assess the impact of several variables on failure risk. The deformed cross-section was then used for a collapse analysis to determine critical loads at maximum operational depth. The study showed that modeling parameters such as soil strength and internal diameter might significantly influence pipe failure and the risk of collapse.

1. Introduction

Multi-layer flexible pipelines are extensively used in the offshore oil and gas industry. While these pipelines may be buried in the soil or trenched for physical protection [1], they remain vulnerable to significant ground deformations due to landslides, fault crossing, ice gouging in Arctic regions (Figure 1a), and other natural events. When subjected to large soil displacements, flexible pipes are also prone to large deformations. Failures can occur under bending situations, leading to extreme tension and compression within their structure, putting the pipeline’s integrity at risk. In contrast to rigid pipes, multi-layer flexible pipes have a low bending stiffness, allowing them to deform considerably and gain strength from the surrounding remolded soil and adjacent undisturbed soil. Thus, it is essential to incorporate the lateral nonlinear soil interaction effects into the performance analysis of buried multilayer flexible pipes subjected to significant lateral soil displacements. However, a coupled analysis incorporating the complexities of both the flexible pipe structure, involving contacts between several different layers, and the advanced constitutive soil models make the analysis computationally costly and inefficient. Therefore, in practice, depending on the area of concern, the accurate modeling of one of these two elements, i.e., the pipe or soil, cannot be achieved at a reasonable computational cost. The current study addressed this practical need by developing a decoupled approach to incorporate the uncertainties of both elements at the same time, while maintaining a reasonable computational effort.

Theoretical and experimental models exist to predict the ultimate lateral resistance or force–displacement responses of embedded moving structures [2,3,4,5,6,7,8], with some specifically concentrating on the interactions between trenched/backfilled pipelines and the surrounding soil [9,10,11,12,13]. All these research works assume the pipeline is rigid, with the emphasis placed on the soil properties and the extent of soil resistance against the moving pipe.

The pipe’s cross-sectional deformation was considered in the numerical models presented in a few studies [14,15,16]. The majority of these studies focus on onshore water transportation systems, mainly concrete pipes. In [14], the concrete pipeline is considered static and a continuum FE model is used to examine the impact of external pressure on the pipeline section. The behavior of fiber-reinforced pipelines while under the effect of surface loads was studied in [15]. Although the pipelines considered are onshore water transportation pipes, the soil modeling and the centrifuge tests used to validate the models are quite like the ones used in the offshore industry. Regarding composite pipes, more specifically, very flexible pipes, a numerical study can be found in [16]. The authors developed a more sophisticated FE model using ABAQUS. However, the pipeline was modeled assuming the layers would be perfectly bonded together. The soil itself was modeled using the native ABAQUS soil elements, and the study was able to determine the section deformation and make comparisons between steel and polymeric pipelines.

An FE continuum model was developed by [17] to simulate the impact of fault crossings on a buried pipeline. The authors considered the cross-sectional pipe deformation using the ABAQUS submodeling technique. The study determined the local buckling deformations in the area subjected to bending.

An elastic pipe section was used by [18] to analyze the wave loads’ effects on buried offshore pipelines. The authors developed a finite element model for porous soil media. An explicit continuum FE model was used in [19] to replicate ice gouging conditions. The sub-gouge soil deformations were translated into pressures and were directly exerted on the exterior of a buried subsea pipeline. This model used a single-layer shell deformable pipe, making it possible to simulate the ovalization of the cross-section. However, the authors focused further on the global response of the pipeline compared with a simplified Winkler beam-spring model. A continuum finite element model was also employed by [20] to simulate a buried pipeline exposed to the fault displacement caused by an earthquake. Both the pipeline and the seabed soil were represented as continuum Lagrangian domains in the model. Nevertheless, there was no discussion regarding section ovalization, given that the paper’s goal was to compare the continuum model with the Winkler models.

Continuum soil models with conventional Lagrangian meshes for large deformation suffer from mesh distortion and instability with large displacements. Thus, the current study employs the Coupled Eulerian–Lagrangian (CEL) method for a large deformation analysis of the seabed soil that allows for the material to flow through the mesh. Also, the pipe layers were set to be detachable to capture the local instabilities and separations.

Typically, the multi-layer flexible pipe assembly comprises multiple concentric layers. The operational requirements of each pipe determine the dimensions and quantities of these layers. Conversely, these loads directly result from the environment in which the pipeline is deployed. This study concentrates on the section of the tube that bears radial loads, referred to as the “structural nucleus” or “annulus” (see Figure 1b).

Figure 1. (a) Pipeline submitted to ice gouging (adapted from [21]); (b) a multi-layered flexible pipe and the layers that comprise the structural nucleus (adapted from [22]).

Figure 1. (a) Pipeline submitted to ice gouging (adapted from [21]); (b) a multi-layered flexible pipe and the layers that comprise the structural nucleus (adapted from [22]).

The CEL model for lateral pipeline–soil interaction analysis was developed in ABAQUS, and the structural model of the multi-layer flexible pipe was constructed in ANSYS. This configuration allowed for a cost-effective analysis of the interaction between flexible pipe and soil, incorporating both the soil and pipe models’ complexities. Further details of the analysis sequence will be presented in Section 2.

2. Proposed Methodology and Simulation Flowchart

Overall, the CEL model was developed to simulate the pressure distribution around the rigid pipe section during substantial lateral displacement in clay. Strain-softening effects were also incorporated by programming a user subroutine in ABAQUS. Then, the obtained soil pressure was fed into the flexible pipe model in ANSYS to simulate the structural response of the pipe and the failure analysis.

As per API and DNV regulations [23,24], deformation of the cross-section is not inherently problematic, but it significantly raises the risk of failure due to collapse. Therefore, the suggested approach involves determining the cross-section deformation using soil loads obtained from a CEL model and then employing this deformed cross-section in collapse analyses. Further details of the sequential analysis are illustrated in the flowchart shown in Figure 2, which utilizes four different numerical models:

• Coupled Eulerian–Lagrangian (CEL) pipeline–soil interaction model;
• Transient structural cross-section deformation analysis;
• Wet collapse quasi-static structural analysis;
• Dry collapse quasi-static structural analysis.

The first step is computationally more expensive: the CEL model employs a hybrid mesh and advanced soil formulations to calculate the acting pressure on a pipe section resulting from the pipeline’s interaction with the seabed. A simple representation of the flexible pipe is used at this stage as a rigid pipe section, assuming that the cross-section deformation is too small to influence the soil resistance. During this analysis, the strain-softening effects are incorporated into the soil by programming a user subroutine in ABAQUS.

Once the pressure profiles are obtained from the CEL analysis, the output is fed into the ANSYS as pressure loads on the flexible pipe model. An advanced pipe cross-section analysis is then conducted by ANSYS, allowing for the separation of the pipe layers and by incorporation of elastic–plastic materials and contacts. The equivalent thickness is a consecrated simplification for radial load analysis in flexible pipes [25,26,27]. The typical transient analysis timeframe is about twenty minutes, considerably reducing the analysis time and making it extremely practical for design purposes. This simplification, associated with the separation between the cross-section analysis and the pipe–soil interaction, allowed for several flexible pipe configurations to be tested with different materials or internal layer geometries. Also, aside from the external diameter, every parameter of the flexible pipe may be altered and quickly simulated.

The transient structural analysis provides the deformed cross-sections during the time frame of the soil interaction process. It is possible to calculate the cross-section ovalization with radial deformations, which can be defined according to two different regulations.

As stated in API 17B [28], Equation (1) can be used to define pipe ovalization. The out-of-roundness of the pipe’s cross-section is an essential factor in collapse due to external pressure. The higher the ovalization, the higher the pipe failure risk [27]. Therefore, limits of 0.2% and 3% ovalization were set in the regulation for the unloaded loaded conditions, respectively. This criterion is employed in the crushing analysis during the flexible pipe’s launching stage as well.

$${API}_{oval}={\displaystyle \frac{D_{max}-D_{min}}{D_{max}+D_{min}}}$$

(1)

In Equation (1), $D_{max}$ represents the cross-section’s largest external diameter, and $D_{min}$ is the smallest. The same regulation also notes that ovalization might be a consequence of the excessive bending of a flowline. The section dedicated to bending tests of tensioned lines specifies that the allowed variations in external diameter (as a consequence of bending) are 3% for the loaded condition and 1% for the unloaded condition. One may notice that the word “ovalization” is not mentioned here; instead, the term direct diameter variation is used. Nevertheless, this will end up causing ovalization.

According to the DNV regulations [24], however, the bending is not the cause of failure itself, but the cross-section deformation caused by it. The regulation defines ovalization, calling it “ovality” (see Equation (2)). It is noticeable that a diameter variation, as per DNV’s definition, can result in greater ovalization compared to the API’s definition. Consequently, the DNV regulation is more stringent.

$${DNV}_{oval}={\displaystyle \frac{D_{max}-D_{min}}{D}}$$

(2)

In Equation (2), $D_{max}$ and $D_{min}$ remain consistent with API’s equation, while $D$ represents the nominal diameter of the cross-section. Regarding failure criteria, the initial ovality for fabrication (unloaded condition) must not exceed 0.5%. The ovalization due to bending should be kept under 3% to prevent any buckling. It is worth mentioning that the DNV regulation mainly applies to steel pipelines; nevertheless, it represents a similar concept to API’s equation and it will be used in this work as a comparatively conservative alternative.

Regardless of the definition, the ovalization obtained by the transient analysis may become an input for a collapse analysis. The model employed at this stage was the symmetric ring model, proposed by [29]. The model is quasi-static and uses the Arc-Length (Riks) method for the simulation, enabling the capture of buckling effects such as snap-through and determining the critical loading of the section with more accuracy than linear buckling techniques. This study conducted the analysis until ovalization was obtained, and the results were verified against the regulation limits.

3. Developed Numerical Models

3.1. CEL Model for Lateral Pipeline–Soil Interaction Analysis

A Coupled Eulerian–Lagrangian (CEL) model was created in ABAQUS to investigate the soil resistance mobilized against the pipe section during the large lateral displacement. The numerical model comprises a rigid pipe section embedded in a uniform clay seabed under plane strain conditions. Figure 3 depicts the model configuration and boundary conditions. The coupled Eulerian–Lagrangian model takes advantage of both Eulerian and Lagrangian meshing techniques. The Eulerian mesh is highly efficient in dealing with excessive deformation where there is discontinuity in the material surface across the neighboring elements, while the Lagrangian mesh is good at solving structural problems.

Eight-node linear Eulerian brick elements (EC3D8R) featuring reduced integration and hourglass control were used to represent the soil domain. Although unlikely, hourglass may occur in Eulerian meshes. Besides the hourglass control, the extreme refinement used for the soil mesh compensates for the reduced number of Gauss points in the elements. The contour plots also show no sign of discontinuities or elements with zero-energy modes.

For the material assignment in the soil domain, a discrete field was created to define the space where the soil material should be located, and the Volume Fraction Tool (VFT) was used to assign the material. The mesh size was tested to ensure that the thickness of the propagated shear bands contains at least three layers of elements in order to effectively simulate the failure mechanism, which requires a mesh size of approximately less than 1/30 of the pipe diameter in this research.

A general-purpose linear brick element with reduced integration (C3D8R) was utilized to model the pipeline in the Lagrangian domain, with the pipe section being constrained to a reference point located halfway through of the section to create a rigid body. In that case, the rigid pipeline has a basic carbon steel placeholder material, with 200 GPa of elastic modulus and a Poisson’s ratio of 0.3. Due to the multi-point constraints impeding any relative deformation, there was no deformation or stress in the pipe during the CEL phase. The clay seabed was modeled as a Eulerian domain, and the strain-softening effects were incorporated into the seabed model using a modified Tresca model, which is a simple constitutive model for clays in relation to the corresponding models for cohesive soils [30,31,32,33], programmed into a VUSDFLD subroutine in ABAQUS (Section 1.2.24 of [34]). More details can be found in [35,36]. Within the plastic region, the undrained shear strength (s), taking into account the strain-softening effect, was determined using the following equation:

$$s=s_u+\left(s_{u,r}-s_u\right){\displaystyle \frac{{\gamma }^p}{{\gamma }_r^p}}$$

(3)

In Equation (3), $s_u$ indicates the maximum undrained shear strength, $s_{u,r}$ is the residual undrained shear strength ($s_{u,r}=s_u/S_t$, and $S_t$ represents the sensitivity of soil), ${\gamma }^p$ denotes the accumulated plastic shear strain, and ${\gamma }_r^p$ represents the ${\gamma }^p$ corresponding to the reduction in shear strength from its maximum value to its residual value. The sensitivity of the seabed soil was taken as 1.3, consistent with the research assumptions in [8]. The seabed soil properties in [8] mainly derived from two centrifuge tests by [9] using a mixture of Speswhite Kaolin clay (50% of total weight) and Sil-Co-Silt (50%) and was pre-consolidated to 400 kPa to reach a target undrained shear strength of around 40 kPa.

A void domain above the seabed soil was also included to account for the geometrical variation in the seabed surface and the formation of heave (see Figure 3). To prevent material from escaping the domain, all external surfaces of the Eulerian domain were constrained in their normal directions. The pipeline was extended beyond the front and rear edges of the Eulerian domain to eliminate artificial friction between the soil and the pipeline’s end faces. General contact was applied throughout the entire CEL model for the Eulerian domain. However, the pipe–soil interface was specifically defined using Coulomb’s friction law, with the interface strength limited to half the soil’s undrained shear strength at the pipe’s springline, and an overall friction coefficient of 0.5 [8]. Separation is allowed during the contact between the pipe and seabed soil to avoid unrealistic soil resistance being induced during pipe displacement. The refined mesh used in the CEL model allows for a good convergence and captures the shear bands developed in the soil during significant deformations.

The developed CEL model was benchmarked using the solutions proposed by [7] regarding the bearing capacity factors (see Figure 4 [8]). The bearing capacity factors calculated from the CEL model analysis results fall within the range defined by the smooth and rough interface solutions proposed by [7].

ABAQUS/Explicit was used to carry out a two-step analysis. Initially, a geostatic step was conducted to achieve the desired gravity stress level in the soil domain. Subsequently, the pipe section, constrained at the center of the reference point, was laterally displaced at a constant velocity of 0.15 m/s to simulate undrained conditions. This velocity was selected based on the studies by [9,13]. The pipeline’s vertical degree of freedom was released during the lateral displacement process.

It is worth mentioning that, due to the bilateral nature of the contact between the pipe and soil, this model works with either the soil or the pipe as the moving component. For this particular study, the pipeline was chosen as the moving part due to the fact that it is easier to constrain using MPCs. However, soil pressure could also increase if an acceleration was applied to the soil while the pipeline was kept in place.

3.2. Structural Model of the Multi-Layer Flexible Pipe

The structural model was constructed in ANSYS, initially proposed and validated by [38], and improved and applied by [29]. The model’s premise is to simplify the complex geometries of the interlocked carcass and pressure armor using nonlinear material characteristics and contacts between the layers. Even though flexible pipes are composed of many layers, the model is composed of three concentric arcs, which are, in fact, semi-elliptical shapes (Figure 5). Each arc represents one layer of the structural nucleus, the pressure-resistant part of a flexible pipe. The tensile armor layers and the high-strength tape were deliberately removed since they are designed to withstand tension and protect the armors from local buckling, respectively, without contributing significantly to the radial resistance of the pipe. The model was built using the plane strain model to optimize the simulation times.

The model incorporates material plasticity for the metallic and polymeric layers, frictional contact between the pressure layer and the liner, and bonded contact between the liner and the carcass. The liner was extruded directly on top of the carcass due to the flexible pipe’s manufacturing process.

3.2.1. Thickness of Layers

The complex geometry of the flexible pipe’s annulus was adapted into a series of equivalent rings in a structural model. To calculate the equivalent thickness, a methodology suggested by [25] was used for both the carcass and pressure armor layers. First, the bending stiffness per unit length of the carcass was calculated by the following equation:

$$ei={\displaystyle \frac{E}{\left(1-{\nu }^2\right)}}{\displaystyle \frac{2I_{Gmin}}{b}}$$

(4)

where $E$ denotes the Young’s modulus, ν represents the Poisson’s ratio, $I_{Gmin}$ represents the minimum inertia moment of the profile cross-section, and $b$ denotes the profile length. The stiffness per unit of length of the equivalent tube (of thickness $t_{eq}$) is calculated as follows:

$$ei={\displaystyle \frac{E}{\left(1-{\nu }^2\right)}}{\displaystyle \frac{t_{eq}^3}{12}}$$

(5)

Notably, the factor of 2 in the numerator of Equation (4) represents the interlocked superposition. This factor could vary based on the profile’s superposition from 1 (no superposition) to 2 (total superposition):

$$ei={\displaystyle \frac{E}{\left(1-{\nu }^2\right)}}{\displaystyle \frac{\left(1+\psi \right)I_{Gmin}}{b}}$$

(6)

in which $\psi$ is the superposition factor, given by:

$$\psi ={\displaystyle \frac{sup1+sup2}{pitch}}$$

(7)

A generic profile of the superposed parts and the pitch is presented in Figure 6.

Combining Equations (6) and (7) and rearranging for the thickness, the expression for the equivalent thickness is obtained:

$$t_{eq}=\sqrt[3]{{\displaystyle \frac{12\left(1+\psi \right)I_{Gmin}}{b}}}$$

(8)

Equation (8) was employed to calculate the equivalent thickness to replace the interlocked carcass and pressure armor in the structural model.

3.2.2. Boundary Conditions

Different boundary conditions with distinct configurations were considered to facilitate different models, as follows:

• Configuration I: On the symmetry plane, a frictionless support was used, and one node of each layer was horizontally restricted to avoid rigid body motion (see Figure 7);
• Configuration II: Deformable Multi-Point Constraints (MPCs) were employed, connecting the carcass’s internal nodes to a master node in the center. This condition allows for the representation of the whole circumference of the flexible pipe section to match the CEL model’s geometry.

Configuration I was adopted for both wet and dry collapse quasi-static models, since they use symmetry. The transient model, on the other hand, required a more profound analysis. Figure 8 shows a graphical comparison between the ovalization and stress results (maximum von Mises for the entire model, averaged in the element) for the two versions of the transient model; the only difference between them is the boundary condition.

As seen in Figure 8, there is an observable difference in the results shown by both models. The “conventional” curve represents the restrained nodes and edges that correspond to “configuration I”, whereas the MPC curve represents “configuration II”. With the MPC configuration, the ovalization tends to be larger and the stresses lower than the alternative conventional boundary conditions. This is expected, since the MPC is less restrictive and generates fewer concentrated stress points than nodal restrictions. Another point to be considered is that the CEL model uses a kinematic coupling between the inner edge of the pipe and a central reference node. The MPC reproduces that condition more accurately, with the advantage of being deformable and allowing relative movements between the nodes. For these reasons, configuration II was selected for the transient model.

3.2.3. Material Properties

The multi-layer flexible pipe’s cross-section was selected as per the API 17B guidelines [28].

Table 1. Material properties of the flexible pipe layers.

Layer	Material	Equivalent Thickness (mm)	Young Modulus (GPa)	Poisson’s Ratio	Yield Strength (MPa)	Tangent Modulus (MPa)
Interlocked Carcass	Carbon Steel	6.32	200	0.3	450	1450
Internal Polymeric Layer	HDPE	7	0.571	0.45	20.74	-
Pressure Layer	Carbon Steel	6.52	200	0.3	450	1450

shows the calculated thicknesses of various layers and the corresponding material selection. For metallic layers, a typical carbon steel material available in the ANSYS library was used. Beyond the elastic portion of the curve, the material was represented as bilinear with kinematic hardening, where a tangent modulus was also defined, in addition to the Young’s modulus. This simple model was extracted from the ANSYS material library, and it is intended to represent a medium carbon steel. The kinematic hardening means that the yield surface is displaced, without increasing in size, accounting for the Bauschinger effect where the properties of material alter as a result of the repeated cycles of loading and unloading. Throughout the transient analysis, a change in loading conditions is expected. A hardening model must therefore be considered. The values used to define the material models are presented in Table 1. Figure 9 and Figure 10 illustrate the stress–strain curves of the carbon steel material and the HDPE, respectively.

3.3. Discussion of the Ovalization Effect on Soil Resistance

The lateral resistance of soil is expected to increase as the ovalization occurs in the cross-section of the pipe. It is intuitive to imagine that an increase in the lateral area due to ovalization under the same lateral pressure would increase the overall force acting on the pipe section. Since a two-way analysis with the CEL model and a deformable cross-section is computationally expensive, a numerical experiment was conducted to verify the effect of a section deformation on the CEL model.

For that purpose, two 2.5″ flexible pipe sections, one circular and another with a 10% API ovalization, were submitted to the same lateral movement while buried in the same soil. Figure 11 shows the difference in the lateral soil resistance exerted on the pipe, where the deformed section increases the lateral soil resistance against the moving section of the pipe.

Figure 12 shows the bearing capacity, which is the lateral force normalized by the soil stiffness ($S_u$) and the pipe diameter. One may notice that the effect is significantly reduced when the force results are normalized.

While the ovalization certainly influences the p–y curves, in quantitative terms, a relatively large ovalization would be necessary to cause a significant force increase. With an ovalization factor of around 10%, an increase of 11% was observed for the maximum force in the load–displacement curve in comparison with the original circle section, as shown in Figure 11. The influence was significantly decreased by normalizing the resistance with soil properties and the diameter, due to the fact that the ovalization alters the diameter, increasing it in the direction perpendicular to the lateral movement. For both cases, the largest difference occurs near a lateral pipe displacement of 0.13 m; it is close to 8.3% for Figure 11, and for Figure 12 the difference is about 5.3%. Further analyses are required for future studies to accurately simulate the ovalization effect on lateral soil resistance. However, this brief study showed that the assumption on which the current methodology is based is reasonable.

4. Parametric Study

As an immediate application of the methodology, a parametric analysis was performed to examine the effects of different parameters, namely the flexible pipe’s internal diameter, seabed soil strength, and burial depth. As shown in

Table 2. Different cases for the parametric study.

Case Name	Internal Diameter (inches)	Seabed Soil Strength (kPa)	Burial Depth Ratio (H/D)
CS-01	22	50	1.92
CS-02	22	25	1.92
CS-03	22	5	1.92
CS-04	6	50	1.92
CS-05	16	50	1.92
CS-06	22	50	2.92
CS-07	22	50	3.92

, these parameters were analyzed individually and compared with the base case (CS-01) using the cross-section ovalization and stress on the metallic layers as output parameters.

Since the two models are independent, only the external diameter of the flexible pipe must remain fixed to guarantee a geometrical match between the models. Every parameter of the structural model can be easily changed. If the external diameter of the flexible must be changed, then multiple CEL analyses will be required.

The CEL analyses were conducted on a workstation with 24 Intel Xeon Cores, 128 GB of RAM, and a Solid State Drive (SSD), taking around six days to complete. The structural FE analyses, on the other hand, were conducted on a notebook with an i7 Intel processor, 8GB of RAM, and an SSD hard drive and took around 15 min for a 32 s transient analysis.

5. Results of CEL Analysis of the Lateral Pipe Displacement

In the primary case, CS-01, the soil stiffness of 50 kPa could represent overconsolidated clay seabed soil. With the significant lateral movement of the pipe section, the shear bands gradually form in front of the pipe and at the rear side to the seabed soil surface, causing the soil resistance on the pipe section to decrease (see Figure 13).

Figure 14 illustrates the extracted soil pressure profiles on the pipe’s cross-section. Each curve represents a snapshot in time during the lateral displacement, where δ is the displacement normalized by the diameter D. The horizontal axis shows the point on the pipe’s surface on which the pressure is acting, where the range from 0 to 0.5 corresponds to the leading edge of the pipeline, and 0.5 to 1 is the trailing edge.

As the pipe section is laterally displaced, the soil pressure on the pipe increases in comparison to the initial starting stage (see Figure 14a,b). As displacement continued, a cavity gradually developed in the rear side of the pipe, and strain-softening enhanced the formation of shear bands in the soil, extending from the front bottom of the pipe to the soil surface. This, in turn, reduces contact pressure on the pipe front upon passing the maximum resistance (see Figure 14b,c). At a later stage, the soil pressure remains at a relatively constant magnitude (see Figure 14d,e).

The soil’s undrained shear strength (s_u = 50 kPa in CS-01, 25 kPa in CS-02, and 5 kPa in CS-03), the pipe section diameter (D = 22 inches in CS-01, 6 inches in CS-04, 7 inches in CS-05), and the burial depth ratio of the pipe (H/D = 1.92 in CS-01, 2.92 in CS-06, 3.92 in CS-07) can substantially impact the formation of shear bands in the soil during deformation. Figure 15 presents a comparison of the shear bands formed in the soil across different case studies.

Comparing Figure 15a,b, the general distribution of shear bands is similar, but a different number of shear bands could be observed in front of the pipe (more shear bands in CS-02). The shear bands in CS-03 differed a lot from those in CS-01 and CS-02, with extended spiral shear bands initiated from the rear bottom of the pipe and reaching the soil surface at the rear side of the pipe. Comparing Figure 15a,d,e, with different pipe diameters, the distribution pattern of the plastic strain is quite similar. Comparing Figure 15a,f,g, which have an increasing burial depth ratio, significant differences can be observed in both the number and distribution of shear bands.

The pipe section’s diameter, the ratio of burial depth, and the soil’s undrained shear strength can also significantly influence the stress profiles around the pipe section. The comparison of stress profiles around the pipe among different study cases is shown in Figure 16.

As indicated in Figure 16a, with a decreased seabed soil stiffness, the magnitude of soil pressure on the front side of the pipe significantly decreases. In the softest seabed (CS-3), the soil surrounded the pipe and a positive soil pressure was retained over a larger area (i.e., regarding the normalized path ratio) on the pipe. The influence of pipe diameter on soil pressure shows scattered results in Figure 16b. Significant increases could be observed in soil pressure magnitude on the pipe with an increase in the burial depth ratio (see Figure 16c).

6. Structural Analysis Results of the Pipe’s Cross-Section

6.1. Dynamic Effects and Time Frame

The soil pressures obtained in the previous sections were applied to the structural transient model of the pipeline cross-section. Each analysis was broken down into several time steps and further split into substeps for load application. At the end of each substep, a specific load increment was added and balanced by the FE program’s solving routine. Each simulation ran from 0 to 60 s, with frames captured at intervals of one-tenth of a second from 0 to 1 s, one frame per second from 1 to 40 s, and every 2 s in the last 20 s.

Since each change in the load condition required multiple substeps to converge, the results were available for hundreds of instants of time. However, the data were filtered to present the results in this article, providing smoother curves. Figure 17 demonstrates this process, with the raw data shown as a thin red curve and the filtered results depicted as a dashed black curve. Additionally, instead of time, the horizontal axis displays the lateral displacement imposed on the flexible pipe, normalized by the pipe diameter.

Due to the use of normalized lateral displacement, the first stage of the simulation (from zero to one second) will be compacted in the subsequent results. This stage includes the geostatic step and the beginning of the pipe’s lateral displacement, which is a relatively unstable process. These initial stages are illustrated in Figure 17. Since the goal of this section is to compare the results for each case, and every case starts similarly, no significant information will be lost. The upcoming sections will present the outcomes from the parameter variations listed in

Table 3. Collapse analysis results.

Soil Strength (kPa)	API Ovalization	DNV Ovalization	Wet Collapse Pressure (MPa)	Dry Collapse Pressure (MPa)
5	0.0539	0.1079	1.2175	3.6800
25	0.4749	0.9497	1.1669	2.2384
50	0.9538	1.9076	1.1191	1.4455
Diameter (inches)	API Ovalization	DNV Ovalization	Wet Collapse Pressure (MPa)	Dry Collapse Pressure (MPa)
6	0.0226	0.0452	23.0325	39.4190
16	0.3709	0.7419	2.8420	6.3686
22	0.9538	1.9076	1.1191	1.4455
Burial Depth Ratio (h/D)	API Ovalization	DNV Ovalization	Wet Collapse Pressure (MPa)	Dry Collapse Pressure (MPa)
1.92	0.9538	1.9076	1.1191	1.4455
2.92	0.9536	1.9069	1.1194	1.3891
3.92	0.8537	1.7070	1.1301	1.4716

.

6.2. Dynamic Effects and Time Frame

The first parametric analysis involved cases CS-01, CS-02, and CS-03. The soil stiffness of the seafloor increased from 5 kPa, representing pre-excavated remolded soil as a backfilling material used in a sufficiently wide trench, to 50 kPa, representing over-consolidated clay. To simplify the problem in this research, it was assumed that the flexible pipe is buried in uniform soil. As a suggestion for future studies, it would be interesting to incorporate the trenching and backfilling effects along with the pipeline–backfill–trench wall interaction mechanisms.

Three different response parameters, namely ovalization, stress in the metallic layers, and reaction force on the boundary condition area were used as the basis for comparing the results. Figure 18 compares the results during the entire simulation time frame. It is important to note that the results were normalized for clarity and potential generalization. Instead of displaying the time, the horizontal axis represents the lateral displacement of the cross-section, divided by the flexible pipe’s external diameter, which is the normalized pipeline displacement.

The type of result variable determines the scale of the vertical axis. The ovalization (API and DNV) discussed at the beginning of this article is a normalized parameter shown as a percentage. At any specific moment, the maximum amount of von Mises equivalent stress is divided by the metallic material’s yield stress, which is 450 MPa. Finally, the reaction force is normalized using the highest value observed across all cases, and is shown separately in Figure 19.

The initial noticeable aspect of the analysis of the ovalization graphs is that all cases complied with the regulation limits. The API (Figure 18a) and DNV (Figure 18b) ovalization limits were 3% for the loaded normal service condition. The maximum value for API was 0.95 and that of DNV was 1.91 for the CS-01 (50 kPa) case. Another key observation is that as soil stiffness increases, overall ovalization also increases. The spacing between the curves becomes more pronounced, and the oscillations are more significant. The amplitude of these oscillations is higher with elevated soil stiffness, which is almost unnoticeable for the 5 kPa case.

The equivalent stress curves in the metallic layers exhibit a similar trend, with values remaining quite low, reaching up to 0.16, i.e., 16% of the yielding limit. This indicates the elastic and fully recoverable nature of all strains. As the soil stiffness increases, the stresses also increase, which is an expected outcome since a higher soil stiffness leads to greater resistance to displacement, resulting in increased loads on the pressure-resistant layers. Figure 18c displays the stresses in the interlocked carcass, while Figure 18d illustrates stresses in the pressure armor, which are comparatively lower than those in the carcass.

The graphs in Figure 19 confirm the trends observed in the previous result parameters. The reaction forces are high with a higher soil stiffness. A particular difference between this and the previous curves is that the oscillation is small, with the 5 kPa curve seeming almost constant.

All cases start with zero ovalization or stress, which is visible at the beginning of the curves. This initial slope is more perceivable in the graph of Figure 17, while in Figure 18a,b, the slope is compacted on the left side of the plot area (as a direct result of the normalization process). As the simulation progresses, the values tend to stabilize and oscillate around a straight line. This stabilization can be attributed to the soil accommodations during the lateral displacement. One of these oscillations calls for attention, the “camel hump” in the middle of the displacement, with two peaks, one at δ/D = 1.462 and the other at δ/D = 2.297.

To determine the cause of that phenomenon, the CEL model results for the soil were revisited, and those peaks were likely to be caused by a combination of shear bands and soil berm formations. Looking at Figure 20, Point 1 marks δ/D = 1.462, and it is possible to notice an increase in the height of the soil berm compared with δ/D = 1.358. This is even more noticeable for the second peak at δ/D = 1.985. The soil berm looks higher than that at δ/D = 1.771 (see zone “A” in Figure 20). On the other hand, there was the rapid formation and propagation of new shear bands (see zone “B” in Figure 19). New shear bands can be seen at δ/D = 1.567 and δ/D = 2.297 for the first and second peaks, respectively. The propagation of new shear bands can be associated with the sudden drop in soil stress in the region immediately in front of the pipe.

The previously mentioned combined effects of the increase in the height of the soil berm and formation of new shear bands is likely to be influenced by the allowance for the separation in the pipe–soil interface (since the cavity formed behind the pipe can influence the deformation of soil). The cavity formation is a feature of the CEL model, and in the real scenario, the cavity space would be occupied by collapsing soil material. As observed in experimental tests [9], the cracks tend to show on the soil surface at the front top of the pipe, and the soil blocks would naturally fall to fill the cavity, which seems exceptionally challenging to simulate with finite element analysis. Furthermore, the cavity formed in the model could be due to the native soil’s high strength and the inaccuracy introduced by the fraction volume computation regime in the CEL model (fraction of material filled in by the Eulerian element).

6.3. Effect of the Flexible Pipe Diameter

The second parametric analysis verified the influence of the flexible pipe diameter. The flexible pipes were manufactured onshore and stored in large reels. Their main characteristic is that they can withstand high axial and pressure loads with a relatively small bending stiffness. Those features directly result from a complex, multilayered structure built under strict tolerances. Consequently, flexible pipes do not have extremely large diameters like rigid pipelines. The commercial brochure for a major manufacturer [39] states that their current range for flexible pipes is between 2 and 22″. Therefore, the diameters selected for this study were 6, 16, and 22″. Figure 21 shows a comparison of the results.

From the curves shown in Figure 21, the first noticeable aspect is that as the internal diameter increases, the ovalization and equivalent stress increase as well. A larger ovalization will lead to larger strains, which will cause higher stresses. The ovalization (both API and DNV) is smaller than the regulation limit of 3% (Figure 21a,b) and a significant difference is observed between the results for medium and large diameters in comparison with small-diameter cases. For the 6″ flexible pipe, the ovalization is almost nonexistent, with a maximum value of 0.0226% for API and 0.1079 for DNV. Similarly to larger diameters, the 6″ flexible pipe ovalization also oscillates with time. Still, it becomes imperceptible when put on the same scale as the larger-diameter results. For the 22″ flexible pipe, the “camel hump” observed in the previous analysis is recurrent. For both ovalization and stresses, it shows up with its distinctive shape, as shown in Figure 20, and the adjacent paragraphs provide a more detailed discussion.

Regarding the stresses in the metallic layers (Figure 21c,d), the soil displacement does not cause any significant stresses on the interlocked carcass or the pressure armor. The maximum values were approximately 16 and 14% of the yielding stress for the carcass and the pressure armor, respectively. Higher stresses on the carcass are typical for rough bore reinforced flexible pipes, with the carcass and the pressure armor sharing the radial loads and generating stresses. The slenderness of the carcass, associated with its complex geometry, is the cause of the higher stresses in comparison with the pressure armor [27]. The stresses in the 6″ flexible pipe were almost negligible, with maximums between 1 and 2% of the yielding stress.

The final comparison between the reaction forces in the different cases is shown in Figure 22. With an increase in the internal diameter, a larger pipe surface is exposed to the soil pressure, resulting in a higher lateral reaction force. There is also a vertical reaction force whose value is relatively small compared to the lateral one. The results are given in the form of a total equivalent force, considering all directions.

From these results, it is reasonable to assume that the diameter plays a fundamental role in flexible pipe damage due to soil displacements. The smaller diameter of the pipe suffered from a small ovalization, as well as almost negligible stresses. Although small, these stresses might be relevant in the long term when considering repetitive loads. During its operational cycle, the flexible pipe is submitted to various cyclic loads, namely internal pressure variations, slugging, and vessel movements. These loads might induce lateral displacement on the buried sections, leading to these small stresses, which may eventually contribute to fatigue damage.

6.4. Effect of the Burial Depth

The third parametric analysis involved a verification of the burial depth’s impact on the flexible pipe’s ovalization, reaction forces, and stresses. Burial depth is crucial in pipeline trenching as it helps mitigate issues such as upheaval buckling and lateral movement, and protects the pipeline against impacts from sea anchors, icebergs, and similar hazards. Since a deeper burial increases the installation costs, the objective is to optimize the burial depth, ensuring the flexible pipe is buried as shallowly as possible while still benefiting from the trench’s protective features.

As shown in Figure 23, the trench depth is presented as a factor of pipe diameter, with three values selected as a basis for comparison: 1.92D, 2.92D, and 3.92D. There are some significant movements between 1.5 and 3 δ/D, with the “camel hump” oscillation (from the CS-01 case) and a single lump or “dromedary hump” on the CS-06 case. As explained in previous sections, this is due to the soil berm formation and the new shear bands propagation (Figure 20).

Unlike the previous parametric analyses, there is no significant spacing between the curves, and they all seem to oscillate around the same straight line. In terms of the maximums, the 1.92D and 2.92D cases practically share the same value of 0.95% API and 1.91% DNV ovalization, whereas the 3.92D case shows a maximum of 0.85% API and 1.7% DNV ovalization (Figure 23a,b). In terms of stresses, the same pattern can be observed. The 1.92D and 2.92D cases have a maximum of 0.16% for the interlocked carcass yield stress and 0.14% for the pressure armor yield stress, whereas the 3.92D case showed 0.15% and 0.13% for the carcass and pressure armor, respectively (Figure 23c,d). So far, these are the most significant values for stress and ovalization observed in the entire study, but they are still under the allowable limits.

The curves in Figure 23 show similar responses during certain intervals of time, allowing for a pattern to be identified and broken into four stages of lateral flexible pipe responses to soil interaction:

• Stage 1 (between 0 and 0.1 δ/D)—Mount from zero, when the initial resistance of the soil is broken and the flexible pipe begins the lateral displacement.
• Stage 2 (between 0.1 and 1.3 δ/D)—The fist soil berms start forming and the sear bands remain the same in number but increase in value, causing minor oscillations.
• Stage 3 (between 1.3 and 2.7 δ/D)—The soil berms reach a maximum height and new shear bands start to propagate, which is seen as the cause of the “camel hump” and “dromedary hump” oscillations.
• Stage 4 (beyond 2.7 δ/D)—Response stabilization, with smoother oscillations.

Provided that there is a long simulation time, it is possible to speculate that stages 2 to 4 will repeat until the movement ceases, simply because the formed soil berms will be left downwind of the moving pipeline and new ones will be formed upwind.

Figure 24 illustrates the curves of the total force reaction for the three cases, with an increasing force reaction for deeper burial depths. The response is also considerably smoother, not showing the oscillations reported in the ovalization and stress graphs, with a slight tendency to decrease.

Regarding the cross-section flexible pipe deformation and stresses (not expenditure-wise), the burial depth does not play a significant role. The only variable that shows a visible change with burial depth is the reaction force; however, this does not appear to result in a notable difference in stress levels within the interlocked carcass and pressure armor.

6.5. Extreme Values Comparison

In terms of the design, the extreme values, including the maximum ovalization, force reaction, and stresses for each case, are particularly important. As shown in Figure 25, each row of plots corresponds to one of the three parametric analyses performed in this study, as distinguished by the labels on the horizontal axes. It should be noted that the curves are not fitted to the results; they merely connect the data points.

With the change in the soil strength between 5 kPa and 50 kPa, the API ovalization maximums increased from 0.05 to around 0.95%, showing a 17 times’ increase in magnitude. Due to DNV’s more conservative approach, the ovalization increased significantly from 0.1 to 1.9% (19 times larger). The stresses also exhibit a tendency to rise with increasing soil stiffness. In the interlocked carcass, the stresses are more significant and increase from 1.14 to nearly 16% of the yield stress, representing a 14 times’ increase. At the same time, the pressure armor showed a similar response, with an overall increase of 14.7 times. Finally, the maximum value of the total reaction force increased sixfold with the rise in soil stiffness.

The flexible pipe internal diameter is also shown to significantly impact the result parameters. From the second row of graphs in Figure 25, the API ovalization goes from 0.02 to 0.95% and DNV goes from 0.04 to 1.9% with a diameter increase from 6 to 22 inches. This represents an increase of 42 times for both API and DNV ovalization. The stresses showed an increase of 11.5 times for the interlocked carcass and an increase of nearly 15 times for the pressure armor. Despite the relatively large difference, none of the values surpassed the yield limit of 450 MPa or the maximum allowed ovalization of 3%. The force reaction value increases by about 3.5 times, with a diameter increase of 3.7 times.

In the previous section, alterations in the burial depth did not appear to significantly affect the selected parameters of the results for this analysis. At first glance, a comparison of the extreme values might suggest otherwise, as the curves in the third row of Figure 25 indicate some variation with increased burial depth. However, this impression is misleading and is the result of the graphs’ scale. The ovalization parameter shows a 12% overall decrease (API and DNV), which is quite small compared to the effect of the internal diameter, for instance. The stresses show a decrease of 7.2 and 8.7% for the interlocked carcass and pressure armor, respectively, which is also small compared to the 11.5 and 15 times increase observed in the internal diameter study. The force reaction is the parameter that altered the most with increased burial depth. With a depth increase of almost two times, the maximum total reaction force increased by 32%. This is the only parameter that shows some alteration compared to the ovalization and the stresses. However, the scale of the changes is too small for any conclusion regarding the effect of the burial depth on the first two parameters.

6.6. Collapse Analysis

The ovalization results obtained from the previous steps may be used, with the structural model converted to configuration 01, for a collapse analysis (Figure 7). Depending on the type of collapse, the pressure load is exerted either on the surface of the pressure armor or the internal liner. Since the Newton–Raphson method is limited by the requirement of a monotonically increasing force application, the solver is manually switched to the Arc-length method. Despite its being more unstable, this solver allows for negative load increments, which usually appear during the instability analyses. An example of the final results is shown in Figure 26.

Besides the visualization of the flexible pipe’s collapsed shape, the primary objective of the analysis is to investigate the critical load, since this marks the moment of structural failure. In the finite element analysis, this event is related to a sudden decrease in the stiffness of the structure. The critical load is the moment at which a stable equilibrium state becomes unstable. A perturbation or an initial imperfection of any kind triggers a large, sudden displacement of the structure, usually catastrophic.

Table 4. Collapse analysis results.

Soil Strength (kPa)	API Ovalization	DNV Ovalization	Wet Collapse Pressure (MPa)	Dry Collapse Pressure (MPa)
5	0.0539	0.1079	1.2175	3.6800
25	0.4749	0.9497	1.1669	2.2384
50	0.9538	1.9076	1.1191	1.4455
Diameter (inches)	API Ovalization	DNV Ovalization	Wet Collapse Pressure (MPa)	Dry Collapse Pressure (MPa)
6	0.0226	0.0452	23.0325	39.4190
16	0.3709	0.7419	2.8420	6.3686
22	0.9538	1.9076	1.1191	1.4455
Burial Depth Ratio (h/D)	API Ovalization	DNV Ovalization	Wet Collapse Pressure (MPa)	Dry Collapse Pressure (MPa)
1.92	0.9538	1.9076	1.1191	1.4455
2.92	0.9536	1.9069	1.1194	1.3891
3.92	0.8537	1.7070	1.1301	1.4716

shows the critical load results, which is also called “collapse pressure” in this article.

The columns labeled API ovalization and DNV ovalization show the results of the previous analyses with the CEL and structural transient models. These are the initial imperfections imposed on the collapse models. First, an increase in soil stiffness leads to higher levels of ovalization, reducing the critical load of the structure. This means the flexible pipe will collapse with a smaller pressure load after a lateral displacement on high-strength soil. A similar response can be observed from the diameter results: larger diameters suffer from a more significant deformation during the lateral displacement, reducing the collapse pressure value. The burial depth, however, had minimal influence on the overall deformation of the flexible pipe cross-section. This also reflects the small change observed in the wet and dry collapse pressure values.

The wet and dry collapse of flexible pipes is widely covered in the literature. The results presented in Table 3 are in accordance with the findings of [27,29], whose results showed that increasing the initial imperfection results in a reduction in the collapse pressure. Using the present work, it is possible to directly associate the lateral movement of buried flexible pipes with subsequent collapse failure.

7. Conclusions

A methodology was developed to determine the level of damage inflicted on a buried multilayer flexible pipe due to the large lateral displacements. The methodology employs a CEL model for lateral soil resistance and a transient structural model for the flexible pipe. The CEL pipe–soil interaction model was created, utilizing the modified Tresca model [30,31] to account for the strain-softening effects of the seabed soil. Subsequently, the pressure profiles around the pipe section, calculated in the CEL model, were integrated into the FE model featuring simplified metallic interlocked layers. Both models were individually verified against analytical solutions [40], fully tridimensional numerical models [27], and centrifuge tests [8,13].

The decoupled numerical models address the challenges associated with simulating the large lateral displacements of buried flexible pipes with complex multiple layers, and the use of the equivalent ring for the flexible pipe’s complex layers significantly reduces the simulation time. The CEL model provides a more precise modeling of the soil material properties compared to conventional Lagrangian models and is more suitable for large displacement analyses. Soil stiffness significantly affects the flexible pipe’s ovalization, stresses, and reaction forces, with stiffer soils leading to larger fluctuations due to the soil berm formations and new shear band propagation. The pipe diameter has a profound impact on response to lateral movements, with maximum ovalization values differing by a factor of 42 and stresses being 15 times higher in larger diameters. Burial depth, however, has a less significant impact, showing minimal differences in ovalization and stresses across various depths, although the reaction forces increase significantly with burial depth. None of the studied cases exceeded the limits for ovalization and stresses according to API and DNV regulations, reducing the risk of collapse due to high external pressure. A collapse analysis corroborates the ovalization and stress findings, showing that soil strength and pipe diameter substantially influence the collapse pressure, while burial depth has little effect. Although validated, the methodology requires further verification to confirm the validity of the decoupled model assumption, which may include testing different levels of ovalization or employing a two-way soil–structure interaction model despite the higher computational costs. The methodology does not account for voids left by pipe movement, although the boundary conditions mitigate the impacts on ovalization and stress outcomes. Future advancements could improve the safeguarding of submarine communication cables, which share structural similarities with offshore multilayered flexible pipes.

Even though both models employed in this study were previously validated, further work is needed to verify the range of validity of the primary assumption of this model, which is the decoupling between the pipe–soil interaction analysis. This can be achieved in different manners, including using cross-sections with different levels of ovalization on the pipe–soil analysis or using a two-way soil–structure interaction model, with the downsides of a significantly higher computational cost and convergence issues. Another limitation of the model is that the void left in the wake of the pipe movement is left unoccupied by the soil particles. However, the boundary conditions impede the pipe from going in the opposite direction, minimizing the effects on the ovalization or stress outcomes. Finally, submarine communication cables share many structural features with offshore multi-layered flexible pipes. Further developments in this methodology and models might be of strategic value, improving the safeguarding of vital communication lines.