Meshfree Model for Wave-Seabed Interactions Around Offshore Pipelines

Abstract

The evaluation of the wave-induced seabed instability around a submarine pipeline is particularly important for coastal engineers involved in the design of pipelines protection. Unlike previous studies, a meshfree model is developed to investigate the wave-induced soil response in the vicinity of a submarine pipeline. In the present model, Reynolds-Averaged Navier-Stokes (RANS) equations are employed to simulate the wave loading, while Biot’s consolidation equations are adopted to investigate the wave-induced soil response. Momentary liquefaction around an offshore pipeline in a trench is examined. Validation of the present seabed model was conducted by comparing with the analytical solution, experimental data, and numerical models available in the literature, which demonstrates the capacity of the present model. Based on the newly proposed model, a parametric study is carried out to investigate the influence of soil properties and wave characteristics for the soil response around the pipeline. The numerical results conclude that the liquefaction depth at the bottom of the pipeline increases with increasing water period (T) and wave height (H), but decreases as backfilled depth ($H_b$), degree of saturation ($S_r$) and soil permeability (K) increase.

1. Introduction

Offshore pipelines have been a commonly used facility for transportation of offshore oil and gas. In addition to construction causes, another key failure mode is the wave-induced seabed instability in the vicinity of pipelines [1,2]. Therefore, the evaluation of seabed stability around the pipeline is one of key factors that needs to be considered in an offshore pipeline project.

In general, ocean waves will exert fluctuations of dynamic pressures over the sea floor, which will further induce excess pore pressures and effective stresses within the seabed. The shear resistance in the vicinity of pipelines may be loss due to the liquefaction of surrounding soil, when the excess pore pressure increases. Thus, it is particularly important to understand the process of the wave-pipeline-soil interactions for the design of submarine pipelines [3]. The mechanisms of the wave-induced soil liquefaction can be classified into two categories, residual and oscillatory, in accordance with the way how the excess pore pressure is generated [4]. The residual liquefaction mechanism is resulted from the build-up of pore pressure induced by volumetric contraction under cyclic loading [5]. Momentary liquefaction usually appears in the seabed under wave troughs where the pore pressure is accompanied with some damping and phase lag [6]. This study focuses on the second mechanism.

Numerous investigations for the wave-induced soil response around submarine pipelines by adopting conventional numerical methods have been carried out since the 1980s, such as finite-element method (FEM), finite difference method (FDM) and boundary element method (BEM). Among these, Cheng and Liu [7] introduced the Boundary Integral Equation to obtain the distribution of pore-water pressure around a pipe fully buried in a sediment-filled trench. In their model, the inertia terms were considered, i.e., it is a $u-p$ approximation. Thomas [8,9] established a one-dimensional finite-element model to analyze the wave-induced soil response in saturated and unsaturated soil. Later, the model was extended to investigate the influences of variable permeability and shear modulus of soil and non-linearity of wave loading [10,11]. Furthermore, by using finite-element formulation, Madga [12,13,14] further estimated the wave-induced pore pressure and uplift force acting on submarine pipelines. Based on the assumption of no slipping at the interface between pipeline and soil, Jeng and Lin [15] developed a finite-element model to examine the wave-pipeline-seabed interactions in an in-homogeneous seabed. Then, the inertial forces and soil-pipe contact effects were involved in Luan et al. [16]’s model. With the commercial software (ABQUS), the effects of combined non-linear wave and current loading were considered by Wen et al. [17]. Recently, by considering pre-consolidation due to self-weight of the pipeline, Zhao et al. [18] investigated the build-up pore pressures. Recently, Zhao and Jeng [19] extended the integrated numerical model to evaluate the influence of backfilled depth of trench layer. Later, Duan et al. [20] proposed a two-dimensional model to investigate the wave and current-induced soil response around a partially buried pipeline in a trench layer.

The aforementioned investigations have employed the conventional approaches with meshes, for example, the principle of FEM is to divide the computational domain into small elements, and these elements do not overlap each other. A field function was established within each element by adopting simple interpolation functions. If the element is severely distorted, the shape function of this element would be of poor quality, which may lead to an unacceptable numerical result. Unlike the finite-element technique, the interpolation functions are established directly on nodes instead of elements, when meshless methods are applied. This could avoid the drawback of the conventional mesh-based techniques such as FEM and FDM. In recent years, meshless methods have attracted increasing attention from numerical modelers due to the faster formulation process with less data storage and no extensive mesh.

