Numerical Modeling of Composite Load-Induced Seabed Response around a Suction Anchor

Abstract

Suction anchors play a crucial role as marine supporting infrastructure within mooring systems. In engineering practice, the composite load comprising nonlinear waves and cyclic pull-out loads can have adverse effects on the seabed soil, posing a threat to the pull-out bearing capacity of the suction anchor. While existing research predominantly focuses on cyclic pull-out loads, the influence of nonlinear wave actions at the seabed surface remains overlooked. This study employs a two-dimensional integrated numerical model to investigate the dynamic soil response around a suction anchor under the influence of both nonlinear waves and cyclic pull-out loads, focusing on the mechanisms that lead to liquefaction and the deterioration of the interfacial friction due to the excess pore pressure buildup. The numerical results reveal that the cyclic pull-out load is the primary factor in the deterioration of the frictional resistance at the suction–soil interface, especially when the pull-out load is inclined with the suction anchor. Parametric studies indicate that the relative difference in frictional resistance deterioration between cases considering and excluding surface water waves becomes more pronounced in soils characterized by a small consolidation coefficient ($C_v$) and relative density ($D_r$).

1. Introduction

Escalating offshore activities have heightened concerns about the stability of floating structures, with a key focus on the critical roles of the mooring system and anchors. Anchors are securely installed at the seabed through a variety of methods, ensuring their stability and resistance to demanding forces. Suction anchors are steel structures designed in the shape of giant buckets that are horizontally oriented, with an open bottom and a closed top. Compared to conventional anchors and drag-embedded anchors, suction anchors possess significant advantages, including ease of handling during field installation, precise positioning, greater pull-out capacity, and adaptability to complex loading conditions [1,2]. As suction anchors are more commonly used in nearshore areas, the challenges of wave actions become a significant concern [3,4,5]. This is mainly because wave-induced cyclic loads can weaken the soil around the anchors, reducing their pull-out capacity. Such effects are especially pronounced when these waves coincide with low-frequency cyclic loads. To address these challenges, a comprehensive understanding of the wave-induced soil response and the frictional behavior at the structure–soil interface is imperative. Armed with this understanding, engineers can devise strategies to fortify anchoring systems, thereby enhancing the stability and performance of offshore structures.

In contemporary research, extensive field studies have been conducted on the bearing capacity of suction caissons, their loading points, and soil characteristics, as documented by scholars such as Andersen et al. [6], Dyvik et al. [7], and Ravichandran et al. [8]. Parallel to these discussions, Burg and Bang [9] highlighted factors affecting the horizontal loading capacity of suction caissons, including the point of horizontal load application, flange and pile diameters, seafloor soil strength, pile diameter-to-length ratio, and properties of layered soils. Apart from field tests, laboratory experiments, particularly centrifuge tests and 1 g model tests, have played a crucial role in gaining a deeper understanding of the mechanical behavior of suction anchors [10,11,12,13,14,15,16,17]. Zhu et al. [18] performed an exhaustive investigation with approximately 10,000 cycles of 1 g laboratory cyclic loading tests, revealing a positive correlation between long-term settlement, cycle count, and cyclic loading amplitude, with the most significant settlement observed in two-way cycling conditions. Bang et al. [19] conducted a series of centrifuge model experiments to assess how the load inclination angle and the mooring line attachment point affect the load-bearing capacity of a suction caisson model, determining that the optimal loading positions are typically between 70% and 75% of the suction caisson’s height. Guo et al. [20] investigated the behavior of suction caissons under vertical loads through a series of model tests and indicated that the passive pore pressure inside the caisson can result in the simultaneous movement of the interior soil plug and the caisson itself. To better understand the anti-uplift behavior of suction caissons, Kelly et al. [21,22] conducted tests with vertical monotonic and cyclic loading on a model and discerned a positive correlation between positive pore pressure and cycling rate.

