Fractional-Order Elastoplastic Modeling of Sands Considering Cyclic Mobility

Abstract

Seabed soil may experience a reduction in strength or even liquefaction when subjected to cyclic loadings exerted by offshore structures and environmental loadings such as ocean waves and earthquakes. A reasonable and robust constitutive soil model is indispensable for accurate assessment of such structure–seabed interactions in marine environments. In this paper, a new constitutive model is proposed by enriching subloading surface theory with a fractional-order plastic flow rule and multiple hardening rules. A detailed validation of both stress- and strain-controlled undrained cyclic test results of medium-dense Karlsruhe fine sand is provided to demonstrate the robustness of the present constitutive model to capture the non-associativity and cyclic mobility of sandy soils. The new fractional cyclic model is then implemented into a finite element code based on a two-phase field theory via a user subroutine, and a numerical case study on the response of seabed soils around a submarine pipeline under cyclic wave loadings is presented to highlight the practical applications of this model in structure–seabed interactions.

1. Introduction

In recent decades, considerable efforts have been devoted to the study of structure–seabed interactions in marine environments [1]. Numerical modeling [2,3,4], physical modeling [5,6,7], and field observation [8,9] are all regarded as effective methods to investigate such complicated problems and are adopted widely by researchers. Among these, when a large number of various conditions are of interest, numerical solutions are preferred owing to their flexibility and low cost and, thus, have been continuously pursued. In general, the seabed is subjected to cyclic loadings induced by both offshore structures built in/on it and marine environment phenomena such as ocean waves and earthquakes [10]. Excess pore pressure and plastic strain would accumulate in the seabed under these cyclic loadings, and a reduction in or even loss of strength (i.e., liquefaction) is expected to take place for the seabed soils. Therefore, although great success has been achieved by considering the seabed as elastic in the simulation of structure–seabed interactions in marine environments [11,12], the development of elasto-plastic constitutive models that can reasonably mimic the behavior of soils under complex cyclic loadings has attracted enthusiastic interest recently [13,14,15,16], especially for sandy soils whose liquefaction behavior under cyclic loading is considered to have a great effect on the stability and serviceability of marine structures.

Numerous constitutive models based on various plastic theories [17,18,19,20] have been proposed to represent drained and undrained cyclic behaviors of sandy soils. Dafalias [17] introduced the well-known bounding surface plasticity model as an efficient tool for capturing the rate-independent mechanical behavior of geomaterials. It was widely accepted and gradually enriched by many other researchers [21,22,23,24] to replicate the static and cyclic behaviors of both loose and dense sandy soils. Pastor [18] developed a simple framework of generalized plasticity in which no plastic potential, yield surfaces, or consistency rules need to be explicitly defined to consider the liquefaction and cyclic mobility behavior of sands under wave or earthquake loading. Due to that, some of the parameters in generalized plasticity models may not have exact physical meanings. Ling and Liu [25] extended the generalized plasticity model with fifteen parameters to include pressure dependency and densification behavior of sand. Wu et al. [19] incorporated the critical state concept into the hypoplastic model and used only four parameters to simulate the behavior of sands with different densities, which was extended by Liao and Yang [26] to describe the fabric effect of sand under monotonic and cyclic loading conditions with the anisotropic critical state theory. Zhang et al. [27] proposed a cyclic mobility model considering the evolution of stress-induced anisotropy, over-consolidation ratio, and structure within the framework of the subloading surface theory [20] to describe the cyclic mobility behavior of medium-dense sand. The cyclic mobility model is able to describe the behavior of sands with different densities using one same set of material parameters and mimic the cyclic mobility (e.g., butterfly-shaped stress path) of sands upon liquefaction under cyclic loadings. In spite of these abundant constitutive models based on different plastic theories, further model evolution for sands under cyclic loadings is still required to reasonably consider the non-associated, state-dependent (pressure, density) behaviors and the cyclic mobility under undrained relaxation of mean effective stress. Recently, a novel fractional plasticity approach emerged and has been gradually adopted to consider the mechanical behaviors of sand and clay [28,29,30]. Compared to traditional bounding surface models, an additional plastic potential surface is no longer required and the fractional gradient at the current yield surface is used to determine the plastic flow direction. However, most current fractional plasticity models are limited to modeling the monotonic behaviors of soils, and only a few of them have incorporated the cyclic behaviors of soils [30].

