Accumulation of Pore Pressure in a Soft Clay Seabed around a Suction Anchor Subjected to Cyclic Loads

Abstract

A suction anchor is an appealing anchoring solution for floating production. However, the possible effects of residual pore pressure can be rarely found any report so far in term of the research and design. In this study, the residual pore pressure distribution characteristics around the suction anchor subjected to vertical cyclic loads are investigated in a soft clay seabed, and a three-dimensional damage-dependent bounding surface model is also proposed. This model adopts the combined isotropic-kinematic hardening rule to achieve isotropic hardening and kinematic hardening of the boundary surface. The proposed model is validated against triaxial tests on anisotropically consolidated saturated clays and normally consolidated saturated clays. The analytical results show that the excess pore water pressure accumulates primarily on the outside of the suction anchor, whereas negative pore water pressure mainly on the inside. The maximum values of both sides appear in the lower part of the seabed. According to the distribution characteristics of the residual pore pressure, a perforated anchor is proposed to reduce the accumulation of excess pore water pressure. A comparative study generally shows that the perforated anchor can effectively reduce the accumulation of excess pore water pressure.

1. Introduction

Various types of offshore platforms are applied for energy production in deep-water environments. As the water depth increases, the offshore platform that rests upon the seabed and relies on gravity base foundations or traditional pile foundations is uneconomical and impractical. The suction anchor is widely used as a cost-effective mooring foundation for floating production systems [1,2], as shown schematically in Figure 1. In natural marine environments, suction anchors are subjected to cyclic loads produced by the floating platform. With the increase in the cyclic numbers, excess pore water pressure is developed and accumulated in a soft clay seabed; therefore, the soft clay structure is degraded, followed by degradation in strength and a reduction in stiffness [3]. Such a phenomenon may reduce the uplift capacity of the suction anchor, and it is, thus, essential to assess the accumulation of the pore water pressure within the seabed around the suction anchor subjected to vertical cyclic loads.

In previous research, experimental studies attempted to research the response of suction caissons subjected to cyclic loads. Cheng and Wang [4] conducted a series of small-scale laboratory model tests to predict the stability of the suction anchor under cyclic loading conditions. Wallace et al. [5] studied the response of suction caisson foundations with cyclic loads in soft clay by means of centrifuge testing. Dyvik et al. [6] investigated the cyclic response of the suction anchor under vertical loads by carrying out field tests on a model foundation. In general, the experimental study of the seabed response around a suction caisson primarily contains single-gravity mode tests, centrifuge model tests, and reduced scale field tests. Only a small number of field tests on suction anchors have been reported in the previously published literature [7]. The pore pressure response in soft clay is difficult to measure accurately under cyclic loading conditions. This is because scaling rules for different aspects of behavior are conflicting in small models [8]. This is the reason for the general lack of available data on pore pressure in the existing model experiments [9,10,11,12]. Currently, a direct comparative analysis seems impossible. A reasonable numerical analysis model needs to be developed to capture the accumulation of pore water pressure.

To date, considerable effort has been devoted to investigating the pore pressure response within the seabed around the offshore foundation by numerical methods [13,14]. The first method is the semi-empirical method, which is based on the simplified elastic-plastic model and predicts the pore pressure response around the offshore foundation subjected to cyclic loads [15,16]. The method is simple, which is suitable for the analysis of the oscillatory pore pressure response induced by waves acting on the sandy seabed. The second method is to apply Biot’s poro-elastic theory for soil models to describe the pore pressure response [17,18,19,20,21]. Based on Biot’s poro-elastic theory, Shen et al. [22] first explored the pore-pressure response in the soil around the suction anchor under cyclic loading conditions, but this was limited to the sandy seabed.

For soft clay seabeds, the classical elasto-plastic theory is based on the assumption that the soil is completely elastic within the yielding surface. It fails to describe the accumulation of pore water pressure, nonlinearity of modulus, and other cyclic behaviors of soft clay [23]. The bounding surface theory based on plastic hardening modulus presents a general framework to describe the cyclic behaviors of soft clay. However, the conventional bounding surface models [24,25] are generally based on the assumption that the mapping origin is fixed at the coordinate origin and the unloading stage is completely elastic, thereby failing to capture the real soil behaviors, such as cyclic stiffness degradation and initial anisotropy, under cyclic loads.

This study proposes a bounding surface plasticity model with the mixed hardening rule to predict the development of pore water pressure within the seabed around the suction anchor under vertical cyclic loading conditions. A damage parameter and initial anisotropic tensors are introduced into the bounding surface model, to represent the remolding of the soil structure and initial anisotropy, respectively. The present model is efficient at capturing the development of pore water pressure in the soft clay seabed subjected to cyclic loads. It is validated against available experimental laboratory data. Subsequently, the influences of the load amplitude and soil material on the distribution of the residual pore pressure around the suction anchor are examined. Finally, according to distribution characteristics of the residual pore pressure, an improved rational structure named the perforated suction anchor is proposed.