Compared with field tests and laboratory experiments, numerical simulations have become the predominant method for investigating the dynamic responses of the seabed around suction caissons. In recent years, the finite element method (FEM) has been widely utilized for the study of soil behavior in the vicinity of suction anchors [23,24,25,26]. Specifically, Fu et al. [25] utilized ABAQUS finite element analysis software to assess the impact of tension gaps on the holding capacity of suction anchors in offshore settings, focusing on determining optimal padeye positions and offering practical design strategies to address these gaps. Cheng et al. [26] conducted a detailed study on the performance of caisson anchors in sandy soils under vertical–horizontal loads, utilizing finite element analysis to evaluate the effects of factors like padeye location and soil properties. Andersen and Jostad [27] proposed the limiting equilibrium model of the inclined pull-out capacity of suction anchors under undrained conditions. Ukritchon et al. [28] used the 2D axisymmetric finite element software PLAXIS-2D version 8 to study the undrained pull-out capacity of suction caisson in clays with homogeneous strength and strength linearly increasing with depth. Based on the finite element method (FEM), numerous studies have also analyzed and captured the accumulation of pore water pressure around offshore foundations [29,30,31,32,33,34,35,36,37]. The potential threats arising from the dynamic response of pore fluids surrounding suction anchors have also captured the attention of scholars and prompted extensive discussion. For instance, Thieken et al. [38] presented numerical simulations of a suction caisson, focusing on coupled pore fluid diffusion under fluctuating tension conditions in sandy terrain. Shortly thereafter, Cerfontaine et al. [39] applied the Prevost elastic–plastic model to explore the monotonic and cyclic behaviors of suction caissons under vertical transient loading. Shen et al. [40] proposed a numerical model to predict the oscillation and accumulations of pore pressures around the suction anchor under the influence of long-term vertical cyclic load variations. Recently, Moghaddam et al. [41] employed the CycliqCPSP soil model alongside Biot’s consolidation and linear wave theories to investigate wave-induced liquefaction around suction caissons. However, the aforementioned studies predominantly concentrated on the stability of suction caissons under pull=out cyclic loading conditions, with a notable omission of wave conditions that significantly influence the seabed surrounding these structures. This overlooks a critical aspect of the operational environment, as in reality, multifaceted external loads, particularly under severe meteorological scenarios like storm surges, can intensify destructive impacts and potentially result in unforeseen damage. This gap in research highlights the necessity for comprehensive investigations that incorporate the dynamic interplay of wave conditions and seabed interactions in the assessment of suction caisson stability.

This study extensively investigates seabed responses near the suction anchor under the combined influence of nonlinear wave loading and vertical or inclined pull-out cyclic loading. The presented model employs an one-way coupling, integrating an OpenFOAM v8 flow model with a COMSOL multiphysics 5.4 seabed model. A comprehensive parametric study was carried out to assess how varying wave and soil conditions impact liquefaction susceptibility and to evaluate frictional interactions at the suction–soil interface.

2. Theoretical Formulations

This section outlines the theoretical framework and the governing equations central to the multiphysics model, encompassing the flow sub-model, seabed sub-model, and structure sub-model. Each sub-model is developed to capture the distinct physical phenomena and interactions within its domain. To illustrate the integrated fluid–soil model graphically, the schematic sketch of wave–seabed–structural interactions is provided in Figure 1. This schematic illustration methodically demarcates the wave generation and absorption zones, the air–wave interface, and the intricacies of the structural and seabed sub-models. The governing equations for each domain, such as the Reynolds-averaged Navier–Stokes (RANS) equations for fluid dynamics and Biot’s consolidation equation for seabed mechanics, are also included in the figure.

2.1. Flow Sub-Model

With the premise of incompressible and continuous fluid dynamics in three-dimensional space, the Eulerian approach is utilized to model the scenario in which gravity acts along the negative direction of the y-axis, and the wave profile propagates across the $x-y$ plane. The following governing equations are used to simulate the motion of the two-phase incompressible flow, encapsulating the conservation laws of mass and momentum, which are written as

$$\frac{\partial \langle u_i\rangle }{\partial x_i}=0,$$

(1)

$$\frac{\partial \rho \langle u_i\rangle }{\partial t}+\frac{\partial }{\partial x_j}\left[\frac{1}{n}\rho \langle u_i\rangle \langle u_j\rangle \right]=-n\frac{\partial {\langle p^{\ast }\rangle }^f}{\partial x_i}+n{g}_i{X}_j\frac{\partial \rho }{\partial x_i}+\frac{\partial }{\partial x_j}\left[{\mu }_{eff}\frac{\partial \langle u_i\rangle }{\partial x_j}\right]-\left[CT\right],$$

(2)

where $u_i$ is the velocity vector; t is the time; “$\langle \rangle$” and “${\langle \rangle }^f$” denote Darcy’s volume averaging operator and the intrinsic averaging operator, respectively; n represents the porosity, while $p^{\ast }$ is the pseudo-dynamic pressure; and g and X stand for the acceleration of gravity and the position vector, respectively. ${\rho }_f$ is the density of the fluid, which is calculated by ${\rho }_f=\alpha {\rho }_{water}+(1-\alpha ){\rho }_{air}$, where $\alpha$ is the indicator function that represents the quantity of water per unit volume in each cell. In Equation (2), ${\mu }_{eff}$ defines the efficient dynamic viscosity. Thereafter, the last term in Equation (2) regarding the resistance of the porous media can be written as

$$\left[CT\right]=a\langle u_i\rangle +b\left|\langle u\rangle \right|\langle u_i\rangle +c\frac{\partial \langle u_i\rangle }{\partial t},$$

(3)

in which $CT$ is commonly known as the hydraulic gradient; c is typically assigned the value 0.34 kg/m${}^3$ [42]; $a=\alpha \frac{{(1-n)}^3}{n^2}\frac{\mu }{D_{50}^2}$ and $b=\beta (1+\frac{7.5}{KC})\frac{(1-n)}{n^2}\frac{\rho }{D_{50}}$ are empirical drag force coefficients, as characterized by Engelund [43]’s formulas, where $D_{50}$ denotes the mean nominal diameter of the granular material; and $KC$ stands for the Keulegan–Carpenter number.