In this paper, a modified subloading model for sands with a new fractional-order plastic flow rule is proposed by which the state dependency, non-associativity, and cyclic mobility behavior of sands can be well captured. The experimental database of Karlsruhe fine sand under two-way stress- and strain-controlled cyclic loading is used to validate the performance of the fractional cyclic model at the element level. Following that, the new fractional cyclic model, coded into a user subroutine, is implemented into a finite element code based on a two-phase field theory [31], and the response of soils around a submarine pipeline under cyclic wave loadings is studied through large-scale numerical simulation.

2. New Fractional Cyclic Model

In this section, a detailed description of the newly proposed fractional cyclic model is presented, including the elastic strain part, the subloading yield surface, the fractional-order plastic flow rule, multiple hardening rules, and the incremental stress–strain relationship.

2.1. Elastic Strain

Based on the elastoplastic theory, the total strain ${\epsilon }_{ij}$ and its increment ${d\epsilon }_{ij}$ are decomposed into the elastic strain component and the irreversible plastic strain component:

$${\epsilon }_{ij}={\epsilon }_{ij}^p+{\epsilon }_{ij}^p,{d\epsilon }_{ij}{=d\epsilon }_{ij}^e{+d\epsilon }_{ij}^p$$

(1)

For the elastic strain part in Equation (1), Hooke’s law is used:

$${d\epsilon }_{ij}^e=\frac{1+v}{E}d{\sigma }_{ij}-\frac{v}{E}d{\sigma }_{kk}{\delta }_{ij}$$

(2)

where v denotes Poisson’s ratio under drained conditions, E is Young’s modulus, and ${\delta }_{ij}$ represents the Kronecker delta.

Young’s modulus E and the bulk modulus K are determined by the swelling line in the e-lnp′ plane under isotropic loading/unloading conditions:

$$E=\frac{K}{3\left(1-2v\right)},K=\frac{1{+e}_0}{\kappa }p\prime$$

(3)

where κ is the slope of the swelling line in the e-lnp′ plane, e₀ is the initial void ratio, and $p\prime$ is the current mean effective stress.

2.2. Subloading Yield Surface

Two state variables, namely the similarity ratio of the superloading yield surface to the reference yield surface R* and the similarity ratio of the superloading yield surface to subloading yield surface R, are defined. A brief description of three yield surfaces is shown in Figure 1.

$$R^{*}=\frac{\hat{p}}{\overline{p}}=\frac{\hat{q}}{\overline{q}},R=\frac{p}{\overline{p}}=\frac{q}{\overline{q}}$$

(4)

where $\left(p\prime ,q\right)$, $\left(\hat{p}\prime ,\hat{q}\right)$, and $\left(\overline{p}\prime ,\overline{q}\right)$ represent the current stress, the reference stress, and the structured stress, respectively.

The current subloading surface is given in the following two equal forms:

$$f=\ln \frac{p\prime }{{p\prime }_0}+\ln \frac{M^2-{\zeta }^2{+\eta }^{\ast }{}^2}{M^2-{\zeta }^2}+\ln R^{*}-\ln R-\frac{{\epsilon }_v^p}{C_p}=0$$

(5)

$$\overline{f}=S^2+\left(M^2-{\zeta }^2\right)p\prime \left(p\prime -\frac{{Rp\prime }_0}{R^{*}}\exp \left(\frac{{\epsilon }_v^p}{C_p}\right)\right)=0$$

(6)

where $p\prime$ is the current mean effective stress and ${p\prime }_0$ = 98 kPa is the reference mean effective stress; ${\beta }_{ij}$ is the anisotropic stress tensor, $\zeta =\sqrt{\frac{3}{2}{\beta }_{ij}{\beta }_{ij}}$ is the anisotropic state variable, and ${\eta }^{\ast }=\sqrt{\frac{3}{2}\left(\frac{{\sigma }_{ij}}{p\prime }{-\delta }_{ij}{-\beta }_{ij}\right)\left(\frac{{\sigma }_{ij}}{p\prime }{-\delta }_{ij}{-\beta }_{ij}\right)}$ is the anisotropic stress state difference between the stress ratio tensor ${\eta }_{ij}=\frac{{\sigma }_{ij}}{p\prime }{-\delta }_{ij}$ and ${\beta }_{ij}$; ${\epsilon }_v^p$ is the volumetric plastic strain; $C_p=\frac{\lambda -\kappa }{1{+e}_0}$ is the dilatancy parameter which can be expressed by $\lambda$ and $\kappa$, the compression and the swelling index, respectively; and $S=\sqrt{\frac{3}{2}{S}_{ij}{S}_{ij}}$ is the deviatoric stress with the deviatoric stress tensor $S_{ij}={\sigma }_{ij}-{p\prime \delta }_{ij}-{p\prime \beta }_{ij}$. It should be noted that the current stress point always lies on the subloading yield surface by definition.