2. Theoretical Formulations and Numerical Approach

2.1. Constitutive Model

2.1.1. Anisotropic Bounding Surface

Based on the two-dimensional bounding surface model developed by Hu and Liu [26] for the cyclic dynamic analysis of saturated clay, a three-dimensional form is established by the generalized Mises criterion method, $q=\sqrt{3\left({\overline{s}}_{ij}-p{\alpha }_{ij}\right)\left({\overline{s}}_{ij}-p{\alpha }_{ij}\right)/2}$. The mathematical equation of the bounding surface is described as:

$$F_m=(\overline{p}-{\xi }_p^{(m)}{)}^2-(\overline{p}-{\xi }_p^{(m)})p_c^{\left(m\right)}+\frac{3}{2(M-{\alpha }_0^2)}\left[{\overline{s}}_{ij}-{\xi }_{ij}^{\left(m\right)}-(\overline{p}-{\xi }_p^{\left(m\right)}){\alpha }_{ij}^0\right]·\left[{\overline{s}}_{ij}-{\xi }_{ij}^{\left(m\right)}-(\overline{p}-{\xi }_p^{\left(m\right)}){\alpha }_{ij}^0\right]$$

(1)

in which $\left(\overline{p},{\overline{s}}_{ij}\right)$ denotes the image stress point on the bounding surface, as shown in Figure 2; $p_c$ is the reference size of the bounding surface; $\left({\xi }_p^{\left(m\right)},{\xi }_{ij}^{\left(m\right)}\right)$ is the coordinate of the endpoint; M represents the slope of the critical state line (CSL) in triaxial space; and ${\alpha }_0$ denotes the inclination of the bounding surface, which represents the degree of soil anisotropy:

$${\alpha }_0=\sqrt{\frac{3}{2}{\alpha }_{ij}^0{\alpha }_{ij}^0}.$$

(2)

For consolidated samples, the initial stress ratio is $K_0={\sigma }_3^0/{\sigma }_1^0$. The anisotropic tensors can be defined as follows:

$${\alpha }_{11}^0=\frac{2(1-K_0)}{1+2K_0}, {\alpha }_{22}^0={\alpha }_{33}^0=\frac{(K_0-1)}{1+2K_0}$$

(3)

2.1.2. The Evolution of the Boundary Surface

The generalized isotropic hardening rule is applied in the present model, which assumes that the boundary surface isotropically hardens around the discrete mapping center in the stress space. When the loading path changes its direction, the boundary surface should translate along the direction of stress reversal point to the image stress point. The mapping center can be expressed as:

$$\left(o_{p,n+1},o_{ij,n+1}\right)=\{\begin{array}{cc}\left(o_{p,n},o_{ij,n}\right)& \frac{\partial F}{\partial {\overline{\sigma }}_{ij,n}}d{\sigma }_{ij,n+1}\ge 0\\ \left(p_n,s_{ij,n}\right)& \frac{\partial F}{\partial {\overline{\sigma }}_{ij,n}}d{\sigma }_{ij,n+1}<0\end{array},$$

(4)

in which $\left(o_{p,n},o_{ij,n}\right)$ and $\left(o_{p,n+1},o_{ij,n+1}\right)$ are the coordinates of the mapping center at incremental steps of n and n+1, respectively. The location of the boundary surface depends on the stress path direction.

1) When the stress path changes its direction:

$$\begin{array}{c}{\xi }_p^{(m+1)}={\xi }_p^{\left(m\right)}(p-\overline{p})\\ {\xi }_{ij}^{(m+1)}={\xi }_{ij}^{\left(m\right)}(s_{ij}-{\overline{s}}_{ij})\end{array}\}.$$

(5)

2) When the stress path does not change its direction:

$$\begin{array}{c}{\xi }_{p,n+1}^{(m+1)}=o_p^{(m+1)}+({\xi }_{p,n}^{(m+1)}-o_p^{(m+1)})\frac{p_{c,n+1}^{(m+1)}}{p_{c,n}^{(m+1)}}\\ {\xi }_{ij,n+1}^{(m+1)}=o_{ij}^{(m+1)}+({\xi }_{ij,n}^{(m+1)}-o_{ij}^{(m+1)})\frac{p_{c,n+1}^{(m+1)}}{p_{c,n}^{(m+1)}}\end{array}\},$$

(6)