The volume of fluid (VOF) approach [44,45] is frequently used to track the free surface in most maritime engineering situations where water and air are the only two fluid phases involved and can be written as follows:

$$\frac{\partial \alpha }{\partial t}+\frac{1}{n}\frac{\partial \langle u_i\rangle \alpha }{\partial x_i}+\frac{1}{n}\frac{\partial \langle u_{ci}\rangle \alpha (1-\alpha )}{\partial x_i}=0,$$

(4)

where ${\mathbf{u}}_c$ stands for the compression velocity, defined as $|{\mathbf{u}}_c|=min\left[c_{\alpha }\left|\mathbf{u}\right|,max\left(\right|\mathbf{u}\left|\right)\right]$, and the parameter $c_{\alpha }$ is assigned a default value of 1, in accordance with established practices in maritime engineering simulations [46], effectively balancing numerical stability with the physical accuracy of the interface tracking process. Notably, a larger value of $c_{\alpha }$ enhances the compression of the interface.

2.2. Seabed Sub-Model

The relationship between the effective stress and the pore pressure in the saturated soil domain is expressed as

$${\sigma }_{ij}={\sigma }_{ij}^{\prime }+{\delta }_{ij}{u}_e,$$

(5)

where ${\sigma }_{ij}$ presents total stress; ${\sigma }_{ij}^{\prime }$ is effective stress; $u_e=u_e^{\left(1\right)}+u_e^{\left(2\right)}$ is the excess pore pressure, where $u_e^{\left(1\right)}$ is the oscillatory component and $u_e^{\left(2\right)}$ is the residual component; and ${\delta }_{ij}$ defines the Kronecker delta. It is noted that tension is taken as negative in this work.

Based on the conservation of mass, Biot’s consolidation equation is employed in the soil model, assuming that the soil is a hydraulically isotropic porous medium. In this context, the soil skeleton is regarded as an elastically isotropic material, and the pore fluid is assumed to be compressible and governed by Darcy’s law. As such, the governing equation can be expressed as follows:

$$k_s{\nabla }^2{u}_e^{\left(1\right)}-{\gamma }_w{n}_s{\beta }_s\frac{\partial u_e^{\left(1\right)}}{\partial t}={\gamma }_w\frac{\partial }{\partial t}\left(\frac{\partial u_s}{\partial x}+\frac{\partial v_s}{\partial y}\right),$$

(6)

$$G{\nabla }^2{u}_s+\frac{G}{\left(1-2{\mu }_s\right)}\frac{\partial }{\partial x}\left(\frac{\partial u_s}{\partial x}+\frac{\partial v_s}{\partial y}\right)=\frac{\partial u_e^{\left(1\right)}}{\partial x},$$

(7)

$$G{\nabla }^2{v}_s+\frac{G}{\left(1-2{\mu }_s\right)}\frac{\partial }{\partial y}\left(\frac{\partial u_s}{\partial x}+\frac{\partial v_s}{\partial y}\right)=\frac{\partial u_e^{\left(1\right)}}{\partial y},$$

(8)

where ${\nabla }^2$ is the Laplace operator; $n_s$ is the porosity of the soil; $k_s$ denotes Darcy’s permeability coefficient; ${\gamma }_w$ represents the unit weight of the pore fluid; and $u_s$ and $v_s$ are the soil displacements in the x- and y-directions, respectively.

In (6), according to Yamamoto et al. [47], the compressibility of the pore fluid ${\beta }_s$ can be defined as

$${\beta }_s=\frac{1}{K_f}=\frac{1}{K_w}+\frac{1-S_r}{P_{w0}},$$

(9)

where $K_f$ is the apparent bulk modulus of the pore fluid; $K_w$ is the true bulk modulus of the elasticity of water; $P_{w0}={\gamma }_{w}d$ is the absolute static pressure of pore water; and $S_r$ is the degree of saturation.

In response to cyclic shearing under the action of ocean waves, residual pore pressure can be generated when the drainage condition is impeded. In this study, a simple but workable method is employed by adding a source term to the Biot’s consolidation equations. This source term is directly associated with the phase-resolved shear stress [48] and is responsible for the generation of residual pore pressure, i.e.,

$$\begin{array}{cc}\hfill \frac{k_s}{{\gamma }_w}\frac{G}{\left(1-2{\mu }_s\right)}& \left(\frac{{\partial }^2{u}_e^{\left(2\right)}}{\partial x^2}+\frac{{\partial }^2{u}_e^{\left(2\right)}}{\partial y^2}\right)+\frac{{\sigma }_0^{\prime }}{T}{\left[\frac{\left|{\tau }_{ins}(x,y,t)\right|}{{\alpha }_r{\sigma }_0^{\prime }}\right]}^{-\frac{1}{{\beta }_r}}=\frac{\partial u_e^{\left(2\right)}}{\partial t},\hfill \end{array}$$

(10)