2.3. Fractional-Order Plastic Flow Rule

In order to better replicate the non-associated behavior of sandy soils, a fractional-order plastic flow rule based on the subloading yield surface is adopted in the proposed constitutive model, which can be expressed as:

$${d\epsilon }_{ij}^p=\Lambda \frac{{\partial }^{\mu }\overline{f}}{\partial {\sigma }_{ij}^{\mu }}=\Lambda \frac{{\partial }^{\mu }\overline{f}}{\partial S_{kl}^{\mu }}\frac{\partial S_{kl}}{\partial {\sigma }_{ij}}$$

(7)

where $\frac{{\partial }^{\mu }\overline{f}}{\partial S_{kl}^{\mu }}$ is the fractional derivative of the subloading yield surface at the point $S_{kl}$.

Since the subloading yield surface function $\overline{f}$ is made up of two stress invariants $S$ and $p\prime$, the fractional-order plastic flow rule can be obtained:

$${d\epsilon }_{ij}^p=\Lambda \frac{{\partial }^{\mu }\overline{f}}{\partial {\sigma }_{ij}^{\mu }}=\Lambda \left(\frac{{\partial }^{\mu }\overline{f}}{\partial S_{}^{\mu }}\frac{\partial S}{\partial {\sigma }_{ij}}+\frac{{\partial }^{\mu }\overline{f}}{\partial {p\prime }_{}^{\mu }}\frac{\partial p\prime }{\partial {\sigma }_{ij}}\right)$$

(8)

The plastic flow direction is highly relevant to the definition of the adopted fraction derivative and the Caputo fraction derivative operator is selected here for simplicity, which is given in the following form:

$$\left(\frac{{\partial }^{\mu }\overline{f}}{\partial {p^{\prime }}^{\mu }},\frac{{\partial }^{\mu }\overline{f}}{\partial S^{\mu }}\right)=\left(\frac{\mu \left(M^2-{\alpha }^2\right){p^{\prime }}^{2-\mu }-S^2\left(2-\mu \right){p^{\prime }}^{-\mu }}{\Gamma \left(3-\mu \right)},\frac{2S^{2-\mu }}{\Gamma \left(3-\mu \right)}\right)$$

(9)

with the famous Euler Gamma function $\Gamma \left(z\right)={\displaystyle {\int }_0^{\infty }{e}^{-\tau }}{\tau }^{z-1}d\tau ,(Re\left(z\right)>1)$.

Thus, the fractional order $\mu$ controls the plastic flow direction by determining the ratio of volumetric plastic strain increment to deviatoric plastic strain increment. The larger the fractional order $\mu$ is, the larger the dilatancy ratio is (the ratio of volumetric plastic strain increment to deviatoric plastic strain increment). It should be noted that the value of fractional order $\mu$ is between 0 and 2. If the fractional order $\mu$ is equal to 1, the plastic flow rule would turn into the associated flow rule and the model would degenerate into the original cyclic model by Zhang et al. [27].

2.4. Hardening Rules

In order to consider the cyclic mobility of sands, multiple hardening rules are defined in the proposed model, namely the evolution rule for the anisotropic stress tensor, the evolution rule for the degree of structure, and the changing rate of over-consolidation ratio.

The evolution rule for the anisotropic stress tensor is defined as:

$$d{\beta }_{ij}=\frac{M}{C_p}{b}_r(b_{l}M-\zeta )\sqrt{\frac{3}{2}}d{\epsilon }_d^p\cdot \frac{{\eta }_{ij}-{\beta }_{ij}}{\eta \ast }$$

(10)