in which $\left({\xi }_p^{(m+1)},{\xi }_{ij}^{(m+1)}\right)$, $\left({\xi }_{p,n}^{(m+1)},{\xi }_{ij,n}^{(m+1)}\right)$, and $\left({\xi }_{p,n+1}^{(m+1)},{\xi }_{ij,n+1}^{(m+1)}\right)$ represent the coordinates of endpoints at incremental steps of 0, n, and (n + 1) in the (m + 1)th loading events, respectively. To describe the degradation of soft clay under cyclic loading conditions, a damage parameter is introduced into the isotropic hardening rule. The evolution of the size of the boundary surface depends on plastic volumetric strains and the damage parameter, which can be defined as:

$$p_{c,n+1}=p_{c,n}exp(\frac{1+e_0}{\lambda -\kappa }d{\epsilon }_{v,n+1}^p){\omega }_{n+1},$$

(7)

$${\omega }_{n+1}=exp(-\beta {\epsilon }_A),$$

(8)

where

$${\epsilon }_A={\displaystyle \int \sqrt{2d{e}_{ij,n+1}^{p}d{e}_{ij,n+1}^p/3}},$$

(9)

in which λ and κ denote the compression index and swelling index in the $e-lnp$ space, respectively; $e_0$ represents the void ratio after consolidation under $p=p_c$. The state variable $\omega$ is a function of the deviatoric plastic strain $e_{ij}^p$, which is decreasing with the accumulated $e_{ij}^p$. The model parameter $\beta$ controls the rate of damage accumulation. The decrease in $\omega$ represents the degradation in stiffness of the clay structure.

2.1.3. Mapping and Flow Rules

The classical radial mapping criterion developed by Dafalias [27] is used in the model, due to its effectiveness. The conventional radial mapping criterion was modified for soft clay in this research. The mapping rule is shown in Figure 3 and can be expressed as:

$$\begin{array}{c}\overline{p}=b(p-o_p)+o_p\\ {\overline{s}}_{ij}=b(s_{ij}-o_{ij})+o_{ij}\end{array}\},$$

(10)

where b is a scalar factor, that can be expressed as:

$$b=\frac{{\delta }_0}{{\delta }_0-\delta },$$

(11)

in which ${\delta }_0-\delta$ denotes the distance between the current stress point and the mapping center, and ${\delta }_0$ indicates the distance from the current stress point to the image stress point. The loading index is calculated by imposing the consistency condition to its corresponding bounding surface equation:

$$\Lambda =L=\frac{1}{{\overline{K}}_p}\left(\frac{\partial F}{\partial \overline{p}}d\overline{p}+\frac{\partial F}{\partial {\overline{s}}_{ij}}d{\overline{s}}_{ij}\right),$$

(12)

in which ${\overline{K}}_p$ is the plastic modulus at image stress states.

2.1.4. Incremental Equations

The elastic volumetric strain increment and elastic shear strain increment can be given by:

$$d{\epsilon }_v^e=\frac{dp}{K}, d{e}_{ij}^e=\frac{d{s}_{ij}}{2G}$$

(13)

It is assumed that the bulk modulus and shear modulus can be obtained by:

$$K=\frac{1+e_0}{\kappa }p, G=\frac{3K(1-2\nu )}{2(1+\nu )},$$

(14)

in which, ν denotes Poisson’s ratio. Here, an associated flow rule is used, i.e., the plastic strain increment vector is always normal to the yield surface, with Л and L coinciding. The plastic strain direction is defined by the boundary. The plastic constitutive relations can be defined as:

$$d{\epsilon }_v^p=\Lambda \frac{\partial F}{\partial \overline{p}},$$

(15)

$$d{e}_{ij}^p=L\frac{\partial F}{\partial {\overline{s}}_{ij}},$$

(16)

2.1.5. Hardening Modulus

$K_p$ and ${\overline{K}}_p$ are related to the actual plastic modulus and the bounding plastic modulus at the image stress ${\overline{s}}_{ij}$, respectively. According to the consistency condition, the following equation can be obtained:

$$\frac{\partial F}{\partial \overline{p}}d\overline{p}+\frac{\partial F}{\partial {\overline{s}}_{ij}}d{\overline{s}}_{ij}+\frac{\partial F}{\partial {\epsilon }_v^p}d{\epsilon }_v^p+\frac{\partial F}{\partial \omega }d\omega =0.$$

(17)

Substituting the hardening rules, the loading index and the equation of the bounding surface into the consistency conditions, the bounding plastic modulus at the image stress can be written as:

$${\overline{K}}_p=(\overline{p}-{\xi }_p)p_c{\overline{n}}_p({\chi }_0-\frac{2\beta }{M^2-{\alpha }_0^2}\frac{\overline{q}}{{\overline{n}}_p}),$$

(18a)

in which:

$${\chi }_0=\frac{1+e_0}{\lambda -\kappa },$$

(18b)

$$\overline{q}=\sqrt{\frac{3}{2}{\widehat{s}}_{ij}{\widehat{s}}_{ij}},$$

(18c)