where the quantity $\left|{\tau }_{\mathrm{ins}}(x,y,t)\right|$ represents the instantaneous phase-resolved shear stress, which is intensively related to the incidence of wave phases; T is the wave period; and ${\sigma }_0^{\prime }$ is the mean normal effective stress of seabed foundation before experiencing cyclic shearing. The relative density ($D_r$) is defined by

$$D_r=\frac{{\left(e_s\right)}_{max}-e_s}{{\left(e_s\right)}_{max}-{\left(e_s\right)}_{min}},$$

(11)

where $e_s$ is the void ratio, and ${\left(e_s\right)}_{max}$ and ${\left(e_s\right)}_{min}$ are the maximum and minimum void ratios, respectively. The coefficients ${\alpha }_r$ and ${\beta }_r$ are a first approximation in (10), which can be defined from the empirical expressions [49]

$${\alpha }_r=0.34D_r+0.084,{\beta }_r=0.37D_r-0.46.$$

(12)

2.3. Structure Sub-Model

Normally, steel is selected as the primary material of choice for suction anchors, and its physical properties are isotropic and linearly elastic. Therefore, this study uses the elasticity theory to simulate the suction anchor. Under the conditions of plane strain, the stress–strain relationships are given by Hooke’s Law,

$$\sigma =\mathrm{C}\epsilon ,$$

(13)

The elasticity matrix C in (13) is defined as

$$\mathrm{C}=\frac{E_p}{\left(1+{\mu }_p\right)\left(1-2{\mu }_p\right)}\left(\begin{array}{ccc}1-{\mu }_p& {\mu }_p& 0\\ {\mu }_p& 1-{\mu }_p& 0\\ 0& 0& \frac{1-2{\mu }_p}{2}\end{array}\right),$$

(14)

where $E_p$ and ${\mu }_p$ are the suction anchor Young’s modulus and Poisson’s ratio, respectively. The transient response of the suction anchor is governed by the equations of motion, namely,

$$G_p{\nabla }^2{u}_p+\frac{G_p}{\left(1-2{\mu }_p\right)}\frac{\partial }{\partial x}\left(\frac{\partial u_p}{\partial x}+\frac{\partial v_p}{\partial y}\right)=0,$$

(15)

$$G_p{\nabla }^2{v}_p+\frac{G_p}{\left(1-2{\mu }_p\right)}\frac{\partial }{\partial y}\left(\frac{\partial u_p}{\partial x}+\frac{\partial v_p}{\partial y}\right)=0,$$

(16)

where $u_p$ and $v_p$ are the pile displacement in the x- and y-directions, respectively, and $G_p$ is the suction anchor shear modulus.

3. Numerical Model Setup

3.1. Integration of Sub-Models

In this study, a one-way coupling process is adopted to integrate the hydrodynamic model (within the framework of OpenFOAM v8) and the structural and geotechnical model (COMSOL multiphysics 5.4). The primary focus is on time-matching schemes, nonmatching mesh strategies, and ensuring the continuity of the interface pressure. The one-way coupling technique has been utilized in previous studies for wave–seabed–structure interactions [50,51,52], yielding notable results in achieving a balance between computation cost and accuracy. Figure 2 provides the schematic representation of the one-way coupling process between the two models. Note that as we transition into shallower waters, the interaction between sea water and the porous seabed, influenced by ocean waves, becomes more pronounced. In such scenarios, the significance of the velocity exchange necessitates the adoption of a two-way coupled model [53].

The present study employed a computational model with the domain configuration depicted in Figure 1, comprising suction anchors and seabed soil in a rectangular area with dimensions of 60 m × 120 m. It is noteworthy that to mitigate any computational error arising from the reflecting waves at the lateral boundary, the lateral extent of the domain was set to three times the wavelength [54]. Additionally, sensitivity analyses were performed to optimize the mesh discretization, time steps, and model solver. The anchor body and surrounding soil were encrypted to enhance the accuracy and precision of the results. The mesh around the anchor body is shown in Figure 3, where the geometry of the anchor is defined based on the stocky suction anchors utilized in the Timor Sea [55]. Accordingly, the depth was set to five times the suction anchor length to minimize the impact of boundary conditions on the soil behavior near the suction anchor.

3.2. Boundary Conditions

In the presented model, the following boundary conditions are adopted:

(1)

In the fluid domain, a proactive wave absorption theory is applied to mitigate the re-reflection of incoming waves at the outlet. This is achieved by imposing a reasonable velocity profile on the absorbent boundaries. For a more comprehensive understanding of the specific boundary conditions utilized in Olaflow, a thorough and detailed description can be found in [56].

(2)

Concerning the instantaneous response of the seabed, it is assumed that both the vertical effective normal stress and shear stress are zero at the seabed surface. Additionally, the vertical flow gradient and soil displacement, adhering to the rigid seabed bottom, are set to zero. At the lateral boundaries of the soil domain, the soil skeleton is allowed to slide, and the normal pore pressure gradient is maintained at zero, indicating no flow.

(3)