where $b_r$ is a parameter that controls the changing rate, $b_l$ is a fixed constant of 0.95, suggested by Zhang et al. [27], and $d{\epsilon }_d^p$ is the deviatoric plastic strain increment. It could be found that when ${\eta }_{ij}$ is equal to ${\beta }_{ij}$, the development of anisotropy would stop.

The hardening rule for the degree of structure $R^{\ast }$ is defined as follows:

$$d{R}^{\ast }=\frac{aM}{C_p}{R}^{\ast }(1-R^{\ast })d{\epsilon }_d^p$$

(11)

where $a$ is the parameter that controls the collapsing rate of soil’s structure during shearing. The value of structure degree $R^{\ast }$ is between 0 and 1 according to the definition. When the value of structure degree $R^{\ast }$ is equal to 1, the evolution rule would have no effect on the hardening process. At the same time, the value of structure degree would increase monotonically upon shearing.

The changing rate of over-consolidation ratio is controlled by two factors, namely the plastic component of stretching and the increment in anisotropy. It is suggested by Zhang et al. [27] as follows:

$$dR=-\frac{mM}{C_p}\left(\frac{{\left(p^{\prime }/{p^{\prime }}_0\right)}^2}{{\left(p^{\prime }/{p^{\prime }}_0\right)}^2+1}\right)\ln R\Vert d{\epsilon }^p\Vert +\frac{R}{C_p}\frac{\eta }{M}\frac{\partial f}{\partial {\beta }_{ij}}d{\beta }_{ij}$$

(12)

$$\Vert d{\epsilon }^p\Vert =\sqrt{d{\epsilon }_{{}^{ij}}^{p}d{\epsilon }_{{}^{ij}}^p}=\Lambda \sqrt{\frac{{\partial }^{\mu }\overline{f}}{\partial {\sigma }_{ij}^{\mu }}\frac{{\partial }^{\mu }\overline{f}}{\partial {\sigma }_{ij}^{\mu }}}$$

(13)

where $\eta =\sqrt{\frac{3}{2}\left(\frac{{\sigma }_{ij}}{p\prime }-{\delta }_{ij}\right)\left(\frac{{\sigma }_{ij}}{p\prime }-{\delta }_{ij}\right)}$ is the stress ratio.

2.5. Incremental Stress–Strain Relationship

Within the framework of the plasticity theory, the consistency condition of the subloading yield surface is adopted herein:

$$\frac{\partial f}{\partial {\sigma }_{ij}}{d\sigma }_{ij}+\frac{\partial f}{\partial {\beta }_{ij}}{d\beta }_{ij}+\frac{1}{R^{*}}{dR}^{*}-\frac{1}{R}dR-\frac{1}{C_p}{d\epsilon }_v^p=0$$

(14)

By submitting the multiple hardening rules and fractional-order plastic flow rule into Equation (14), the plastic multiplier $\Lambda$ can be obtained:

$$\Lambda =-\frac{\frac{\partial f}{\partial {\sigma }_{ij}}{d\sigma }_{ij}}{\frac{\partial f}{\partial {\beta }_{ij}}\frac{\partial {\beta }_{ij}}{\partial S}\frac{{\partial }^{\mu }\overline{f}}{\partial S_{}^{\mu }}+\frac{1}{R^{*}}\frac{\partial R^{*}}{\partial S}\frac{{\partial }^{\mu }\overline{f}}{\partial S_{}^{\mu }}-\frac{1}{R}\left(\frac{\partial R}{\partial S}\frac{{\partial }^{\mu }\overline{f}}{\partial S_{}^{\mu }}+\frac{\partial R}{\partial p\prime }\frac{{\partial }^{\mu }\overline{f}}{\partial {p\prime }_{}^{\mu }}\right)-\frac{1}{C_p}\frac{{\partial }^{\mu }\overline{f}}{\partial {p\prime }_{}^{\mu }}}$$

(15)

3. Validation

There are nine independent material parameters in the newly proposed fractional cyclic model, i.e., $\lambda$, $\kappa$, M, N, ν, m, a, b_r, and µ. The five parameters $\lambda$, $\kappa$, M, N, and ν are the same parameters of the modified Cam-Clay model. m controls the degradation of the over-consolidation ratio and can be obtained by isotropic compression test. a and b_r control the degradation of structure and the evolution rate of anisotropy. They can be obtained by undrained cyclic tests. µ is the fractional order that determines the plastic flow direction and can be obtained by the dilatancy ratio.