The most commonly used meshless methods are method of fundamental solution (MFS), method of particular solutions (MPS), smooth particle hydrodynamics (SPH), global radial basis functions (RBF) collocation method (GRBFCM) and local RBF collocation method (LRBFCM). One of the earliest meshless methods is SPH [21,22]. Randles and Libersky [23] was the first group to apply SPH in solid mechanics, and later this method was improved and was adopted in more fields. The main drawbacks of this method are inaccurate results near boundaries and tension instability, which was first studied in 1995. Over the ensuring decades, more meshless methods have been proposed. Karim et al. [24] presented a two-dimensional model using the element-free Galerkin method to investigate transient response of saturated porous elastic soil under cyclic loading system. A radial point interpolation meshless method (radial PIM) was developed by Wang et al. [25] to avoid the occurrence of singularity associated with only polynomial basis. Later, the radial PIM was applied to solve Biot’s consolidation problem [26] and wave-induced seabed response [27]. Existing meshless models have investigated water-soil interaction without any structure. Thus, the present study firstly establishes a meshless model to investigate the wave-seabed-pipeline interactions under various wave loading.

Among meshfree methods, the RBF are commonly used. The global RBF collocation method was first proposed for multivariate data interpolation and partial differential equations [28,29]. This method was further modified [30,31,32]. The GRBFCM can be used to deal with arbitrary and complex domains, but it usually leads to an ill-conditioned system matrix when high resolutions are required. Later, a localization procedure (LRBFCM) was proposed to overcome prescribed difficulty, to transform the dense system matrices into sparse ones [33]. Based on the multi-quadric RBF [34], Lee et al. [33] first proposed the local RBF collocation method (LRBFCM), and then this method was applied to various problems [35,36,37].

In this study, LRBFCM is employed to investigate the wave-induced oscillatory liquefaction around a pipeline in a trench layer. The proposed seabed model is validated with the analytical solution [38], experimental data [39,40] and numerical models [7,41]. Then, a parametric study is conducted to evaluate the influence of pipeline configuration, wave characteristics, and seabed properties for the wave-induced pore pressure around a partially buried pipeline in a trench.

2. Theoretical Models

In this study, an impermeable pipeline is considered with a radius of R which is partially buried in a trench with a finite thickness (h), as shown in Figure 1. The propagation direction of waves is along the positive x-direction.

The present model consists of two sub-models: wave and seabed sub-models. By using open-source CFD toolbox OpenFoam (Open Field Operation and Manipulation), the wave model proposed by Higuera et al. [42] is adopted for generating various wave loadings. A new seabed model is established based on LRBFCM.

By means of Volume-Averaged Reynolds-Averaged Navier-Stokes (VARANS) equations, the wave model was developed for coastal engineering applications by dealing with three-dimensional two-phase flow which was based on a solver in OpenFOAM, IHFOAM [42]. More detailed information about IHFOAM and its applications in coastal engineering, readers can refer to the publications of Higuera and his co-workers [42,43,44,45].

2.1. Boundary Value Problem for the Seabed Model

By neglecting the inertial effect, a quasi-static seabed model is established for the wave-induced seabed response based on the assumption of homogeneous seabed and compressible pore fluid. The effects of inertial terms on the wave-induced soil response has been reported in Jeng et al. [46]. For a two-dimensional problem, the governing equation for compressible homogeneous soil and compressible pore fluid can be represented as [47]:

$$K{\nabla }^{2}p-{\gamma }_w{n}^{\prime }\beta \frac{\partial p}{\partial t}-{\gamma }_w\frac{\partial ϵ}{\partial t}=0,$$

(1)

where p is pore-water pressure, ${\gamma }_w$ is the unit weight of water, $n^{\prime }$ is soil porosity, t is the time. In (1), the volume strain ($ϵ$) and compressibility of pore fluid ($\beta$) are defined as

$$ϵ=\frac{\partial u}{\partial x}+\frac{\partial w}{\partial z},\mathrm{and}\beta =\frac{1}{K_w}+\frac{1-S_r}{P_{wo}},$$

(2)

where u and w are the soil displacements in the x- and z-direction, respectively; $K_w$ is the true modulus of elasticity ($K_w=2\times {10}^9\mathrm{N}/{\mathrm{m}}^2$, [6]), $S_r$ is the degree of saturation and $P_{wo}$ is related to the absolute water pressure.

Based on Newton’s second law, the force balance for the porous flow in $x-$ and $z-$ directions can be expressed respectively as

$$\frac{\partial {\sigma }_x^{\prime }}{\partial x}+\frac{\partial {\tau }_{xz}}{\partial z}=\frac{\partial p}{\partial x},$$

(3)

$$\frac{\partial {\tau }_{xz}}{\partial x}+\frac{\partial {\sigma }_z^{\prime }}{\partial z}=\frac{\partial p}{\partial z},$$