As for residual seabed response, residual pore water pressure at the seabed surface is set as zero. Additionally, at both the lateral sides and bottom of the seabed, the gradient of residual pore water pressure is set to zero as well.

(4)

The suction anchor is treated as an impermeable structure. This means that gradients of both the residual and oscillatory components at the interface between the suction anchor and the porous seabed are assumed to be zero. Furthermore, a nonslip boundary condition is applied at the structure–soil interface to investigate the structure–soil interactions, which can be mathematically expressed as

$$u=u_p,v=v_p,\frac{\partial u_e^{\left(1\right)}}{\partial x}=\frac{\partial u_e^{\left(2\right)}}{\partial y}=0,$$

(17)

where $u_p$ and $v_p$ are the displacements of the anchor in the horizontal and vertical directions, respectively.

3.3. Model Validation

As outlined previously, the model comprises two primary components: a flow model and a seabed model. The following section details the validation case adopted for each of these components.

• Validation #1: comparison with the laboratory research conducted by Mattioli et al. [57], which investigated near-bed dynamic interactions between a submarine pipeline and regular wave patterns.
• Validation #2: comparison with the geotechnical model from the numerical research conducted by Cuéllar et al. [58] to assess residual pore pressure around the pile structure in the sandy seabed.

Mattioli et al. [57] undertook a comprehensive experimental study focused on the near-bed hydrodynamics in proximity to a submarine pipeline resting a marine seabed. These experiments were conducted within a wave flume, boasting dimensions of 50 m × $1.3$ m × 1 m. Within the flume, the local water level ($d_w$) was consistently held at $0.3$ m, while the wave height ($H_w$) and wave period ($T_w$) were set at $0.1\mathrm{m}$ and $2\mathrm{s}$, respectively. To discern the flow characteristics in the vicinity of the pipeline, Particle Tracking Velocimetry (PTV) was employed, with data further validated and calibrated using an Acoustic Doppler Velocimeter (ADV). Figure 4 displays the experimental and simulated values of the dimensionless horizontal velocity $u^{\ast }$ above the pipeline at depths represented by $z/D$ for wave phases $\varphi$ of ${90}^{\circ }$ and ${180}^{\circ }$. In the presented findings, deviation between the computed and experimental results can be observed, becoming more apparent with ascending values of the non-dimensional parameter. This can be attributed to the use of the $k-ϵ$ model, which might be less accurate at lower Reynolds numbers where the flow is not fully turbulent. Additionally, integrating equations within the viscous sub-layer with a damping function can introduce numerical stiffness, especially with steep velocity gradients. The $\frac{dk}{dz}=0$ boundary condition, recommended for smooth beds [59], aligns better with near-wall physics compared to the $k=0$ wall boundary condition considered in this study. Despite these issues, it is noteworthy that the numerical results demonstrate a significant correlation with the experimental data reported by Mattioli et al. [57], highlighting the flow model’s proficiency in capturing the near-bed flow velocity amidst the interplay of wave actions and submerged structures.

Since only limited laboratory studies on pore pressure responses associated with seabed soil and suction anchors area available in the literature, the present model was compared with that of Cuéllar et al. [58], who considered pile rocking effects. Note that Cuéllar et al. [58]‘s model focuses on pore pressure accumulation under specific loading conditions (cyclic lateral loading) and explicitly excludes the influence of wave pressure on the seabed. The parameters underpinning our simulations are detailed in

Table 1. Input data for Validation #2.

	Value	Unit
Seabed characteristics
Seabed thickness ($h_s$)	40	(m)
Shear modulus ($G_s$)	1.0 $\times {10}^7$	(N/m${}^2$)
Porosity ($n_s$)	0.46	(-)
Soil density $\left({\rho }_s\right)$	2020	(kg/m${}^3$)
Submerged weight of soil $\left({\gamma }^{\prime }\right)$	8.1	(kn/m${}^3$)
Coefficient of lateral earth pressure $\left(K_0\right)$	0.42	(-)
Permeability ($k_s$)	2.5 $\times {10}^{-4}$	(m/s)
Poisson’s ratio (${\mu }_s$)	0.29	(-)
Degree of saturation ($S_r$)	1	(-)
Poisson’s ratio (${\mu }_s$)	0.35	(-)
Relative density $\left(D_r\right)$	0.92	(-)
Monopile parameters
Young’s modulus $\left(E_p\right)$	3.0 $\times {10}^{10}$	(N/m${}^2$)
Embedded depth $\left(l\right)$	30	(m)
Shear modulus ($G_s$)	1.25 $\times {10}^{10}$	(N/m${}^2$)
Pile length $\left(S\right)$	60	(m)
Pile diameter $\left(D\right)$	8	(m)
Poisson’s ratio (${\mu }_p$)	0.2	(-)
Cyclic displacement parameters
Peak displacement $\left(x_p\right)$	65	(mm)
Lateral load on top of monopile $\left(H\right)$	5	(MN)
Period of motion $\left(T_p\right)$	10	(s)

.