In this section, two-way stress- and strain-controlled undrained cyclic test results of Karlsruhe fine sand are simulated using the proposed model to show its validity and performance in capturing the state dependency, non-associativity, and cyclic mobility of sands. More details on the test material and test setup can be found in previous studies [32,33]. All the fractional cyclic model parameters of Karlsruhe fine sand are calibrated and listed in

Table 1. Material parameters of Karlsruhe fine sand.

Properties	Values
Compression index λ	0.050
Swelling index κ	0.012
Shear stress ratio at critical state M	1.33
Void ratio N (p = 98 kPa on N.C.L.)	1.103
Poisson’s ratio ν	0.05
Degradation parameter of over-consolidation state m	0.1
Degradation parameter of structure a	5
Evolution parameter of anisotropy br	5
Fractional order µ	0.9/1/1.1/1.2/1.3

.

3.1. Stress-Controlled Cyclic Test

The undrained stress-controlled cyclic test results of the medium-dense sand sample with the initial void ratio (e = 0.83) and initial isotropic pressure (p₀ = 300 kPa) were selected to validate the feasibility of the proposed constitutive approach. The stress cycles were applied in an amplitude/pressure ratio (q^ampl/p₀) of about 0.3, using a constant displacement rate of 5 mm/min. A typical result on the medium-dense sand sample is shown in Figure 2, in which the stress path in the mean effective stress-deviatoric stress plane and the deviatoric stress versus axial strain are presented.

As shown in Figure 2a, the stress path starts from the initial stress point (300 kPa, 0) and then turns into an increasing contractive phase. As the effective stress point approaches zero, the stress path repeatedly exhibits a butterfly loop, which is often referred to as cyclic mobility [34]. That is, in each loop, the stress path turns from a contractive phase to a dilative phase alternately. This behavior is properly incorporated in the present model to capture the realistic soil response upon liquefaction. As shown in Figure 2b, the axial strain progressively accumulates with the increasing loading cycle. During the initial several loops, the accumulated axial strain is quite small (smaller than 0.5%). When the stress point approaches zero and the stress path turns into the butterfly loop, the axial strain increases rapidly and reaches around 4~6% after 15 loading cycles.

The proposed fractional model with different fractional orders µ = 1 (associated flow rule) and µ = 0.9, 1.1, 1.2, and 1.3 (non-associated flow rule) was used to calibrate the results in Figure 2. For the case of µ = 1.1, 1.2, and 1.3, a faster decrease in effective stress was observed in the mean effective stress-deviatoric stress plane compared with the predicted results of the original model with the associated flow rule, although the cyclic mobility is well described in all of these four cases upon liquefaction. However, an opposite trend was found for the case of µ = 0.9, where discernable effective stress still exists at the end of cyclic loading, indicating that the soil sample is not fully liquefied. A more marked difference in axial strain can be found in Figure 2b, where the maximum predicted axial strain in the case of µ = 1.3 (9%) is up to six times that in the case of µ = 0.9 (1.5%), suggesting the important effect of non-associativity on the deformation characteristics of sands subjected to cyclic loadings. Generally speaking, the simulated accumulated axial strain in the case of µ = 1.1 is much closer to the test data.

3.2. Strain-Controlled Cyclic Test

For the strain-controlled cyclic test, the test results of a medium-dense sand sample with the initial void ratio (e = 0.83) and initial isotropic pressure (p₀ = 200 kPa) were adopted. The stress cycles were applied in an axial strain range of about 0.5%, with a constant displacement rate of 5 mm/min.

The typical undrained cyclic test data of Karlsruhe fine sand are shown in Figure 3. Figure 3a presents the stress path that gradually approaches zero from the initial stress point. With the increasing loading cycles, the maximum deviatoric stress rapidly decreases and eventually oscillates with small amplitude around the liquefaction point (0, 0). As the effective stress point approaches zero, the stress path also exhibits a tapering butterfly loop. As shown in Figure 3b, after only 2~3 loading cycles, the deviatoric stress reduces to approximately zero, suggesting the occurrence of liquefaction and following complete loss of shear strength.