(4)

where ${\sigma }_x^{\prime }$ and ${\sigma }_z^{\prime }$ are effective normal stresses; ${\tau }_{xz}$ denotes shear stress component. In this study, tension is determined as positive.

Based on Hook’s law, the effective normal stresses and shear stress can be expressed in term of soil displacements, i.e.,

$${\sigma }_x^{\prime }=2G\left[\frac{\partial u}{\partial x}+\frac{\mu }{1-2\mu }ϵ\right],$$

(5)

$${\sigma }_z^{\prime }=2G\left[\frac{\partial w}{\partial z}+\frac{\mu }{1-2\mu }ϵ\right],$$

(6)

$${\tau }_{xz}=G\left[\frac{\partial u}{\partial z}+\frac{\partial w}{\partial x}\right],$$

(7)

where the shear modulus G is defined with Young’s modulus (E) and the Poisson’s ratio ($\mu$) in the form of $E/2(1+\mu )$.

Substituting (5)~(7) into (3) and (4), the force equilibrium can be represented as

$$G{\nabla }^{2}u+\frac{G}{1-2\mu }\frac{\partial ϵ}{\partial x}-\frac{\partial p}{\partial x}=0,$$

(8)

$$G{\nabla }^{2}w+\frac{G}{1-2\mu }\frac{\partial ϵ}{\partial z}-\frac{\partial p}{\partial z}=0.$$

(9)

To solve the pore pressures and soil displacements in (1), (8) and (9), the following boundary conditions are required.

• At seabed surface ($z=0$) and trench surface, the vertical effective stress and shear stress vanish, and the pore pressure is equal to dynamic wave pressure.

  $${\sigma }_z^{\prime }={\tau }_{xz}=0,\mathrm{and}p=P_b,$$

  (10)

  where $P_b$ is the dynamic wave pressure at the seabed surface, which is obtained from the wave model (IHFOAM).
• At the impermeable seabed bottom ($z=-h$), zero displacements and no vertical flow are specified, i.e.,

  $$u=w=0,\mathrm{and}\frac{\partial p}{\partial z}=0,$$

  (11)
• The pipeline surface is assumed to be impermeable wall. Thus, there is no flow through the pipeline surface, i.e.,

  $$\frac{\partial p}{\partial n}=0,$$

  (12)

  where n denotes normal vector of the pipe surface.

2.2. Meshfree Model for the Seabed Domain

In this study, a rigid pipeline is considered to be partially buried in a trench. The computational domain is discretized into N nodes non-uniformly. Therefore, a linear equation of the following form is required to be established:

$${\left[A\right]}_{N\times N}{\left[\mathsf{\Phi }\right]}_{N\times 1}={\left[B\right]}_{N\times 1},$$

(13)

where ${\left[\mathsf{\Phi }\right]}_{N\times 1}$ is the sought solution, ${\left[B\right]}_{N\times 1}$ is a column vector, and ${\left[A\right]}_{N\times N}$ is a sparse system matrix. Similar structures of ${\left[A\right]}_{N\times N}$ can be found in the FDM and the finite-element method.

For constructing a linear equation for each node ${\mathbf{y}}_n$ in the computational domain, $\mathsf{\Phi }$ in (13) is assumed as $\mathsf{\Phi }\left(\mathbf{x}\right)$ by RBFs:

$$\mathsf{\Phi }\left(\mathbf{x}\right)\approx \sum _{m=1}^{\overline{K}}{\alpha }_m\chi \left(r_m\right),$$

(14)

where $\mathsf{\Phi }$ denotes either p or $u_i$ in the governing equations, ${\alpha }_m$ refers to the corresponding undetermined coefficient and $r_m=\parallel \mathbf{x}-{\mathbf{x}}_m\parallel$ is the Euclidean distance from $\mathbf{x}$ to ${\mathbf{x}}_m$. The group of ${\mathbf{x}}_m$ denote the locations of the $\overline{K}$ nearest neighbor nodes surrounding the prescribed center ${\mathbf{x}}_1$. In this study, the kd-tree algorithm is applied to search the $\overline{K}$ nearest neighbor nodes efficiently [48]. Furthermore, the multi-quadric RBF is expressed as

$$\chi \left(r_m\right)=\sqrt{{r_m}^2+c^2},$$

(15)

with the shape parameter (c) [34].

A localization process [33,35,49] is presented here for the sake of preventing unnecessary ill-conditioned system matrix. Firstly, the expression of $r_m$ is substituted into (14) as

$$\mathsf{\Phi }\left({\mathbf{x}}_n\right)=\sum _{m=1}^{\overline{K}}{\alpha }_m\chi (\parallel {\mathbf{x}}_m-{\mathbf{x}}_n\parallel ),$$