In Figure 5, we present the accumulation of pore pressure under the influence of pile rocking motions, comparing the simulated results of the present model with those of Cuéllar et al. [58]. Here, ‘excess pore pressure’ is an aggregate of oscillatory and residual pore pressures. Figure 5a–d illustrate the temporal shifts in pore pressure at depths of 3 m, 4.3 m, 5.1 m, and 7 m below the seabed. The figure reveals a consistent oscillatory pattern in pore pressure, accompanied by a gradual increase in residual pore pressure over time. While both models generally exhibited similar trends in pore pressure, some disparities were apparent, especially at peak and nadir values. These divergences may be attributable to the distinct approaches each model takes for the pile–seabed interface. Nevertheless, the overall similarities emphasize the fidelity and precision of the proposed geotechnical model, especially when simulating the dynamics of residual pore pressure in seabed regions near pile foundations.

4. Results and Discussion

Throughout their service life, suction anchors are subjected to various external loads. Apart from the typical pre-tension in the mooring line and steady loads from environmental factors such as wind, waves, and currents, they are also influenced by low-frequency cyclic pull-out loads. To more accurately simulate the soil response around a suction anchor in its operational state, this section introduces two specific load combinations as computational examples. In case 1, wave load is integrated with vertical cyclic pull-out load, a scenario commonly observed in Tension Leg Platform (TLP) floating production systems, where cyclic loads are transmitted directly to the top of the suction anchors through their mooring lines. Case 2 combines wave load with oblique cyclic pull-out load, typical for semi-submersible floating production platforms. In practical engineering applications, the cyclic loads generated by these types of floating platforms are directly transmitted to the anchor eye on the external wall of the suction anchor via mooring lines. A schematic representation of the two types of pull-out cyclic load is depicted in Figure 6.

4.1. Combined Wave and Cyclic Pull-out Load

This study designates 3600 s as the duration for the composite cyclic load, which is representative of a typical storm’s duration. Based on the presented model, the seabed response of a suction anchor in a shallow water environment, specifically at a water depth of 20 m, was simulated under these stormy conditions. The input parameters for the suction anchor, seabed soil, and waves are listed in

Table 2. Input data of this numerical example.

	Value	Unit
Wave characteristics
Wave height (H)	5 or various	(m)
Wave period (T)	6 or various	(s)
Water depth (d)	20 or various	(m)
Seabed characteristics
Porosity ($n_s$)	0.425	(-)
Permeability ($k_s$)	1.7 $\times {10}^{-5}$ or various	(m/s)
Poisson’s ratio (${\mu }_s$)	0.33	(-)
Seabed thickness ($h_s$)	60	(m)
Seabed length ($L_s$)	240	(m)
Shear modulus (G)	9 $\times {10}^6$ or various	(N/m${}^2$)
Degree of saturation ($S_r$)	0.985 or various	(-)
Relative density of soil ($D_r$)	0.28	(-)
Frictional angle ($\delta$)	30	(${}^{\circ }$)
Suction anchor characteristics
Suction anchor diameter ($D_p$)	5	(m)
Wall thickness ($t_f$)	0.2	(m)
Anchor length ($L_p$)	12	(m)
Shear modulus of suction anchor ($G_p$)	2.1 $\times {10}^{11}$	(N/m${}^2$)
Poisson’s ratio (${\mu }_p$)	0.35	(-)
Submerged self-weight of suction anchor ($W^{\prime }$)	3.1 $\times {10}^5$	(N)
Cyclic load ratio (CLR)	12 %	(-)
Period of Load $\left(T_1\right)$	20	(s)

. In Case 1, the calculated static pull-out capacity ($V_s$) was 7.1 $\times {10}^6$ N, nearly 23 times the submerged weight. It is noteworthy that during the simulation, the uplift force applied to the suction anchor consistently exceeded its submerged weight yet remained below its uplift bearing capacity in a fully drained state. To simplify computations and modeling, the cyclic pull-out load was assumed to follow a sine waveform with a 20 s period. Simultaneously, the cyclic load ratio (CLR) exerted on the suction anchor, defined as the ratio of load amplitude (Ac) to the static pull-out resistance of the suction anchor ($V_s$) [40], was set at 12%. A schematic of the composite load can be found in Figure 7 for reference.

4.2. Dynamic Soil Response around the Suction Anchor

Figure 8 shows the oscillatory pore pressure around the anchor under the combined action of wave load and vertical pull-out cyclic load (Case 1) over a typical wave period. Notably, significant negative pore pressure was observed within the anchor caisson, with peaks near the tip of the anchor skirt. However, in shallow soil layers, oscillatory pore pressure appeared to be predominantly influenced by cyclic wave loads, fluctuating in accordance with the trough, nodes (occurring when the wave elevation intersected the mean water level), and crest phases.

Figure 9 illustrates the variations in soil shear stress ${\tau }_{xz}$ along the inner and outer walls of the caisson under the influence of composite loading, consisting of wave load and vertical poll-out load, over two loading cycles. The data were collected from twelve monitoring points, strategically placed at depths ranging from 1/6L to L beneath the seabed surface, with intervals of 1/6L, as depicted in Figure 3.