The test results were also calibrated using five different fractional orders, µ = 0.9, 1, 1.1, 1.2, and 1.3, similar to those in the case of the stress-controlled cyclic test. As shown in Figure 3, both cyclic mobility and decreasing deviatoric stress are reasonably captured by the proposed model with different fractional orders. Furthermore, with the increasing fractional order, the simulation results also show a faster decreasing rate of deviatoric stress, and again, the results in the case of µ = 1.1 are generally closer to the test data.

4. Application

In this section, the fractional cyclic constitutive model described above is adopted in a finite element analysis to investigate the liquefaction potential and dynamic response of seabed soils around a fully buried submarine pipeline under wave loadings, and the performance and characteristics of the model are further discussed to highlight its practical applications for the assessment of structure–seabed interactions in marine environments.

4.1. Model Setup

A schematic diagram of the plane strain numerical model with the configuration and dimensions illustrated is shown in Figure 4. The computational domain of the poro-elastoplastic seabed is 120 m in length and 10 m in thickness. A pipe section with a diameter D of 1 m is fully buried at a burial depth b of 2 m in a backfilled trench layer, whose bottom width, depth, and slope are 1 m, 3.5 m, and 1, respectively. Progressive waves are assumed to propagate over the seabed surface in the x-direction. The wave parameters used in the numerical simulation are as follows: water depth d = 10 m, wave period T = 5 s, wave height H = 1 m, and wave length L = 36.6 m, calculated using the wave dispersion equation. The seabed length is more than three times larger than the wave length, which should effectively eliminate the boundary effects according to Ye and Jeng [35].

The following boundary conditions were set: Firstly, the seabed was considered to be laid upon an impermeable rigid bottom, so no displacement and flow occur at the bottom boundary. Secondly, lateral boundaries were also set as impermeable, with the displacement in the horizontal direction fixed. Thirdly, the rigid pipeline was assumed impermeable and there was no relative displacement between the soil and the pipeline surface. Finally, no displacement constraint was set at the seabed surface, and the pore pressure $u_w$ at the seabed surface was equal to the wave-induced pressure calculated using second-order Stokes wave theory:

$$u_w=\frac{1}{2}{\gamma }_{w}H\left\{\frac{\cos \left(kx-\omega t\right)}{\cosh \left(kd\right)}+\frac{3\pi }{2}\left(\frac{H}{L}\right)\frac{\cos 2\left(kx-\omega t\right)}{\sinh \left(2kd\right)}\left[\frac{1}{{\sinh }^2\left(kd\right)}-\frac{1}{3}\right]\right\}$$

(16)

where ${\gamma }_w$ denotes the unit weight of water; $k$ and $\omega$ are the wave number (defined as $2\pi /L$) and the wave angular frequency (defined as $2\pi /T$), respectively. Here, the waves are simplified to a hydro-mechanical boundary at the seabed surface. Although this treatment has been widely adopted in previous related numerical simulations on wave–seabed–pipeline interactions (e.g., Chen et al. [16]; Qi et al. [36]) and deemed workable (e.g., Chen et al. [37]; Zhu et al. [38]), an advanced computational fluid dynamics (CFD) model and a wave-seabed coupled scheme can better capture the interaction between seabed and water at the interface and benefit the understanding of the problems investigated [39].

The proposed fractional cyclic model was implemented into an existing finite element code [40] via a user subroutine to solve the boundary value problem defined above. Based on a two-phase field theory within the framework of Biot formulation, this hydro-mechanical model is able to execute both static and dynamic analyses for soil–structure interaction and more details of its proper functioning can be found in Zhu et al. [41] and Ye et al. [15]. Besides, the detailed validation of this finite element (FE) code against experimental results was given by Chen et al. [16] and, thus, is not provided here. The equilibrium equation and the continuity equation were adopted as the governing equations, which can be expressed as:

$${\rho }_s{\ddot{u}}_i^s=\frac{\partial {\sigma }_{ij}}{\partial x_j}+{\rho }_s{b}_i$$

(17)

$${\rho }_f{\ddot{\epsilon }}_{ii}^s-\frac{{\partial }^{2}p}{\partial x_i\partial x_i}-\frac{{\gamma }_f}{k_s}\left(\frac{{\dot{\epsilon }}_{ii}^s}{n}-\frac{1}{k_f}\dot{p}\right)=0$$

(18)