(16)

or in matrix-vector form as

$${\left[\mathsf{\Phi }\right]}_{\overline{K}\times 1}={\left[\chi \right]}_{\overline{K}\times \overline{K}}{\left[\alpha \right]}_{\overline{K}\times 1},$$

(17)

where

$${\left[\mathsf{\Phi }\right]}_{\overline{K}\times 1}=\left[\begin{array}{c}\mathsf{\Phi }\left({\mathbf{x}}_1\right)\\ \mathsf{\Phi }\left({\mathbf{x}}_2\right)\\ ⋮\\ \mathsf{\Phi }\left({\mathbf{x}}_{\overline{K}}\right)\end{array}\right],$$

(18)

$${\left[\chi \right]}_{\overline{K}\times \overline{K}}=\left[\begin{array}{cccc}\chi (\parallel {\mathbf{x}}_1-{\mathbf{x}}_1\parallel )& \chi (\parallel {\mathbf{x}}_1-{\mathbf{x}}_2\parallel )& \cdots & \chi (\parallel {\mathbf{x}}_1-{\mathbf{x}}_{\overline{K}}\parallel )\\ \chi (\parallel {\mathbf{x}}_2-{\mathbf{x}}_1\parallel )& \chi (\parallel {\mathbf{x}}_2-{\mathbf{x}}_2\parallel )& \cdots & \chi (\parallel {\mathbf{x}}_2-{\mathbf{x}}_{\overline{K}}\parallel )\\ ⋮& ⋮& \ddots & ⋮\\ \chi (\parallel {\mathbf{x}}_{\overline{K}}-{\mathbf{x}}_1\parallel )& \chi (\parallel {\mathbf{x}}_{\overline{K}}-{\mathbf{x}}_2\parallel )& \cdots & \chi (\parallel {\mathbf{x}}_{\overline{K}}-{\mathbf{x}}_{\overline{K}}\parallel )\end{array}\right],$$

(19)

and

$${\left[\alpha \right]}_{\overline{K}\times 1}=\left[\begin{array}{c}{\alpha }_1\\ {\alpha }_2\\ ⋮\\ {\alpha }_{\overline{K}}\end{array}\right].$$

(20)

Then, (17) can be inverted as

$${\left[\alpha \right]}_{\overline{K}\times 1}={{\left[\chi \right]}_{\overline{K}\times \overline{K}}}^{-1}{\left[\mathsf{\Phi }\right]}_{\overline{K}\times 1}.$$

(21)

Now, $L\mathsf{\Phi }\left(\mathbf{x}\right)$ is considered to replace $\mathsf{\Phi }\left(\mathbf{x}\right)$ defined in (14), where L is a linear differential operator related to both the governing equation and the boundary conditions. The collocation of $L\mathsf{\Phi }\left(\mathbf{x}\right)$ on ${\mathbf{x}}_1={\mathbf{y}}_n$ gives

$$L\mathsf{\Phi }\left({\mathbf{y}}_n\right)=\sum _{m=1}^{\overline{K}}{\alpha }_{m}L\chi \left(r_m\right){\mid }_{\mathbf{x}={\mathbf{x}}_1},$$

(22)

or in matrix-vector form as

$$L\mathsf{\Phi }\left({\mathbf{y}}_n\right)={\left[L\chi \right]}_{1\times \overline{K}}{\left[\alpha \right]}_{\overline{K}\times 1}.$$

(23)

In (23), the existence of $L\chi \left(r_m\right)$ is as a result of the influence of operator L on the RBF $\chi \left(r_m\right)$. Thus, (17) and (23) can be combined as

$$L\mathsf{\Phi }\left({\mathbf{y}}_n\right)={\left[C\right]}_{1\times \overline{K}}{\left[\mathsf{\Phi }\right]}_{\overline{K}\times 1},$$

(24)

with

$${\left[C\right]}_{1\times \overline{K}}={\left[L\chi \right]}_{1\times \overline{K}}{{\left[\chi \right]}_{\overline{K}\times \overline{K}}}^{-1},$$

(25)

and

$${\left[L\chi \right]}_{1\times \overline{K}}=\left[\begin{array}{cccc}L\chi \left(r_1\right){\mid }_{\mathbf{x}={\mathbf{x}}_1}& L\chi \left(r_2\right){\mid }_{\mathbf{x}={\mathbf{x}}_2}& \cdots & L\chi \left(r_{\overline{K}}\right){\mid }_{\mathbf{x}={\mathbf{x}}_{\overline{K}}}\end{array}\right].$$

(26)