Due to the combined effects of wave loading and cyclic pull-out loading, the shear stress experienced by the anchor walls reflected a complex pattern of superposition, deviating from a simple periodic waveform, as observed in Shen et al. [40]. This complexity emerged from the varying frequencies and phases of the two contributing loads, which could either reinforce or counteract each other at particular instances, leading to soil behavior diverging from that caused by a single load. Notably, the peak shear stress was observed at a depth of approximately 1 L on both sides of the anchor wall. Inside the anchor, particularly at the tip, the shear stress remained comparatively low, except for a concentration of significant static and dynamic shear stresses and strains induced by soil–structure interaction at the skirt tip of the suction anchor.

Figure 10 presents snapshots of the accumulated residual pore pressures for case 1 after the 5th, 90th, and 180th load cycles. For comparison, snapshots at the same loading cycles without considering surface wave loading are also included. In all cases, the residual pore pressure was more likely to occur in the soil on the external side of the anchor wall, where larger shear deformations were experienced (i.e., Figure 9). A notable observation was that, in comparison to conditions under an isolated tensile load, the residual pore pressure could accumulate to a significant extent under combined loading scenarios. With an increasing number of loading cycles, the difference became more apparent, particularly in the shallow soil layer at the incipient loading stage, spreading both vertically and horizontally along the external side. In contrast, there was no conspicuous accumulation of pore pressure on the internal side of the anchor wall. This was attributed to the generation of negative pore pressures inside the suction anchor when subjected to a pull-out load. These negative pore pressures influenced the accumulation of pore water pressure in the anchor body, as detailed in Luan et al. [60].

Figure 11 displays snapshots of accumulated residual pore pressures for case 2 at four different loading cycles. In this scenario, the oblique pull-out cyclic loads were applied to the anchor eye on the outer right wall of the suction anchor model at 315 degrees. This point was located $\frac{2}{3}{L}_p$ below the seabed. In contrast to case 1, case 2 exhibited clear asymmetric behavior in the residual pore pressure distribution. Furthermore, as the combined cyclic load application persisted, there was a clear accumulation of pore pressure observed on the left internal wall of the anchor. This phenomenon was attributed to the inclined pull-out load applied to the right external wall, inducing excessive deformations in the soil on the left internal wall and thereby facilitating the development of pore pressure in this specific area.

4.3. Frictional Resistance and Liquefaction

When exposed to cyclic pull-out loads, a suction anchor is prone to sliding failure, representing the minimum pull-out capacity of the suction anchor. In drained conditions, the analysis of the static pull-out capacity resistance ($V_s$) of the suction anchor takes into account the self-weight of the submerged suction anchor ($W^{\prime }$) and the soil friction on both the outer and inner sides of the caisson wall:

$$V_s=W^{\prime }+F_{total},$$

(18)

Here, $F_{total}$ is the friction force, which can be decomposed into internal and external components:

$$F_{total}=F_{int}+F_{ext},$$

(19)

$$F_{int}={\int }_0^{L_p}{\sigma }_y^{\prime }dy·(Ktan\delta )\left(\pi D_i\right),$$

(20)

$$F_{ext}={\int }_0^{L_p}{\sigma }_y^{\prime }dy·(Ktan\delta )\left(\pi D_e\right),$$

(21)

where $L_p$ is the anchor length; ${\sigma }_y^{\prime }$ is the effective vertical stress; $D_i$ and $D_e$ represent the internal and external diameters of the anchor, respectively; K stands for the lateral pressure coefficient; and $\delta$ is the mobilized friction angle between the anchor wall and the soil.

In reality, the residual pore pressure accumulating around the anchor can reduce the effective stress of the soil. This, in turn, lessens the friction between the anchor wall and the seabed soil. In order to display the dynamic change in friction force more intuitively, the frictions are normalized by the initial values.

$$F_{total}\ast =\frac{F_{total}}{F_{total}^{\prime }},F_{ext}\ast =\frac{F_{ext}}{F_{ext}^{\prime }},F_{int}\ast =\frac{F_{int}}{F_{int}^{\prime }},$$

(22)

where $F_{total}\ast$, $F_{int}\ast$, and $F_{ext}\ast$ are total, internal, and external friction after normalization; $F_{total}$, $F_{int}$, and $F_{ext}$ are the total, internal, and external friction after applying the wave load; and $F_{total}^{\prime }$, $F_{int}^{\prime }$, and $F_{ext}^{\prime }$ are the total, internal, and external friction before applying the wave load, respectively.

In Figure 12, the time variations of friction along both sides of the suction anchor are illustrated. As anticipated, the friction experienced a decline owing to the accumulation of residual pore pressure, leading to a decrease in the initial effective stresses. Notably, the rate of friction reduction was rapid during the initial stages, followed by a more gradual decrease. Beyond the 3500 s mark, a substantial reduction in friction along the outer wall, amounting to 58%, was observed in contrast to the 20% reduction along the inner wall. Consequently, there was an overall reduction in friction around the caisson wall of 39%.