$${\overline{n}}_p=\frac{\partial F}{\partial \overline{p}},$$

(18d)

$${\widehat{s}}_{ij}={\overline{s}}_{ij}-{\xi }_{ij}-(\overline{p}-{\xi }_p).$$

(18e)

Here, the interpolation is adopted to calculate the plastic-hardening modulus at the current stress state.

$$K_p={\overline{K}}_p+H(\overline{p},\overline{q},{\epsilon }_v^p,\omega )(\frac{\delta }{{\delta }_0-\delta }{)}^r,$$

(19)

in which $H(\overline{p},\overline{q},{\epsilon }_v^p,\omega )$ is the shape hardening function. Different shape hardening functions are adopted in the first loading, reloading, and unloading stages, respectively.

$$H(\overline{p},\overline{q},{\epsilon }_v^p,\omega )=\{\begin{array}{c}\left|K_m-{\overline{K}}_p\right|\begin{array}{cccc}& \mathrm{for}& \mathrm{first}& \mathrm{loading}\end{array}\\ \left|{\zeta }_r{K}_m-{\overline{K}}_p\right|\begin{array}{ccc}& \mathrm{for}& \mathrm{reloading}\end{array}\\ \left|{\zeta }_u{K}_m-{\overline{K}}_p\right|\begin{array}{ccc}& \mathrm{for}& \mathrm{unloading}\end{array}\end{array}$$

(20a)

in which:

$${\zeta }_u=(1+\frac{\partial F}{\partial \overline{p}}/\eta ){\zeta }_r,$$

(20b)

$$K_m=8{\chi }_0{\left(p_c\right)}^3,$$

(20c)

where $\gamma$ and $\eta$ are model parameters. The material parameter ${\zeta }_r$ controls the reloading events, of which a detailed account of physical meaning can be found in the paper by Hu et al. [26].

2.1.6. Implicit Integration Algorithm

In this section, the implicit integration algorithm is used for the implementation of the models. The calculation steps are as follows:

1)

Assume that strain increment $\Delta {\epsilon }_{n+1}$ is complete elastic increment at incremental steps of n+1. For initial iteration count k = 0, these variables are defined as follows:

$$\begin{array}{ccc}{\Lambda }_{n+1}^{\left(0\right)}=0,& \Delta {\omega }_{n+1}^{\left(0\right)}=0,& \Delta {\epsilon }_{v,n+1}^{\begin{array}{cc}p& \left(0\right)\end{array}}=0\\ \Delta e_{ij,n+1}^{\begin{array}{cc}p& \left(0\right)\end{array}}=0,& b_{n+1}^{\left(0\right)}=b_n,& p_{c,n+1}^{\left(0\right)}=p_{c,n}\\ {\overline{K}}_{p,n+1}^{\left(0\right)}={\overline{K}}_{p,n},& {\xi }_{ij,n+1}^{\left(0\right)}={\xi }_{ij,n},& o_{ij,n+1}^{\left(0\right)}=o_{ij,n}\end{array}\},$$

(21)

2)

Non-linear elastic predictor:

$$\begin{array}{l}{p}_{n+1}^{(0)}=p_{n}exp\left(\frac{1+e_0}{\kappa }\Delta {\epsilon }_{v,n+1}\right)\\ s_{ij,n+1}^{(0)}=s_{ij,n}+2G_{n+1}^{(0)}\Delta e_{ij,n+1}\end{array}\},$$

(22)

$$\begin{array}{cc}{K}_{n+1}^{(0)}=p_n\frac{1+e_0}{\kappa },& G_{n+1}^{(0)}\end{array}=\frac{3K_{n+1}^{(0)}(1-2\nu )}{2(1+\nu )}.$$

(23)

3)

Distinguish the unloading process from the loading event, according to Equation (4). (a) reloading, homological center remains constant. (b) unloading, update the homological center and the bounding surface, according to Equations (5)–(6). Then evaluate the following residuals:

$$r_{l,n+1}^{\left(k\right)}=\{\begin{array}{l}{p}_{n+1}^{\left(k\right)}-p_{n}exp\left[\frac{1+e_0}{\kappa }\left(\Delta {\epsilon }_{v,n+1}-\Delta {\epsilon }_{v,n+1}^{\begin{array}{cc}p& \left(k\right)\end{array}}\right)\right]\\ \Delta {\epsilon }_{v,n+1}^{\begin{array}{cc}p& \left(k\right)\end{array}}-{\Lambda }_{n+1}^{\left(k\right)}{\overline{n}}_{p,n+1}^{\left(k\right)}\\ s_{ij,n+1}^{\left(k\right)}-s_{ij,n}-2G_{n+1}^{\left(k\right)}(\Delta e_{ij,n+1}-\Delta e_{ij,n+1}^{\begin{array}{cc}p& (k)\end{array}})\\ \Delta e_{ij,n+1}^{\begin{array}{cc}p& (k)\end{array}}-{\Lambda }_{n+1}^{\left(k\right)}{\overline{n}}_{s,n+1}^{\left(k\right)}\\ p_{c,n+1}^{\left(k\right)}-p_{c,n}exp\left(\frac{1+e_0}{\lambda -\kappa }\Delta {\epsilon }_{v,n+1}^{\begin{array}{cc}p& \left(k\right)\end{array}}\right){\omega }_{n+1}^{\left(k\right)}\\ {\omega }_{n+1}^{\left(k\right)}-{\omega }_{n}exp\left(-\beta \sqrt{\frac{2}{3}\Delta e_{ij,n+1}^{\begin{array}{cc}p& (k)\end{array}}\Delta e_{ij,n+1}^{\begin{array}{cc}p& (k)\end{array}}}\right)\\ b_{n+1}^{\left(k\right)}-b_n-{\Lambda }_{n+1}^{\left(k\right)}\frac{{\overline{K}}_{p,n+1}^{\left(k\right)}-b_{n+1}^{\left(k\right)}{K}_{p,n+1}^{\left(k\right)}}{A_{n+1}^{\left(k\right)}}\\ {\left({\overline{p}}_{n+1}^{\left(k\right)}-{\xi }_{p,n+1}^{\left(k\right)}\right)}^2-\left({\overline{p}}_{n+1}^{\left(k\right)}-{\xi }_{p,n+1}^{\left(k\right)}\right)p_{c,n+1}^{\left(k\right)}+\frac{3}{2\left(M^2-{\alpha }_0^2\right)}{\widehat{s}}_{ij,n+1}^{(k)}{\widehat{s}}_{ij,n+1}^{(k)}\\ {\overline{K}}_{p,n+1}^{\left(k\right)}-\left({\overline{p}}_{n+1}^{\left(k\right)}-{\xi }_{p,n+1}^{\left(k\right)}\right)p_{c,n+1}^{\left(k\right)}{\overline{n}}_{p,n+1}^{\left(k\right)}\left({\chi }_0-\frac{2\beta }{\left(M^2-{\alpha }_0^2\right)}\frac{{\overline{q}}_{n+1}^{\left(k\right)}}{{\overline{n}}_{p,n+1}^{\left(k\right)}}\right)\end{array},$$

(24)

where l is the number of nonlinear equations. Variable A can be presented as:

$$\begin{array}{l}A=p_c\left(p-o_p\right)+2\left(p-o_p\right)\left({\xi }_p-o_p\right)+\frac{3}{\left(M^2-{\alpha }_0^2\right)}\left[\left(s_{ij}-o_{ij}\right)-\left(p-o_p\right){\alpha }_{ij}^0\right]\\ \left[\left({\xi }_{ij}-o_{ij}\right)-\left({\xi }_p-o_p\right){\alpha }_{ij}^0\right]\end{array}$$

(25)

If the $\Vert r_{l,n+1}^{\left(k\right)}\Vert <$ tolerance (taken as ${10}^{-8}$), THEN EXIT

Else GO TO step 4

4)

Solve the linear equations:

$${\left(\frac{\partial r}{\partial U}\right)}_{n+1}^{(k)}\delta U_{n+1}^{(k)}=-r_{l,n+1}^{(k)},$$

(26)

where $\delta U=\left\{\delta p,\delta s_{ij},\delta \Delta {\epsilon }_v^p,\delta \Delta e_{ij}^p,\delta p_c,\delta b,\delta \Lambda ,\delta {\overline{K}}_p,\delta \omega \right\}$.

5)

Update stresses and internal variable $U_{n+1}^{(k+1)}=U_{n+1}^{(k)}+\delta U_{n+1}^{(k)}$, $k=k+1$ and GO TO step (3).

6)

Satisfy the convergence condition, END.

2.2. Numerical Scheme

Based on the proposed model, simulations were carried out with the finite element software ABAQUS (Version 6.14) (Dassault Systemes, Paris, France). The 3D finite element model for the analysis of a suction anchor embedded in a soft clay seabed was established to predict the accumulation of pore water pressure under cyclic loads. The uniform effective unit weight of the seabed was 7.17 kN/m³, representing a typical average value for seabed conditions. The parameters for seabed soil are listed in

Table 1. Properties of seabed adopted in the case studies.

Effective Unit Weight (kN/m3)	Compression Index (λ)	Swelling Index (κ)	Poisson’s Ratio (ν)	Permeability Coefficient (m/s)	$\mathbf{Void} \mathbf{Ratio} \left({\mathit{e}}_0\right)$
7.17	0.17	0.34	0.3	$1\times {10}^{-9}$	0.62
$\gamma$	${\varsigma }_r$	$\eta$	β
1.72	3.5	120	0.5