From (24)–(26), it can be found that the row vector ${\left[C\right]}_{1\times \overline{K}}$ can be obtained if all the values of L, $\chi$ and ${\mathbf{x}}_j$ are known. These equations can be assembled into the system matrix, and finally the resultant sparse system is solved by using the direct solver of SuperLU in this study, which finished the procedure of LRBFCM.

Please note that the radial PIM was adopted to solve Biot’s consolidation problem [26] and wave-induced soil response [27], while the present model uses LRBFCM. The mixed bases of polynomial and radial bases are needed in the radial PIM for the accuracy of polynomials [25]. Compared with the radial PIM, the choice of basis functions in LRBFCM is easier. Furthermore, no submarine structure was included in their model [27]. Thus, the present seabed model is the first model by applying LRBFCM to investigate the wave-seabed interactions around a structure such as pipelines.

2.3. Effects of Lateral Boundary Conditions

This section presents two ways to handle the lateral boundary conditions: periodic and fixed. Generally speaking, the horizontal and vertical displacements and the pore pressure do not vanish at lateral boundaries. To deal with the problem of boundaries for wave-seabed interactions, Jeng et al. [50] applied the principle of repeatability [51]. However, the condition of employing periodic boundary condition is that the length of seabed must be an integer number of wavelength. Moreover, periodic boundary condition is not applicable for seabed with structures. Thus, Ye and Jeng [52] suggested another method by which employing a large computational domain and meanwhile fixing both the lateral boundaries in the horizontal direction, namely considering the boundary as impermeable. This method is under the assumption that the fixed lateral boundary only influences the region nearby. Similar with Ye and Jeng [52], in this section, both periodic and fixed lateral boundary conditions by LRBFCM are examined.

Theoretically, the larger the computational domain, the smaller effects of lateral boundaries. However, a large computational domain requires more computational resources. Thus, the length of computational domain is assumed as 3 times of the periodic wavelength in this study. The input data used is listed in

Table 1. Input data for numerical examples demonstrating effects of lateral boundary conditions.

Wave Characteristics
Wave period $\mathit{T}$	8.0 s
Water depth $\mathit{d}$	20 m
Wave length $\mathit{L}$	88.88 m
Soil Characteristics
Thickness of seabed $\mathit{h}$	30 m
Poisson’s ratio $\mathit{\mu }$	0.33333
Soil porosity $\mathit{n}$	0.3
Soil permeability $\mathit{K}$	${10}^{-2}$ m/s
Degree of saturation ${\mathit{S}}_{\mathit{r}}$	0.98
Shear modulus $\mathit{G}$	${10}^7$ N/m${}^2$

. The maximum pore pressure and effective stresses of soil of these three sections are depicted in Figure 2.

As depicted in Figure 2, sections $A-A^{\prime }$ ($x=50$ m), $B-B^{\prime }$ ($x=125$ m) and $C-C^{\prime }$ ($x=200$ m) are in the range of the first wave length, the second wavelength and the third wavelength, separately. Solid lines and dashed lines represent the soil response of the case with periodic and fixed boundaries, respectively. As shown in Figure 2a,c, the effect of fixed boundary condition is minor for vertical effective stress and pore pressure in sections $A-A^{\prime }$ and $C-C^{\prime }$, but considerable difference can be observed from the horizontal effective stress, which leads to the conclusion that fixed lateral boundaries affect the soil response of the soil region near lateral boundaries significantly. Furthermore, the seabed response at section $B-B^{\prime }$ under periodic boundary conditions is in a complete agreement with that under fixed boundary conditions, which means that the influence of fixed lateral boundaries vanish in the section far away from the boundary. Thus, fixed lateral boundary condition is employed in this study for investigating the wave-induced soil response.

2.4. Convergent Tests

The present seabed model is a new model established by employing LRBFCM for the wave-induced soil response in the vicinity of submarine pipelines, it is necessary to check its convergence. Figure 3 presents three tests for model convergence which is with respect to node number of the whole computational domain, a model parameter c (where c is equal to “$30\times$ the maximum distance between each two nodes in the local region”) and the value of $\overline{K}$ (where $\overline{K}$ indicates the node number of the local region). The pipeline is considered to be impermeable, and input data are as follows: $d=0.533$ m, $L=1.25$ m, $h=0.826$ m, $\mu =0.33$, $n=0.42$, $K=0.0011$ m/s, $S=0.997$, $R=0.084$ m, $b=0.167$ m, $L_x=4.57$ m.