In Figure 13 and Figure 14, the temporal variations of the normalized total frictional force are depicted under various loading conditions, encompassing both vertical pull-out (case 1) and inclined pull-out load (case 2). Evidently, the primary factor contributing to the reduction in frictional forces in both cases was the cyclic pull-out load, as opposed to surface wave loading. For example, under a vertical pull-out load (case 1) and inclined pull-out load with the same loading amplitude (case 2), there was a reduction of approximately 30% and 36%, respectively. In contrast, wave loading induced a comparatively minor reduction of only 5%. Indeed, it is intriguing that in both cases, a simple mathematical superposition resulted in an underestimation of frictional forces, with reductions of approximately 3% and 15% for case 1 and case 2, respectively. This underscores the imperative of enhancing the current model, particularly when dealing with scenarios involving the simultaneous application of wave and pull-out loads, and especially when the load is inclined with the anchor.

Figure 15 and Figure 16 show the distributed pattern of liquefaction potential at four different loading cycles for the loading scenarios of cases 1 and 2, respectively. The liquefaction potential area was observed to continuously spread outward, albeit at a decreasing rate, with an increasing number of loading cycles. In contrast to the loading scenario of case 1, where liquefaction potential was confined to the external wall of the caisson, in the case 2 loading scenario, liquefaction potential could also manifest on the inner wall of the suction anchor’s left side. This phenomenon was interpreted, as mentioned previously, as a consequence of the inclined pull-out load applied to the right external wall, causing excessive deformations in the soil on the left internal wall of the anchor. Interestingly, for both loading scenarios, the liquefaction zone exhibited asymmetric behavior, induced by either nonlinear wave loads or inclined pull-out loads [61]. While asymmetrical behavior was evident under the case 2 loading scenario, it was less apparent under the case 1 loading scenario. This difference was primarily attributed to the effects of nonlinear wave loads, which play a less significant role in influencing soil behavior compared to the cyclic pull-out load, as discussed earlier.

4.4. Parametric Analysis

The consolidation coefficient ($C_v$) delineates how quickly soil consolidates under applied loads, a critical factor in shaping the residual soil response to ocean waves [49]. Figure 17 illustrates the time variations of the normalized frictional force at four different consolidation coefficients. It was found that the frictional force is more likely to decrease at small consolidation coefficients, where the pore water pressure is less likely to drain out, leading to its accumulation and a substantial decrease in the initial effective stress. Simultaneously, with smaller consolidation coefficients, the difference between composite load and pure pull-out load becomes more significant. This suggests that wave loads are more necessary to consider when the soil has poor drainage conditions and low stiffness.

The relative density of soil ($D_r$) is a parameter that indicates the initial state of the soil, significantly affecting its mechanical behavior when subjected to cyclic loads [62,63]. Figure 18 demonstrates how relative density influences the normalized frictional force under pure pull-out load or composite load. As anticipated, there was a more pronounced reduction in the frictional force observed in soil with a small relative density. This is because soil in a loose state is more susceptible to liquefaction when subjected to cyclic shearing. In such a state, the difference between cases with or without considering the effect of surface wave loading also becomes more significant, as indicated in Figure 18b.

5. Conclusions

In this study, a two-dimensional integrated numerical model was developed to explore the dynamic soil response around a suction anchor under the simultaneous influence of combined wave and cyclic pull-out load. The investigation began with an analysis of the development of oscillatory and residual pore pressure around the caisson under two different composite loading scenarios involving vertical or inclined pull-out load. Subsequently, the study delved into the patterns of liquefaction potential and the reduction of frictional resistance. Parametric studies were conducted to examine the influence of soil parameters on the frictional force under pure pull-out load or composite load. The primary conclusions drawn from these analyses were as follows:

(1)

The oscillatory pore pressure around the caisson was primarily affected by the cyclic pull-out load, with significant negative pore pressure being observed near the tip of the anchor skirt. Meanwhile, surface water waves impacted the oscillatory pore pressure in shallow soil layers, exhibiting fluctuations in accordance with the trough, nodes, and crest phases.

(2)

The residual pore pressure on the external side of the anchor wall was significantly more pronounced than on the inner side, displaying apparent asymmetric behavior in the distributed pattern when the pull-out load was inclined with the anchor. The difference between cases with and without considering surface water waves was apparent in the shallow soil layer at the incipient loading stage, spreading both vertically and horizontally along the external side as the number of loading cycles increased.

(3)

Compared to surface water waves, cyclic pull-out load was the primary factor contributing to the reduction in frictional resistance. The loading scenario wherein pull-out load was inclined with the suction anchor could reduce the frictional force more than the scenario with a vertical pull-out load. A simple mathematical superposition resulted in an overestimation of frictional resistance, particularly for the loading scenario with an inclined pull-out load.

(4)

Parametric analysis indicated that a more pronounced reduction in the frictional force was observed in soil characterized by a small consolidation coefficient and relative density. In these conditions, the relative difference between cases with and without considering surface water waves became significant.