for clarity, in accordance with Dafalias [28], which were simulated by Manzari et al. [23] and Tao and Messiner [29]. The anchor length to diameter ratio was h/D = 3, where D is the diameter of the anchor. The suction anchor analyzed was 2 m in diameter, and the penetration depth was 6 m below the seabed floor. The suction anchor was expected to be anchored in a normally consolidated clay seabed. To reduce the influence of the boundary conditions, a seabed diameter of 10 times the suction anchor diameter and a seabed depth of 3.3 times the anchor wall was used. This model was built in terms of a half portion of the anchor and seabed, with the symmetry axis through the center of the anchor.

2.3. Meshing and Boundary Conditions

A graded mesh was used in the simulation. Denser meshes and sparser meshes were adopted in the regions close to the anchor and further away from the anchor, respectively. The finite element mesh used for the analysis of the anchor is illustrated in Figure 4. The suction anchor was modeled as rigid bodies, since the anchor is much stiffer than soft clay. The cyclic load amplitude $F_{cz}$ was 20 kN, with a typical cyclic period (T) of 10 s. The total calculation time was one hour (360T). A drainage boundary at the seabed surface was assumed, that is to say, the upper surface of the seabed layer was allowed to drain freely (i.e., p = 0 at z = 0). The displacements were fixed horizontally on the periphery and in both directions at the bottom of the model domain. The conventional Coulomb friction law was adopted to simulate the friction at the interface between the anchor wall and the surrounding soil. The coefficient of wall friction ($tan\delta$) was set to 0.42 [30].

3. Verification of the Model

In this section, we describe the calibration of the proposed bounding surface model with the combined isotropic-kinematic hardening rule against available existing published experimental data, including those of Stipho [31], Tao et al. [29], and Zhong et al. [32]. The triaxial tests included monotonic loading tests on anisotropically consolidated clay, pore pressure response tests under cyclic loads, and tests on the plastic strain accumulation of soft clay under cyclic loads, respectively. The validity of the proposed model was evaluated based on comparisons between numerical results and experiments. The properties of the soil adopted in this study were the same as those used in the experiments, as tabulated in

Table 2. Model parameters in verification cases.

Parameters	Kaolin Clay (Stipho)	Kaolin Clay (Tao)	Soft Clay (Zhong)
Slope of critical state line (M)	1.12	1.1	1.15
Compression index (λ)	0.14	0.17	0.25
Swelling index (κ)	0.05	0.34	0.05
Effective Poisson’s ratio (ν)	0.2	0.3	0.25
$\gamma$	2	1.72	1.5
${\varsigma }_r$	-	3.5	6
$\eta$	-	120	40
β	-	0.5	1.5

.

The first validation case involved monotonic loading tests on anisotropically consolidated clay, which was performed by Stipho [31]. The tests were performed under strain-controlled conditions. The axial strain was applied to the top surface of the specimen with a magnitude of 12% for compression and extension, respectively. The specimen was consolidated under $p_0=204$ kPa, with the initial void ratio of $e_0=1.1$, with the value of $K_0$ being about 0.8. Figure 5a,b show the normalized $q/p_0$ versus the axial strain ${\epsilon }_1$ and the normalized excess pore water pressure $u/p_0$ versus the axial strain ${\epsilon }_1$, respectively. Figure 5c illustrates a comparison between model simulation and experimental data in terms of stress paths in the normalized $q/p_0$ versus $p/p_0$. As shown in Figure 5c, when it reaches the critical state line(CSL), p from the compression test is different from that obtained by extension test. That is to say, the proposed model can capture a non-unique critical state line (CSL) in the $e-lnp$ space. As one can see, there is an excellent agreement between the model simulation and the experiment.

The second validation case involved cyclic stress-controlled loading tests on isotropically consolidated clay. Tao et al. [29] carried out a series of cyclic three-axis undrained shearing tests on normally consolidated clays, with an initial void ratio of $e_0=0.62$. The triaxial test loading frequency was 0.1 Hz. Figure 6a presents the experimental results and model predictions for a confining pressure of $p_0=450$ kPa, and a cyclic stress amplitude of $q_d=116$ kPa, and Figure 6b presents the experimental results and model predictions for a confining pressure of $p_0=350$ kPa, and a cyclic stress amplitude of $q_d=130$ kPa. It is found that the pore water pressure versus the number of cycles can be predicted well by the proposed model.