where ${\rho }_s$ and $b_i$ denote the soil density and body force, respectively; ${\gamma }_f$ is the density and $k_f$ is the specific gravity of the fluid. $k_s$, $n$, and $k_f$ represent the coefficient of permeability, the soil porosity, and the volumetric compressibility of the fluid, respectively. The same parameters for Karlsruhe fine sand as listed in Table 1 were used in calculations for consistency, except that only µ = 1, 1.1, and 1.2 were considered for brevity.

4.2. Homogeneous Seabed

In this section, the soils inside the trench layer are same as that in the seabed (i.e., fine sands); hence, the seabed considered is homogeneous. The soil density was set as 1700 kg/m³ and the permeability coefficient was set as 1 × 10⁻⁵ m/s.

4.2.1. Wave-Induced Liquefaction

The evolution of the wave-induced liquefaction area within the seabed at four different moments (t/T = 15, 30, 45, and 60) is shown in Figure 5. Note that the pipeline is under a wave crest at the four moments selected and an area of 20 m in length and 10 m in thickness with the pipeline at the center is presented. Two criteria were used to evaluate the degree of liquefaction: liquefaction was considered to take place when the excess pore pressure $u_w$ reached the initial overburden effective stress ${\sigma \prime }_{v0}$ [42]; alternatively, liquefaction was considered to take place when the generalized deviatoric strain ${\epsilon }_d$ reached the level of 5% [34].

The distribution of normalized excess pore pressure (defined as $u_w/{\sigma \prime }_{v0}$) is shown in Figure 5a. A progressive pattern was found for the development of liquefaction front, which is in accordance with the observations in previous physical and numerical modelling on wave-induced liquefaction [35,42]. At t/T = 30, the soils around the pipeline have been completely liquefied, and it seems that the presence of a pipeline moderately accelerates the accumulation of excess pore pressure around it. However, the pipeline may act as a “barrier” when the liquefaction front far from the pipeline moves to a slightly deeper depth than that just under the pipeline. It could be observed that after t/T = 45, the development rate of liquefaction around the pipeline was slightly slower than other areas in the same vertical plane.

Figure 5b shows the distribution of deviatoric strain ${\epsilon }_d$ at four different moments (t/T = 15, 30, 45, and 60). The development of liquefaction determined using the criterion of deviatoric strain was much slower than that using the criterion of excess pore pressure, especially at early stages. At t/T = 15, only seabed soils above the pipeline undergo discernable shearing with a deviatoric strain level of about 2.5%. From t/T = 30 to 60, the area of seabed that reaches an extreme deviatoric strain of 5% grows rapidly, with the exception of soils underneath the pipeline. While the liquefaction front extends to a depth of about 4 m within the region far away from the pipeline at t/T = 60, only the soils around the top half of the pipeline reach the critical deviatoric strain level (5%).

4.2.2. Excess Pore Pressure

The accumulation of wave-induced excess pore pressure within the seabed is closely related to the plastic volumetric contraction of sandy soils. Therefore, different flow rules adopted in the constitutive model would have an important influence on the simulation results. Time traces of excess pore pressure at two typical locations, namely the top (point A) and the bottom (point B) of the pipeline (see Figure 4), are further examined here, and results derived using both the associated flow rule (μ = 1) and non-associated flow rule (μ = 1.1, 1.2) are provided for comparison.

As shown Figure 6a, the excess pore pressure rises rapidly in the beginning stage (t < 100 s) and then reaches a plateau when it rises to the initial overburden effective stress, thus suggesting the occurrence of liquefaction. Similar trends can be observed for all three kinds of soils with μ = 1, 1.1, and 1.2, although a slightly higher accumulation rate of excess pore pressure is observed with the increasing fractional order. On the other hand, as shown in Figure 6b, at point B, where a submarine pipeline is lying above, excess pore pressure would hardly reach ${\sigma \prime }_{v0}$ due to the presence of the pipeline as a barrier. Meanwhile, the results with the fractional order μ = 1.1 and 1.2 exhibit a considerably higher accumulation rate of excess pore pressure during t = 50–150 s and reach a higher level after entering the plateau phase than that of μ = 1, indicating that the non-associated behavior of sandy soils may have an important effect on the build-up of wave-induced excess pore pressure subjected to certain conditions.