In principle, instability of the trend of soil response around pipeline occurs at the beginning stage of node number increase. The results should remain unchanged after the node number is increased to a certain extent. With fixed node numbers, the numerical results should not be changed when the value of c is in a reasonable range, which can prove that the model is convergent and reliable. As presented in Figure 3a, the wave-induced pore pressure keep changed when the node number varies from 16,558 to 45,000 approximately, but the values maintain a steady state in the process of the node number increases from 45,375 to 53,351 and even 65,231, which verifies the stability of the model. c is one of coefficients of the present model. In Figure 3b, the node number is determined as 45,375, then the trend of pore pressure can be observed through changing the value of c. There is almost no change for the pore pressure in the vicinity of the pipeline when c is equal to 0.3, 0.548 and 0.8, respectively, which can be evidence of the model convergence. $\overline{K}$ refers to the number of the nearest neighbor nodes of unknown node x. Usually, the value of $\overline{K}$ can be regarded as 5, 9, and 13. From Figure 3c, it can be found that 9 or 13 is applicable for the present model, and the result looks more smooth when $\overline{K}$ is equal to 9. Thus, the convergence of present model is verified from these three cases. Combined Figure 3a–c, it can be concluded that the numerical result of this case scale is satisfactory when the node number, and the value of c and $\overline{K}$ are determined as 45375, 0.548 and 9, respectively.

3. Model Validation

3.1. Comparison with the Analytical Solution for Wave-Seabed Interactions

For a homogeneous seabed, the previous analytical solution for the wave-induced soil response [38] will be compared with present numerical results. This comparison is to confirm the capacity of the present model.

In this comparison, the following input data are used: wave period $T=15.0$ sec, water depth $d=70$ m, wavelength $L=311.59$ m, thickness of seabed $h=25$ m, Poisson’s ratio $\mu =0.333$, soil porosity $n=0.3$, soil permeability $K={10}^{-4}$ m/s for fine sand and ${10}^{-2}$ m/s for coarse sand, degree of saturation $S=0.932$ for unsaturated soil, shear modulus $G={10}^7{\mathrm{N}/\mathrm{m}}^2$. The numerical results of the comparison are presented in Figure 4. In the figure, the present results are presented by lines and the analytical solution [38] is denoted as circles. The vertical distributions of the maximum amplitude of the wave-induced pore pressure ($\left|p\right|/p_0$) and effective stresses ($|{\sigma }_x^{\prime }|/p_0$, $|{\sigma }_z^{\prime }|/p_0$), and shear stress ($|{\tau }_{xz}/p_0|$) versus $z/h$ are presented. In the figure, $p_0$ is the amplitude of linear wave pressure at the seabed surface, which is defined as $p_0={\gamma }_{w}H//2coshkd$. It is found that the present results are in complete accordance with the analytical solution of Hsu and Jeng [38]. The difference between the analytical solution and the present model is less than ${10}^{-3}$ for both fine and coarse sands.

3.2. Comparison with Experimental Data and FEM Results for Wave-Pipeline-Seabed Interaction

The second validation is the comparison between present numerical results and experimental results [39] with respect to the linear wave-induced soil response around a fully buried pipeline. Considering an impermeable pipeline with a radius of R is fully buried within a porous elastic seabed with a finite thickness (h). The propagation direction of waves is regarded along the positive x-direction. The input data employed in this validation is same as Section 2.4. Turcotte et al. [39] reported seven experiments with different wave period and wave height, and only three typical comparisons are presented here which are that with the longest wave length ($L=4.91$ m, $T=2.3$ s, $H=0.0302$ m), the medium wave length ($L=3.536$ m, $T=1.75$ s, $H=0.143$ m) and the shortest wave length ($L=1.25$ m, $T=0.9$ s, $H=0.0524$ m). The water depth was a constant of 0.533 m for all tests. In the comparison, the results from previous boundary element model [7] and the COMSOL finite-element model are also included. The COMSOL model was based on the one proposed by Jeng and Zhao [41] and applied to the case with a buried pipeline.

From Figure 5a, slight difference can be found from the comparison between the numerical results and experimental data. For the intermediate wave period illustrated in Figure 5b, the present model is the closest to the experimental data than other two numerical models. From the third circumstance shown in Figure 5, it can be seen that the error of the finite-element model with the experimental data is relatively higher than the error of present results with the experimental data. Furthermore, the length of the computational domain was fixed as 4.57 m. Hence, the fixed lateral boundary condition is not applicable for the situation of $L=4.91$ m. However, Figure 5 still presents that case for a complete comparison with the numerical solution of Cheng and Liu [7]. It can be observed that the amplitude of wave-induced pore pressure increases with the increase of the wavelength.

3.3. Comparison with Experimental Data for Wave-Induced Soil Response Around a Pipeline Buried in a Trench