The third validation case involved plastic strain accumulation tests on isotropically consolidated clay, which were performed by Zhong et al. [32]. The tests were conducted under stress-controlled conditions. The axial stress was applied to the top surface of the specimen with different amplitudes $q_d$. The specimen was consolidated under an initial pressure of $p_0=50$ kPa, with an initial void ratio of $e_0=1.099$. The triaxial test loading frequency was 0.1 Hz. Figure 7 shows a comparison between the model predictions and the experimental data. Though the permanent strain predicted by the proposed model was smaller than that shown in the experiment at a higher stress level, the general trend was consistent. Figure 8 shows the effective stress path under typical conditions. With the increase in the cyclic number, the soft clay specimen finally reached a cyclic steady state, which means that the cyclic shakedown phenomena occur.

Overall, these validation cases demonstrate that cyclic behaviors of soft clay can be well predicted by the present model under cyclic loading conditions. The model is able to capture the build-up of pore water pressure caused by cyclic loads.

4. Numerical Results and Interpretations

4.1. Accumulation of Pore Water Pressure Around the Suction Anchor

Figure 9 indicates the residual pore pressure ($p_s$) distribution within the seabed around a suction anchor at different times. The maximum excess pore water pressure appeared at the base of the suction anchor. The excess pore water pressure accumulated primarily on the outside of the suction anchor, while the negative pore water pressure built up mainly on the inside. With the increase in loading time, the negative pore water pressure inside the suction anchor gradually developed upwards.

Figure 10 indicates the vertical distribution of the residual pore pressure within the seabed along the anchor wall. As can be seen from the figures, with an increase in the loading time, the residual pore pressures on both sides of the anchor wall increased remarkably. The maximum residual pore pressure was distributed in the region of 0.9 < z/h <1. The maximum value of excess pore water pressure was about 3 kPa, and the maximum value of negative pore water pressure was about 4 kPa. The maximum value of residual pore pressure for both sides appeared in the deep layer of the seabed, due to the larger effective stress acting on the anchor wall. It is worth noting that, in this region, the length of the suction anchor accounted for 10% of the total length, but the excess pore water pressure accounted for 27~31% of the total excess pore water pressure, and the negative pore water pressure accounted for 15.6~22% of the total negative pore water pressure. The external wall friction force was reduced, according to the principle of effective stress. Therefore, the uplift capacity of the suction anchor was diminished. The residual pore pressure at the seabed surface was zero. This is because, during cyclic loading, drainage occurs at the seabed surface.

Figure 11 illustrates the time series of residual pore pressure at position z/h = 0.95 within the seabed, where the maximum residual pore pressure occurred. Figure 11a,b show the accumulation of pore water pressure on the inside and outside of the suction anchor, respectively. The values of residual pore pressures increased continuously over time. On the outside of the anchor, the rate of increase in excess pore water pressure was rapid in the initial loading stage. Thereafter, the rate became slower with time. It is noted that the excess pore water pressure did not reach a steady state. On the inside of the anchor, the negative pore water pressure experienced a rapid increase in the initial loading stage, and then gradually reached a relatively steady state. A potential reason that might account for this phenomenon is that the soil inside the suction anchor showed shear dilatation characteristics under cyclic loads. When the plastic deformation reached a certain level, plastic deformation no longer increased due to the limitation of the anchor wall.

The method used to determine the uplift capacity was adopted from the limit equilibrium method for a suction anchor in clay recommended by API. For a fully sealed anchor, the failure mechanism of the suction anchor is usually assumed that the soil plug is broken and the soil plug is pulled out together with the anchor. This failure mechanism has been verified by tests for anchors with aspect ratios between 1.7 and 10 [33,34]. The uplift capacity can be expressed as:

$$R=Q_f+Q_P+W,$$

(27)

where W and $Q_P$ are the soil plug weight and reverse end bearing, respectively. $Q_f$ is the external wall friction force, which can be calculated as:

$$Q_f={\displaystyle {\int }_0^h{{\sigma }^{\prime }}_v}dz\left(Ktan\delta \right)\left(\pi D\right),$$

(28)

in which h represents the embedded depth, K is the lateral pressure coefficient, ${{\sigma }^{\prime }}_v$ denotes the effective vertical stress in soil, and D represents the external diameter of the anchor.

It is acknowledged that the excess pore water pressure in the seabed soil could diminish the effective stress in soil, consequently, reducing the friction along the external wall–soil interface. By assuming that K and δ remain as constants, the interface friction can be calculated according to Equation (28), and normalized by the initial values.

Figure 12 illustrates the relationship between the normalized friction force versus the loading time at the soil–anchor interface. It is shown that the increase in excess pore water pressure led to a slight decrease in the friction force. Under current conditions, the external wall–soil friction is reduced by about 4%.

4.2. Effect of Load Amplitude

To investigate the effect of the load amplitude, four grades of load amplitudes (F_CZ = 10, 20, 30, and 40 kN) were adopted in the simulation. The period of cyclic loads was 10 s. The total calculation time was one hour (360T). The seabed parameters remained constant. The residual pore pressure distributions on both sides of the anchor wall for various load amplitudes are plotted in Figure 13. It is shown that the load amplitude has a remarkable influence on the distribution characteristics of the pore water pressure. The increase in the load amplitude resulted in a remarkable increase in the residual pore pressure. In the region of 0.9 < z/h <1, the excess pore water pressure accounted for 17~31% of the total excess pore water pressure. As the load amplitude increased, the proportion increased gradually.