4.2.3. Effective Stress Path

Figure 7 illustrates the effective stress paths of point A and point B in the mean effective stress-deviatoric stress plane, in which the current stress state and loading history of soil are revealed and the strength and deformation of soil can be inferred. For point A, the initial stress point is around (8 kPa, 6 kPa) as shown in Figure 7a. The mean effective stress and deviatoric stress both decrease sharply with an oscillatory pattern as soon as the wave loadings are applied. As the mean effective stress drops continuously and approaches zero, the cyclic mobility behavior of soil is observed with the characteristic irregular butterfly loop of stress path around the liquefaction point (0, 0). The soils with μ = 1.1 and 1.2 enter the cyclic mobility phase prior to that with μ = 1 and show a lower amplitude of deviatoric stress during this phase. As shown in Figure 7b, the stress path of soil at point B does not reach the liquefaction point with a mean effective stress of about 0.5 kPa remaining, suggesting that complete liquefaction is not achieved. Meanwhile, the stress path is also affected by fractional order μ in the proposed constitutive model. Since a larger fractional order μ represents a larger dilatancy ratio, the stress path shows a similar trend to that in Figure 7a, i.e., the stress path of soil with μ = 1.1 and 1.2 moves more rapidly towards the liquefaction point and the corresponding butterfly loop is much smaller than that in the case of μ = 1.

4.3. Seabed with a Trench Layer

In this section, the seabed soils are fine sands, while the backfilling materials of the trench layer are considered to have relatively high permeability. The parameters for the backfilling materials were as follows: $\lambda$ = 0.06, $\kappa$ = 0.002, M = 1.45, N = 0.8, ν = 0.2, m = 0.1, a = 2.2, and b_r = 1.5. The fractional order µ for both seabed soils and backfilling materials was set as 1.1 to properly capture the non-associated behavior of granular materials. The permeability coefficients of the backfilling materials ranged from 1 × 10⁻⁴ to 1 × 10⁻³ m/s, higher than that of the fine sands (1 × 10⁻⁵ m/s).

As shown in Figure 8a, the wave-induced liquefaction zone reduces remarkably with the increase in permeability coefficient of backfilling materials due to the rapid dissipation of wave-induced excess pore pressure. When k_s is 1 × 10⁻³ m/s, the liquefaction zone merely reaches a depth of 1 m within the trench layer, indicating successful protection for the pipeline. The performance of the trench layer can be further verified by the time traces of excess pore pressure at point A. As shown in Figure 8b, the wave-induced excess pore pressure for the case of k_s = 1 × 10⁻³ reaches a maximum of approximately 6 kPa (50% of the initial vertical effective stress) at 200 s and decreases afterwards, suggesting that no liquefaction occurs around the pipeline.

5. Conclusions

In this paper, a fractional cyclic constitutive model is proposed by incorporating a novel fractional-order plastic flow rule. The proposed model is capable of capturing the state dependency, cyclic mobility, and non-associated behavior of sands, and the performance of it is validated by a comparison with the results of undrained cyclic triaxial tests on Karlsruhe fine sand. After implementation into a finite element program based on a two-phase field theory, this model was adopted in a numerical case study on wave-seabed-pipeline interaction. The main conclusions are drawn as follows:

• The non-associated flow rule was successfully incorporated into the proposed model based on the fractional-order plasticity theory, and the fractional order that determines the plastic flow direction can be obtained by the dilatancy ratio. Combined with the multiple hardening rules, the state dependency, cyclic mobility, and non-associated behavior of sands are reasonably mimicked by the proposed model.
• With the increase in the fractional order, the dilatancy ratio and the critical state ratio would both increase. In undrained cyclic conditions, a larger fractional order would result in a quicker accumulation of excess pore pressure in the soil sample modeled. Accordingly, the stress point would approach the liquefaction state more quickly.
• The proposed model shows good robustness during large-scale numerical simulation of dynamic geotechnical problems. Based on the results of numerical simulation, it was found that non-associativity of sand has an important effect on the accumulation of wave-induced excess pore pressure and plastic strain. Besides, soils at the top of the pipeline are more prone to wave-induced liquefaction than those at other locations within the seabed, and a trench layer of non-liquefiable materials with high permeability is found useful and is, thus, recommended to prevent submarine pipeline from seabed instability under wave actions.