The third validation is to compare the model with experimental data [40] for the case of a homogeneous seabed, in which the pipeline is in a trench. Stoke II wave loading simulated with OpenFoam is employed in this case. Sun et al. [40] conducted a series of laboratory experiments to examine the wave-induced pore-water pressure along the surface of a pipeline partially backfilled in a trench, and for brevity of presentation, only two typical comparisons are shown in this section in which the wave and seabed conditions are listed in

Table 2. Input data for the comparison with experimental data [40].

Wave Characteristics
Water depth $\mathit{d}$	0.4 m
Soil Characteristics
Thickness of seabed $\mathit{h}$	0.58 m
Poisson’s ratio $\mathit{\mu }$	0.32
Soil porosity $\mathit{n}$	0.396
Soil permeability $\mathit{K}$	$3.56\times {10}^{-5}$ m/s
Degree of saturation ${\mathit{S}}_{\mathit{r}}$	0.998
Geometry of the Pipe
Pipe radius $\mathit{R}$	0.05 m

and

Table 3. Wave and seabed conditions for the comparison with experimental data [40].

Case No.	Wave Condition		Seabed Condition
	Wave Height H (m)	Wave Period T (s)	Trench Depth (m)	Backfill Depth (m)
45	0.12	1.6	0.15	0.1
46	0.12	1.6	0.15	0.125

. The corresponding comparisons of wave-induced pore pressure around the pipeline are illustrated in Figure 6. Red lines denote present numerical results and circle denote experimental data. It can be observed that the present numerical model agrees well with the laboratory experiments.

4. Results and Discussion

The aim of this study is to investigate the wave-pipeline-seabed interactions around a trenched pipeline by employing the proposed time-dependent meshless seabed model. In this section, the influence of soil properties, wave characteristics, and pipe configuration on the wave-induced oscillatory liquefaction are examined. Zen and Yamazaki [53] introduced and verified the concept of “oscillatory” excess pore pressure by conducting a series of experiments, and their criterion to determine the soil oscillatory liquefaction is used in this study, which can be expressed as

$${\sigma }_0^{\prime }(z,0)\le u_e(z,t)=-[p(0,t)-p(z,t)]$$

(27)

where ${\sigma }_0^{\prime }(z,0)$ is the initial effective stress, and $u_e(z,t)=-[p(0,t)-p(z,t)]$ means the excess pore pressure. $p(0,t)$ and $p(z,t)$ denote the wave pressure at seabed surface and wave-induced pore pressure, respectively.

As shown in Figure 1, the pipeline is buried in the partially backfilled trench, and the lateral boundaries are considered as impermeable in this case. The direction of wave propagates along the positive x-direction. The Stokes II wave loading is simulated by OpenFOAM. Wave and soil parameters are listed in

Table 4. Input data for the parametric study.

Wave Characteristics
Wave period $\mathit{T}$	10 s
Water depth $\mathit{d}$	8 m
Wave height $\mathit{H}$	3 m
Soil Characteristics
Thickness of seabed $\mathit{h}$	30 m
Poisson’s ratio $\mathit{\mu }$	0.35
Soil porosity $\mathit{n}$	0.425
Soil permeability $\mathit{K}$	$1\times {10}^{-3}$ m/s
Degree of saturation ${\mathit{S}}_{\mathit{r}}$	0.98
Geometry of the Pipe
Pipe radius $\mathit{R}$	0.4 m
Backfilled depth ${\mathit{H}}_{\mathit{b}}$	0.5 m

. When the effect of one parameter for the wave-induced pore pressure is examined, values of other parameters are kept fixed.

4.1. Effects of Soil Characteristics

Soil characteristics are significant factors to affect the wave-induced oscillatory soil response in the vicinity of a partially buried pipeline. In this section, two parameters are examined in detail. They are the degree of saturation $S_r$ and soil permeability K. Figure 7 presents the distribution of liquefaction depth around pipeline in the trench under various soil conditions, and Figure 8 shows the distributions of the maximum excess pore pressures under wave trough in the vertical section through the center of the pipeline. For Figure 7a and Figure 8a, permeability is $4.5\times {10}^{-3}$ m/s. For Figure 7b and Figure 8b, degree of saturation is 0.98. Soil properties of backfills are chosen as same as bottom soil.

Three typical values of degree of saturation are considered in this section, they are: 95%, 97% and 99%, respectively. Figure 7a demonstrates that degree of saturation ($S_r$) significantly affects the liquefaction depth in the trench. The depth is deeper with decreasing degree of saturation. Furthermore, from Figure 8a, it can be found that the soil on the bottom of the pipeline is much easier to be liquefied when degree of saturation is relatively small.