4.3. Effect of λ and κ

The soil properties are key influential factors that determine the development of residual pore pressure within the seabed. Among those soil properties, the compression index (λ) and swelling index (κ) have important effects on the distribution of the residual pore pressure [35]. To evaluate the effect of λ and κ on the distributions of residual pore pressure, κ was kept constant, and λ was changed. Analyses were performed with κ = 0.034 and κ/λ ranging from 0.2 to 0.4 [36]. Figure 14 shows the residual pore pressure distributions along the anchor wall when κ was 0.034 and λ was varied from 0.17 to 0.85. It is shown that κ/λ and the load amplitude had opposite effects on the residual pore pressure distribution. As κ/λ increases, the residual pore pressures on both sides of the anchor wall showed obvious decreases. In the region of 0.9 < z/h <1, the proportion of excess pore water pressure to the total excess pore water pressure increased from 20.1% at κ/λ = 0.4 to 31.3% at κ/λ = 0.2. The negative pore pressures inside of the anchor accounting for the total negative pore water pressures increases from 6.9% to15.6%.

5. Perforated Suction Anchor

5.1. New Anchor Structure Style

The numerical analysis shows that the excess pore water pressure appears outside of the anchor, while negative pore water pressure appears inside. According to the principle of effective stress, the soft clay structure on the external anchor wall is degraded followed by the degradation in strength. The uplift capacity of suction anchor gradually decreases with the development of excess pore water pressure. According to the distribution characteristics of the residual pore pressure, a perforated anchor is proposed based on the traditional suction anchor, as shown schematically in Figure 15. Perforation at the lower part of the anchor wall is helpful to reduce the accumulation of excess pore water pressure. The installation process of the suction anchor is divided into two stages. First, the suction anchor is installed by self-weight to around half the design depth, and then negative pressure is applied until the anchor reaches the final depth [37,38]. Perforation is located at the lower part of the suction anchor, and therefore, the new structure style will not influence the installation process.

5.2. Comparative Study of the Pore Water Pressure Distribution

To investigate the effectiveness of the perforated suction anchor on reducing the development of excess pore water pressure, the distribution characteristics of residual pore pressure around the perforated anchor are analyzed, and a comparison with the traditional anchor is presented. The perforation height is 0.5m, and this is located at the lower part of the suction anchor (0.9 < z/h <1).

Figure 16 shows a comparison of the residual pore pressure distributions between the traditional anchor and the perforated anchor. As shown in Figure 16, the accumulation of excess pore pressure outside of the perforated anchor is reduced remarkably, compared with the traditional structure. The excess pore water pressure around the perforated anchor is reduced by about 27.5%. Therefore, the perforated anchor can effectively reduce the accumulation of excess pore water pressure and increase the effective stress in the soil. In this paper, the perforation height and perforation size were not quantified according to the length ratio of suction anchor and the initial penetration depth, so further research is wanted.

6. Concluding Remarks

In this study, numerical simulations were conducted to investigate the buildup of pore water pressures in a soft clay seabed around a suction anchor under cyclic loading conditions. Based on the existing two-dimensional bounding surface model, a three-dimensional bounding surface model was proposed to describe the cyclic behaviors of soft clay. Based on the analysis results obtained, the most significant conclusions can be summarized as follows:

(1)

A damage-dependent bounding surface model with combined isotropic-kinematic hardening rule was proposed to predict the accumulation of pore water pressure. The proposed model reasonably agreed well with the experiment result, against triaxial tests on anisotropically consolidated saturated clays and normally consolidated saturated clays. Thus the presented model was available to describe the key features of the cyclic behaviors of soft clay under cyclic loading conditions, including the pore water pressure response, accumulation of plastic deformation, and initial anisotropy.

(2)

Under the vertical cyclic loading condition, the excess pore water pressure primarily appeared on the outside of the suction anchor, and negative pore pressure mainly appeared on the inside, respectively. The maximum values on both sides appeared in the lower part of the seabed, due to the larger effective stress in the seabed soil, and these increased gradually with the loading time. The accumulation of excess pore water pressure can decrease the effective stress in the soil, and further reducing the uplift capacity of the suction anchor.

(3)

According to the distribution characteristics of the residual pore pressure around the suction anchor, a new structure can reduce the accumulation of excess pore water pressure in the lower part of the seabed was proposed, which could increase external wall friction force under cyclic loading conditions. For the working conditions adopted in the present study, the new structure can reduce the excess pore water pressure by about 27.5%.