To investigate the influence of soil permeability on the wave-induced soil response, Figure 7b illustrates the distribution of the liquefaction depth in the trench with variable value of permeability, $4\times {10}^{-3}$ m/s, $1\times {10}^{-3}$ m/s and $7\times {10}^{-4}$ m/s, respectively. Figure 8b shows the distribution of the maximum excess pore pressure under wave trough along the vertical section through the center of pipeline with variable permeability. It can be concluded that the liquefaction depth become large with decreasing permeability, and soil around the impermeable pipeline is much easier to be liquefied when permeability is relatively small.

4.2. Effects of Wave Characteristics

In addition to soil characteristics, wave parameters have been found to significantly influence the wave-pipeline-seabed interactions. The influences of two wave parameters, wave height (H) and period (T), are examined in this section. The wave height H can directly affect the magnitude of wave loading exerting on the seabed surface. Wave period T affects the wave-induced oscillatory excess pore pressure by affecting wavelength. Figure 9 illustrates the distribution of oscillatory liquefaction depth in the trench under various wave height (H) and period (T), and the wave period is 10 s for Figure 9a, the wave height (H) is 3 m for Figure 9b. Figure 10 shows the distribution of the maximum excess pore pressure $|u_e|/{\sigma }_0^{\prime }$ along the vertical section under the center of the pipeline.

Three values of wave height are examined in this section: 2 m, 2.5 m and 3.1 m. From Figure 9a, it can be found that the liquefaction depth increases with increasing wave height. Meanwhile, from Figure 10a, it can be seen that soil under the pipeline is much easier to be liquefied when wave height is relatively large.

Figure 9b and Figure 10b present the effect of wave period on the wave-induced oscillatory excess pore pressure around the trenched pipeline. Liquefaction depth from the seabed surface increases when wave period increases, and the soil on the bottom of pipeline is much easier to be liquefied when wave period is large.

4.3. Effects of Backfill

In this section, the effect of backfill depth ($H_b$) on the wave-induced oscillatory liquefaction around a pipeline buried in a trench under the Stokes II wave loading is investigated. Three variable backfill depth are examined: 0.3 m, 0.5 m and 0.9 m. Figure 11 depicts the distribution of oscillatory liquefaction depth around the partially buried pipeline in a trench for various backfill depths, and Figure 12 illustrates the distribution of oscillatory excess pore pressure of the vertical section on the bottom of the pipeline under wave trough for the same four backfill depths. Figure 11 demonstrates that the liquefaction depth is greater with decreasing backfill depth. Similarly, from Figure 12, it can be seen that the maximum excess pore pressure increases as the backfill depth decreases.

5. Conclusions

This study proposes a two-dimensional seabed model by LRBFCM to investigate the wave-induced oscillatory liquefaction around a partially buried pipeline in a trench under non-linear wave loading. The model is validated by comparison with analytical solution, experimental data, and previous numerical results. The effects of wave characteristics, soil properties, and backfill depth in the trench are examined. The following conclusions can be drawn:

(1)

Unlike previous investigations using conventional numerical methods, this study established a meshless seabed model by employing LRBFCM and applied it to examine the wave-induced soil response. The validation with the analytical solution [38] and experimental data [39,40] shows that present model is satisfactory.

(2)

The wave-induced oscillatory excess pore pressure is relatively susceptible to the adjustment of degree of saturation ($S_r$) and permeability (K) of soil. Low values of $S_r$ and K lead to great magnitude of wave-induced excess pore pressure around the pipeline.

(3)

Oscillatory liquefaction depth is influenced significantly by wave characteristics, such as wave height (H) and wave period (T). Figure 9 shows that the liquefaction depth is deeper with increasing wave height (H) and wave period (T).

(4)

Pipe configuration is significantly important for the analysis of wave-pipeline-seabed interaction. In the process of increasing buried depth of pipe, the magnitude of oscillatory excess pore pressure at the bottom of the trenched pipeline decreases, which means that relatively large value of backfill depth can reduce the risk of liquefaction.

When the conventional methods with meshes are applied to analyze the computational domain with irregular boundaries, the elements or meshes may be distorted. Interpolation and re-meshing can be used to solve this problem. However, it requires more intensive work for complicated engineering problems. The meshless model presented in this study is designed for avoiding the poor mesh quality existing in conventional models. However, it needs to be further developed for different engineering problems if a huge number of nodes are required.

This study focuses on the wave-induced soil response under wave loading in two-dimensional. However, in real ocean environments, waves may approach the pipeline from any direction. Therefore, the effect of wave oblique on the soil response in the vicinity of pipelines will be examined in the future.