Response of a Porous Seabed around an Immersed Tunnel under Wave Loading: Meshfree Model

Abstract

Seabed instability surrounding an immersed tunnel is a vital engineering issue regarding the design and maintenance for submarine tunnel projects. In this study, a numerical model based on the local radial basis function collocation method (LRBFCM) is developed to evaluate the seabed behaviour in a marine environment, in which the seabed is treated as the porous medium and governed by Biot’s “$u-p$” approximation. As for the flow field above the seabed, the VARANS equations are used to simulate the fluid motion and properties. The present model is validated with analytical solutions and experimental data which show a good capacity of the integrated model. Both wave and current loading are considered in this study. Parametric studies are carried out to investigate the effects of wave characteristics and soil properties. Based on the numerical results, the maximum liquefaction depth around the immersed tunnel could be deeper under the wave loading with long wave period (T) and large wave height (H). Moreover, a seabed with lower permeability ($K_s$) and degree of saturation ($S_r$) is more likely to be liquefied.

1. Introduction

In recent years, to meet the continual improvement requirements in coastal transportation, the immersed tunnel has become one of the choices to fundamentally transform the transport in the region of oceans and rivers, replacing the conventional methods such as the ferry. The immersed tunnel has a history of about 100 years and it shows a good performance in reliability and applicability under complex natural dynamic loading; for example, the longest immersed tunnel in the world, the Hongkong–Zhuhai–Macao Bridge immersed tunnel. As an alternative to a bridge, the immersed tunnel has advantages in less environmental effects and no obstruction of navigation channels. Compared to bridges, the immersed tunnel is commonly constructed in a soft and loose seabed. Thus, the stability of the seabed soil surrounding the immersed tunnel in complex marine environments becomes one of the main concerns implicated in tunnel design and maintenance.

It has been recognized that the pore water pressures and stresses in seabeds are affected by the water pressures generated by the natural dynamic loading. If the pore water pressure reaches the limit value, the liquefaction could occur with the effective stress in seabed vanishing. To avoid seabed instability around the immersed tunnel, the study of seabed dynamic behaviour is necessary under the real hydrodynamic loading. Two mechanisms of wave-induced liquefaction has been figured out based on a mass of laboratory tests and field exploration [1,2,3], which are transient liquefaction and residual liquefaction. The transient liquefaction is motivated by the oscillatory excess pore water pressures under wave pressure vibration which usually happens with amplitude reduction and phase lag of pore pressure in seabed soil [4]. While the residual liquefaction is on the consequence of the excess pore water pressure build-up under cyclic wave loading [5]. Later, Jeng and Seymour [6] proposed a simplified approximation to predict the liquefaction process in a large seabed, and concluded that the residual mechanism is more essential under large waves while the transient mechanism is dominate in a seabed under small wave loading.

In the past 40 years, various analytical formulas have been developed and verified in regard to the seabed dynamic response [4,7,8]. However, the seabed dynamic response around the structure is difficult to be described by analytical methods. Thus, several numerical model was developed to simulate the soil behaviour around the offshore structures, such as breakwaters [9], and buried pipes [10]. However, the research of a wave-induced response around the undersea immersed tunnel is quite limited. For instance, a real case simulation of the Busan–Geoje fixed link in South Korea was conducted by Kasper et al. [11] under large waves (wave height up to 9.2 m) generated by typhoons. Nevertheless, this study did not consider the impact of the seabed liquefaction around the tunnel. Recently, [12,13] simulated the seabed transient and residual response around the immersed tunnel under wave loading based on Biot’s consolidation equations neglecting the inertial terms for soil skeleton and fluid phase.

In natural ocean environments, current is another crucial component besides wave. For instance, a long-term mentoring data of the Lingding Bay, in which the Hong Kong–Zhuhai–Macau bridge tunnel located, shows the current co-exists with wave varying from bottom to surface as the consequence of the irregular semi-diurnal tide. The maximum velocity of the surface current reaches up to 1.5 m/s [14]. The interaction between current and wave is found to be able to affect the hydrodynamic properties directly and further impact the porous seabed dynamics. Ye and Jeng [15] investigated the transient response of a porous seabed under wave combined current loading firstly. Later, the current effects in the vicinity of a submarine pipeline were examined by Wen et al. [16] based on the commercial software ABAQUS. Lately, Liao et al. [17] simulated the residual seabed liquefaction under the flow field that wave and current generated simultaneously. To date, how current affect the liquefaction of seabed soil surrounding the immersed tunnel has not been examined yet to the author’s knowledge.

The aforementioned numerical models mainly developed adopting the conventional methods such as the finite element method and finite difference method. In recent years, a new numerical technique has come up which uses a set of nodes instead of the conservative meshes to approximate the solution. In consequence of discretization of the partial differential equations directly on nodes instead of meshes, the meshfree method is more qualified in dealing with mesh entanglement problems and constructing the approximations with arbitrary order of continuity than the traditional mash-based numerical methods. The start point of the meshfree method was using the moving least-square (MLS) method to establish shape functions for a set of scatter nodes by Nayroles et al. [18]. After that, Belytschko et al. [19] used the Galerkin weak form to improve the diffuse element method, which formed a new element-free Galerkin method. The element-free Galerkin method has been widely used in soil mechanics problems. In addition, an innovative interpolation scheme to overcome the disadvantage of the MLS was proposed by Liu and Gu [20], which is called the point interpolation method (PIM). The original PIM method picks up the polynomial basis as the basis function which performs well for one dimensional problems. However, it is hard to determine the ranks when this method extended to multi-dimensional range, which is the result of basis function selection. In order to figure this problem, Wang et al. [21] proposed a point interpolation method based on the radial basis functions (RBF). This method maps the multi-dimensional space into one-dimensional space though a radial function, which makes choosing basis function easier. This method has been widely used on the geomechanics problems. For example, Wang and Liu [22] simulated the Biot’s consolidation process by radial PIM method. The displacement and the pore water pressure were discretised by the same shape function in space, while the fully implicit integration scheme was adopted for time discretization to avoid the spurious ripple effect. In the literatures, two different type of the RBFs are commonly used, which are the multiquadrics (MQ) [23], and the Gaussian radial basis functions [24]. In this study, the MQ was used. These functions were first under intensive research in multivariate data interpolation and used to solve the partial differential equations by Kansa [25]. Thus, this method is usually called Kansa’s method, and is known as the global RBF collocation method (GRBFCM). This method has been applied to computational fluid dynamics problems such as solutions of the Navier–Stokes equation [26], natural convections in porous media [27], numerical wave tanks [28], and solid–liquid phase change problems [29]. The GRBFCM has the obvious advantages in dealing with arbitrary and complex domain, which applied in many field, such as surface fitting, turbulence analysis, neural networks, meteorology. However, ill-conditions can occur when the resolution is high. For the purpose of overcoming the difficulties of ill-conditions as mentioned above, a localization procedure to transform the dense matrices into sparse matrices was proposed. On the basis of the multiquadric RBF [30], Lee et al. [31] proposed the local RBF collocation method (LRBFCM) for the first time. This method has been applied in various fields, such as the solutions of diffusion problem [32], Darcy flow in porous media [33], water wave scattering [34], macro-segregation phenomena [35] and so on. Recently, Wang et al. [36] adopted the LRBFCM based on the Biot’s consolidation equation to simulate a wave-induced seabed response around a pipeline.

In this paper, an integrated numerical model is proposed to simulate the sandy seabed dynamic response and transient liquefaction in the vicinity of an immersed tunnel under natural complex loading. The flow model is developed based on IHFOAM [37], while the seabed model is established adopting the LRBFCM with Biot’s “$u-p$” approximation which considered the inertial term of soil skeleton. The verification of the integrated model is carried out by comparisons with the analytical solution [7] and published experimental data [38]. In this study, the effects of the current on the seabed behaviour around the immersed tunnel are examined. Parametric studies are conducted in regard of the different wave characteristics and soil properties.

2. Theoretical Models

In this study, an integrated 2D numerical model is developed to simulate the fluid-structure seabed interaction, which is composed of two sub-models: the flow model and seabed model. The flow model is responsible for simulating the wave motion including the wave pressure and current velocity, while the seabed model evaluates the effective stress, displacement, and pore pressure in the seabed under the dynamic loading. The two sub-models adopt the one-way coupling algorithm which is connected by the continuous water pressure.

2.1. Flow Model

In this paper, the flow model is based on one of the solvers in OpenFOAM^®, IHFOAM. Recently, some other open-source codes have been developed for modeling wave propagation and wave–structure interaction problems [39,40] as well. To simulate the coastal, offshore, and hydraulic engineering process, this model solves the three-dimensional Volume Averaged Reynolds Averaged Navier–Stokes (VARANS) equations with regard to two incompressible phases of air and water. The fluid model adopts the finite volume discretization and the volume of fluid (VOF) method [37]. Including continuity and momentum conservation equations; the VARANS equations as the governing mathematical expressions in this model can be expressed as:

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

(1)

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

(2)

where $\langle \rangle$ is Darcy’s volume averaging operator, while ${\langle \rangle }^f$ is the intrinsic averaging operator; ${\rho }_f$ represents the density which is computed by ${\rho }_f=\alpha {\rho }_{water}+(1-\alpha ){\rho }_{air}$; $\alpha$ is the indicator function defined in (5); $u_i$ is the velocity vector while n is the porosity; $p^{\ast }$ is the pseudo-dynamic pressures; $\mathbf{g}$, $g_i$ is the gravitational acceleration; ${\mu }_{eff}$ is the efficient dynamic viscosity, defined as ${\mu }_{eff}=\mu +{\rho }_f{\nu }_{turb}$, in which $\mu$ is the molecular dynamic viscosity and ${\nu }_{turb}$ is the turbulent kinetic viscosity, given by the chosen turbulence model. The $k-ϵ$ turbulence model is adopted in this study; the last term in (2) corresponds to the resistance of porous media, which is shown as:

$$\left[CT\right]=A\langle u_i\rangle +B\left|\langle u\rangle \right|\langle u_i\rangle +C\frac{\partial \langle u_i\rangle }{\partial t},$$

(3)

compared to factors A and B, factor C is less significant. Moreover, $0.34$ (kg/m${}^3$) is often adopted for the value of C by default [41].

In present fluid model, the values of A and B are derived by Engelund [42]’s formulae, which also employed in Burcharth and Andersen [43].

$$A=E_1\frac{{(1-n)}^3}{n^2}\frac{\mu }{d_{50}^2},\mathrm{and}B=E_2(1+\frac{7.5}{KC})\frac{1-n}{n^2}\frac{{\rho }_f}{d_{50}},$$

(4)

in which $d_{50}$ is the medium grain diameter of the materials; $KC$ is the Keulegan Carpenter number, which is defined as $KC=u_m{T}_0/\left(n{d}_{50}\right)$; $u_m$ is the maximum oscillating velocity while $T_0$ is the period of the oscillation. $E_1$ and $E_2$ are parameters that define the linear and non-linear friction terms. The default values of these parameters ($E_1=50$ and $E_2=1.2$) are used in this study [44].

A two-phased fluid mixture containing air and water is taken into consideration in each cell, which can be controlled by an indicator function $\alpha$. In the present model, $\alpha$ is defined as the quantity of water per unit of volume, which varies from 0 (air) to 1 (water):

$$\alpha =\left\{\begin{array}{cc}1,\hfill & \mathrm{water}\\ 0,\hfill & \mathrm{air}\\ 0<\alpha <1,\hfill & \mathrm{free} \mathrm{surface}\end{array}\right.$$

(5)

As an phase function, $\alpha$ can represent any variation of fluid properties considering the mixture properties, such as density and viscosity:

$$\mathsf{\Phi }=\alpha {\mathsf{\Phi }}_{water}+(1-\alpha ){\mathsf{\Phi }}_{air},$$

(6)

where ${\mathsf{\Phi }}_{water}$ and ${\mathsf{\Phi }}_{air}$ stand for the properties of water and air, separately, such as density of the fluid.

The fluid movement could be traced by solving the advection equation as [44]:

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

(7)

where $|u_c|=min\left[c_{\alpha }\left|u\right|,max\left(\right|u\left|\right)\right]$, in which the default value of $c_{\alpha }$ is 1, however, the user can specify a greater value to enhance the compression of the interface, or zero to eliminate it.

The solving algorithm used in flow model is PIMPLE, which is a combination of PISO (Pressure Implicit with Splitting of Operators) and SIMPLE (Semi-Implicit Method for Pressure-Linked Equations) algorithms. In this study, the $\kappa -ϵ$ RAS (Reynolds-averaged simulation) turbulence model is adopted to model the turbulent viscosity ${\nu }_{turb}$ as:

$${\nu }_{turb}=\frac{C_{\mu }{\kappa }^2}{ϵ}$$

(8)

where $C_{\mu }$ is an empirical constant; $\kappa$ is the turbulence kinetic energy while $ϵ$ is the turbulence energy dissipation rate, separately.

The IHFOAM implements the wave generation and active wave absorption in the fluid domain, which introduce several boundary conditions: (i) the inlet boundary condition allows the generation of wave according to different wave theories as well as adding different steady current flow; (ii) the outlet boundary condition applies an active wave absorption theory to prevent the re-reflection of incoming wave; (iii) the slip boundary condition (zero-gradient) is applied on the bottom of the fluid domain and the lateral boundary of the numerical wave flume; (iv) the top boundary condition is set as the atmospheric pressure. The details of IHFOAM could be found in Higuera et al. [37].

2.2. Seabed Model

It is widely known that saturated soil considered as a multi-phase material is formed by soil particles and the void of the skeleton. The pores of the solid phase are filled with the water and trapped air distributed through the body. Therefore, in order to simulate the interaction between the soil skeleton and pore water in a porous seabed, a seabed model is established based on the partially dynamic Biot’s equation ( also known as “$u-p$” approximation) [45] with consideration of acceleration inertia term of flow and soil particles.

With the assumptions of a homogeneous and isotropic seabed and compressible pore fluid, the mass conservation equation of pore fluid can be expressed as [46]:

$$K_s{\nabla }^2{p}_s-{\gamma }_w{n}_s{\beta }_s\frac{\partial p_s}{\partial t}+K_s{\rho }_f\frac{{\partial }^2{ϵ}_s}{\partial t^2}={\gamma }_w\frac{\partial {ϵ}_s}{\partial t},$$

(9)

where $p_s$ is pore water pressure, ${\gamma }_w$ is the unit weight of water, $K_s$ is the soil permeability, $n_s$ is soil porosity; ${\beta }_s$ is the compressibility of pore fluid while ${ϵ}_s$ is volume strain, which can be expressed as:

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

(10)

where $u_s$ and $w_s$ are the soil displacements in the x- and z- direction, respectively; $K_w$ is the bulk modulus of pore fluid ($K_w=1.95\times {10}^9\mathrm{N}/{\mathrm{m}}^2$ [4]); $S_r$ is the degree of saturation and $P_{wo}$ is absolute static water pressure, which defined as $P_{wo}={\gamma }_{w}d$, in which d is the water depth.

Based on Newton’s second law, the force equilibrium equation in a poro-elastic medium in horizontal and vertical directions can be given as:

$$\frac{\partial {\sigma }_x^{\prime }}{\partial x}+\frac{\partial {\tau }_{xz}}{\partial z}=-\frac{\partial p_s}{\partial x}+\rho \frac{{\partial }^2{u}_s}{\partial t^2},$$

(11)

$$\frac{\partial {\tau }_{xz}}{\partial x}+\frac{\partial {\sigma }_z^{\prime }}{\partial z}=-\frac{\partial p_s}{\partial z}+\rho \frac{{\partial }^2{w}_s}{\partial t^2},$$

(12)

where ${\sigma }_x^{\prime }$ and ${\sigma }_z^{\prime }$ are the effective normal stresses in horizontal and vertical direction respectively; ${\tau }_{xz}$ is shear stress component; $\rho$ is the average density of a porous seabed and can be obtained by $\rho ={\rho }_f{n}_s+{\rho }_s(1-n_s)$, in which ${\rho }_f$ is the fluid density while ${\rho }_s$ is the solid density.

Based on the pore-elastic theory, the effective normal stresses and shear stress can be expressed in term of soil displacements:

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

(13)

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

(14)

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

(15)

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

Substituting (13)–(15) into (11)–(12), we have the governing equations for force balance as

$$G{\nabla }^2{u}_s+\frac{G}{1-2{\nu }_s}\frac{\partial }{\partial x}(\frac{\partial u_s}{\partial x}+\frac{\partial w_s}{\partial z})=\frac{\partial p_s}{\partial x}+{\rho }_s\frac{{\partial }^2{u}_s}{\partial t^2},$$

(16)

$$G{\nabla }^2{w}_s+\frac{G}{1-2{\nu }_s}\frac{\partial }{\partial z}(\frac{\partial u_s}{\partial x}+\frac{\partial w_s}{\partial z})=\frac{\partial p_s}{\partial z}+{\rho }_s\frac{{\partial }^2{w}_s}{\partial t^2},$$

(17)

2.3. Boundary Conditions

At seabed surface ($z=0$), the vertical effective stress and shear stress is assumed to be zero and the water pressure is directly acting on it.

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

(18)

where $P_b$ is the dynamic wave pressure at the seabed surface, which is obtained by the wave model.

For the soil resting on the seabed bottom and lateral boundaries, zero displacements and impermeable are considered at the seabed bottom ($z=-h$), i.e.,

$$u_s=w_s=0,\mathrm{and}\frac{\partial p_s}{\partial z}=0,$$

(19)

In addition, the boundary condition of tunnel surface is assumed to be impermeable with zero displacements, i.e.,

$$u_s=w_s=0,\mathrm{and}\frac{\partial p_s}{\partial n_v}=0,$$

(20)

where $n_v$ denotes normal vector of the tunnel surface.

2.4. Meshfree Model for the Seabed Domain

The LRBFCM is adopted to solve the governing equations listed above. First, an approximation function $\mathsf{\Phi }\left(y_i\right)$ which can either stand for displacement or pore pressure in the “$u-p$” formulation is considered in the computational geometry. This function is composed of an arbitrarily distributed points series $y_j(j=1,2,\cdots ,n)$ located both in the computational domain and on its boundary. Approximate $\mathsf{\Phi }=\mathsf{\Phi }\left(x\right)$ around $y_i$ by the RBF $\chi \left(r_m\right)$ to construct a linear equation for each node $y_n$ as,

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

(21)

where ${\alpha }_m$ is the undetermined coefficient for $\chi \left(r_m\right)$, and $\chi \left(r_m\right)$ is the MQ defined by,

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

(22)

with $r_m$, c, and $x_m$ being the Euclidean distance from x to $x_m$, the shape parameter [30], and the positions of the K nearest neighbour nodes around the prescribed centre $x_1=y_n$, respectively. An algorithm based on the kd-tree is adopted to search the K nearest neighbour nodes [47].

Then, (21) is collocated on the K nearest neighbour nodes, which can be expressed as:

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

(23)

$${\left[\mathsf{\Phi }\right]}_{K\times 1}={\left[\mathsf{\Phi }\left(x_1\right),\mathsf{\Phi }\left(x_2\right),\cdots ,\mathsf{\Phi }\left(x_K\right)\right]}^T,$$

(24)

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

(25)

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

(26)

Inversion (23) gives,

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

(27)

Here, we assume a linear differential operator of both governing equation and the boundary condition which represented by L. If collocated the $L\mathsf{\Phi }\left(x\right)$ on $x_1=y_n$, (23) can be written as:

$$L\mathsf{\Phi }\left(y_n\right)=\sum _{m=1}^K{\alpha }_{m}L\chi \left(r_m\right){|}_{x=x_1}$$

(28)

which can be expressed as the vector form:

$$L\mathsf{\Phi }\left(y_n\right)={\left[L\chi \right]}_{1\times K}{\left[\alpha \right]}_{K\times 1}$$

(29)

In (29) mentioned above, the $L\chi \left(r_m\right)$ is on behalf of the results of the differential operator act on the the RBF $\chi \left(r_m\right)$. Then, by combining (29) and (27), it can be obtained as:

$$L\mathsf{\Phi }\left(y_n\right)={\left[C\right]}_{1\times K}{\left[\varphi \right]}_{K\times 1}$$

(30)

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

(31)

$${\left[L\chi \right]}_{1\times K}=\left[\begin{array}{cccc}L\chi \left(r_1\right){|}_{x=x_1}& L\chi \left(r_2\right){|}_{x=x_2}& \cdots & L\chi \left(r_K\right){|}_{x=x_K}\end{array}\right]$$

(32)

Obviously, the value of the row vector ${\left[C\right]}_{1\times K}$ can be obtained if all variables L, $\chi$ and $x_j$ are known.

For all $y_n$ in computational domain, the linear equations can be obtained by the above-mentioned localization procedure on either the governing equation or boundary conditions. Next, these equations can be assembled into the system matrix, as:

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

(33)

where

$$\begin{array}{cc}\hfill {\left[\mathsf{\Phi }\right]}_{N\times 1}=& \left[u_1\left(y_1\right),u_1\left(y_2\right),\cdots ,u_1\left(y_N\right),u_2\left(y_1\right),u_2\left(y_2\right),\cdots ,u_2\left(y_N\right),\right.\hfill \\ \hfill & {\left.p\left(y_1\right),p\left(y_2\right),\cdots ,p\left(y_N\right)\right]}^T,\hfill \end{array}$$

(34)

which, ${\left[\varphi \right]}_{N\times 1}$ is the sought solution, while ${\left[B\right]}_{N\times 1}$ is a column vector contributed from the external loadings. It is noticeable that (33) is a sparse system matrix which is similar to both finite difference method as well as the finite element method. In the present numerical model, the direct solver of SuperLU [48] is adopted to solve the sparse system, (33). The procedure of the LRBFCM is finished here.

In numerical strategy, it is necessary to integrate the governing equations in the time domain if the boundary conditions or the extrernal loadings are time-dependent. In the present model, the single-step time integration method of the Newmark method [49] is adopted, which could handle each time step independently when the first- and second-order time derivatives exist at the same time. More details can be found in some previous studies [50,51].

2.5. Integration Procedure of Flow Model and Seabed Model

For the integration procedure in this study, the one-way coupling algorithm is adopted for pore pressure delivery from the wave model to the seabed model. As demonstrated in Figure 1, the flow model is responsible to simulate the wave–current interaction which includes wave generation and fluid propagating according to the given wave/current characteristics. The water pressure will be extracted from the results after solving the VARANS equations, then used as the boundary conditions of seabed surface as the external loading in the seabed model, while the soil model can carry out the seabed behaviour under the external wave loading combined with the input parameters of seabed model. The displacement, the pore water pressure and the effective stresses could be determined by solving the Biot’s “$u-p$” approximation.

2.6. Convergence Tests

Before conducting parametric studies of the wave-induced dynamic response in a porous seabed under wave loading, it is necessary to check the convergence of a newly proposed numerical model. The convergence tests are carried out in regards of the nodes distance size ($\Delta x$), shape factor (c) and the node number of the local region (K), which could have an influence on the numerical accuracy and computational efficiency.

Firstly, the small node distance size makes the results more accurate, however, it will result in enormous computational cost. As shown in Figure 2, the $\Delta x$ is equal to $L/50$, $L/100$, and $L/200$ respectively (L is the wavelength). The non-dimensional pore water pressures ($p_s/p_0$) are depicted, $p_0$ represents the amplitude of linear wave pressure at the seabed surface. From the figure, the result for the case of $\Delta x$ equal to $L/50$ is slightly difference from the others, while the results are almost the same for $\Delta x$ equal to $L/100$ and $L/200$, which indicate the model is convergent with a node distance that is smaller than $L/100$.

Next, the convergence test of shape factor is conducted.The shape factor is generally equal to 15–60 times the maximum Euclidean distance between two adjacent nodes by convention. The c adopts 3.975 (15$\times \Delta x$), 7.95 (30$\times \Delta x$), and 15.9 (60$\times \Delta x$) separately. From the Figure 3, it can be conclude that the numerical results are not sensitive to the affect the shape factor on the consequence of the almost the same results obtained for this set of shape factors.

Moreover, the effect of the number of nearest neighbour nodes in a local region is examined. The dynamic pore water responses in a seabed are depicted in Figure 4 regards for the different numbers of neighbour nodes in the local region, which are 5, 9 and 13, respectively. From the figure, the result shows a good tendency during the whole process when $K=9$. It also can be seen that for $K=5$, the amplitude of the result is slightly different with the result of $K=9$. When $K=13$, the value of $|p_s|/p_0$ is even beyond 1 after 2 s, which is obviously wrong. This condition might be the ill-condition in this case. Thus, the number of the nodes located in the local region is 9 in this study.

3. Verification of the Proposed Model

Before conducting the parametric study, it is essential to verify the capability of the newly developed meshfree model. In this section, the proposed model is verified by comparison with analytical solutions and a set of published laboratory experiments from a previous study.

3.1. Comparison of the Present Model and the Analytical Solution for Wave-Seabed Interactions

To verify the proposed meshfree model under the dynamic condition, the present model is used to simulate the saturated isotropic seabed under linear wave loading based on a partially dynamic model. The result of the numerical solution are compared with the analytical solution presented by Hsu and Jeng [7]. The linear wave loading is applied to the seabed surface. The parameters for the comparison are given as: wave period $T=10$ s; water depth $d=20$ m; degree of saturation $S_r=0.975$; Poisson’s ratio ${\nu }_s=0.4$; soil porosity $n_s=0.35$; soil permeability $K_s={10}^{-2}$ m/s; shear modulus $G=5\times {10}^6$ Pa.

The vertical distributions of the wave-induced pore pressure ($|p_s|/p_0$), effective stresses ($|{\sigma }_z^{\prime }|/p_0$) and shear stress ($|{\tau }_{xz}|/p_0$) versus the soil depth ($z/H$) in a porous seabed are plotted in Figure 5. As depicted in the figure, the blue lines represent the results obtained from analytical solution, while red circles denote soil behaviour simulated from the present “$u-p$” time-dependent model. The figure clearly shows that the result obtained by the numerical method are in a good agreement with that of analytical solution, which demonstrates the capability of meshfree model for seabed dynamic response simulation.

3.2. Comparison of the Present Model with the Laboratory Experiments and Fem Results for Combined Waves and Current Loading

It is necessary to verify the performance of the integration model including both fluid and soil models under the circumstance of complex nature loading. There are numerous laboratory experiments related to the seabed response around the marine structures to date. However, it is quite limited in terms of the experimental data available for the case of immersed tunnel. Thus, the verifications of the integration model are carried out by comparison with the laboratory experiments and the FEM (Finite Element Method) results from DIANA-SWANDYNE II [52] for the seabed without the structure instead in this section. Qi and Gao [38] conducted a series of flume tests considering wave and wave combined currents as dynamic loading, respectively.

The first validation of this section is compared to the laboratory experiments conducted under wave loading only [38]. The input data for the first validation are: wave height $H=0.12$ m, water depth $d=0.5$ m, wave period $T=1.4$ s, seabed thickness $h=1.2$ m, degree of saturation $S_r=1.0$, shear modulus $G={10}^7$ N/m${}^2$, Poisson’s ratio ${\nu }_s=0.3$, permeability $K_s=1.88\times {10}^{-4}$ m/s. Figure 6 depicts the wave patterns with corresponding dynamic pore water pressure of the seabed, which are predicted by the present model and obtained from the experiment and FEM model separately. It can be seen that the result obtained from both the wave model and seabed model are in good agreement with the test data, which indicate that the present model is capable for simulating the wave motion in the fluid domain as well as the corresponding soil response of a sandy seabed.

The second validation in this section is to compare with the previous laboratory experiments conducted by Qi and Gao [38]. Unlike the previous case, this test simulates the seabed dynamic response under wave and current, which are generated synchronously. The following current with velocity of 0.05 m/s is adopted. The wave parameters and the soil properties are the same as above. As shown in Figure 7, the fluid pattern tracked by the fluid model matches well with the experiment data, while the pore water pressure simulated by the present model in correspondence with that obtained from the experiment and the FEM model. Thus, the current model performs well for simulating a more realistic marine dynamic elastic behaviour including both the fluid and soil parts.

4. Dynamic Response of the Seabed

In this study, a new meshfree model is developed, based on Biot’s “$u-p$” approximation to simulate the dynamic sandy seabed behaviour around an immersed tunnel under complex natural loading including wave and current. As shown in Figure 8, the external loadings are assumed to be propagating over the porous seabed in which the immersed tunnel is buried. The buried depth is defined as b adopting 0.5 m below the seabed surface in this case. The sandy seabed foundation is treated as an elastic two-phase medium above a rigid impermeable bottom with 200 m for seabed length ($Lx$) and 40 m for seabed thickness (h). The seawater depth is specified as d. The immersed tunnel is assumed to be placed on the trench dredged on the middle of the seabed of the computational domain. The trench is back-filled by the same type of loose sand with the seabed soil. The tunnel geometry, wave profile and seabed profile in this case are roughly the same as the actual conditions of the Hong Kong–Zhuhai–Macao bridge tunnel, which could provide a reference of the sandy seabed dynamic for such a large immersed tunnel under wave loading.

The detailed dimensions of the immersed tunnel are given in the Figure 9. As shown in figure, the immersed tunnel in this study is considered as an elastic material comprising two traffic tubes of 30 m long and 9 m high (cross section). The boundary condition of the immersed tunnel is treated as impermeable with zero pore pressure gradient. No relative displacement is assumed between the seabed soil and the tunnel frame on the consequence of the high fraction exists between the concrete tunnel surface and seabed soil.

The configuration of the fluid domain and seabed domain can be found in

Table 1. Parameters used for standard case in parametric study.

Wave Characteristics	Value	Unit
Wave period (T)	10	s
Wave height (H)	4	m
Wavelength (L)	142.355	m
Water depth (d)	30	m
Soil Properties	Value	Unit
Poisson’s ratio (${\nu }_s$)	0.20	-
Permeability ($K_s$)	$1.0\times {10}^{-5}$	m/s
Porosity ($n_s$)	0.42	-
Degree of saturation ($S_r$)	0.95	-
Shear modulus (G)	$5.0\times {10}^6$	Pa
Density of soil particles (${\rho }_s$)	2650	kg/m${}^3$
Tunnel buried depth (b)	0.5	m
Modelling Parameters	Value	Unit
Number of KNN (K)	9	-
Shape factor (c)	6	-
Nodes number (x direction)	2000	-
Nodes number (z direction)	400	-
$\Delta t$	0.5	s
Time steps in a period ($Ndt$)	20	-

, as well as the wave characteristics, seabed soil properties and modelling parameters. The seabed foundation is assumed to be composed of relative dense deposited sand [53]. Figure 10 shows the applicability range of the different wave theories [54]. The wave characteristics adopted in this study are in the range of Stokes second-order wave ($H=$ 4 m, $T=$ 10 s, $d=$ 30 m, shown as a red star in Figure 10), which is generated and simulated by fluid model. The node number of the local region (K) for local RBF method is 9, while a positive constant c which is known as the shape factor equals 6. The total number of nodes in this case is 771,140. The convergence test has been done to check the stability of the present parameter configuration, which is quite enough to obtain an accuracy and detailed result. The time step $\Delta t$ set in this case is 0.5 s, while 20 time steps are contained in one wave period.

The aim of this section is tracking the dynamic soil response during the wave propagating over the seabed around the immersed tunnel. During the wave propagation from the left to the right of the computational domain continuously, the effective stresses and pore water pressures show a correlation trend with the change of the water pressure acting on the surface of the seabed. As shown in Figure 11, the oscillatory wave-induced pore water pressures, horizontal displacements ($u_s$) and vertical displacement ($w_s$) for the computational domain of the seabed at $t=$ 13 s are presented, respectively. It is found that the seabed dynamic behaviour with immersed tunnel is not periodic symmetry any more under the cyclic wave loading. The figure shows that the existence of the immersed tunnel has an obvious influence on the dynamic behaviour of the sandy seabed soil nearby by comparing the region located on the leftward and the rightward of the tunnel. It can be seen in Figure 11 that the placement of a tunnel weakens the displacement change of local area in a way, while the fluctuation of the dynamic pore water pressure decrease around as well. In Figure 11c, the dynamic pore water pressure of the seabed soil beneath the tunnel bottom shows a different tendency from the surrounding that the positive oscillatory pore pressure occurs on the left corner while the negative occurs on the right corner. In order to figure out a more detailed dynamic soil response in the vicinity of the immersed tunnel, the results of dynamic pore water pressure of some typical locations are analysed.

As shown in Figure 8, points A to D are a set of symmetrical nodes about the tunnel at the depth of −5 m in the seabed ($x=$ 60, 82, 118, 140 m, respectively), while points E to F located at the depth of −10 m ($x=$ 60, 85, 115, 140 m respectively). Points F and G are set on 0.5 m below the two base corners of the immersed tunnel. Figure 12 depicts the time series of the dynamic pore water pressure generated by the two point sets. From the figure, it can be seen that the vibration of the pore water pressure is more violent on the remote seabed from the tunnel, which is consistent with the conclusion mentioned above. Furthermore, the amplitude reduced for point F and G below the tunnel are slightly less than the points B and C on two sides of the tunnel. The pore water pressure vibration of the soil on points F and G is even larger than on points B and C, which indicates that the seabed foundation beneath the immersed tunnel is more likely to be unstable due to transient liquefaction. In addition, a phenomenon occurred which is that there is a phase difference of dynamic pore pressure brought out under the base corners of the tunnel (points F and G), as shown in Figure 12b.

5. Wave-Induced Liquefaction

The stability of the coastal structures and its seabed foundation is one of the main concerns for engineering design procedure. The wave-induced liquefaction in a porous seabed is one of the most significant unstable factors. Zen and Yamazaki [55] pointed out that the liquefaction of the porous seabed is responding to the variation of the ocean wave, which is actually caused by the periodic upward seepage force. Thus, we proposed to estimate the liquefaction state in two-dimensions, i.e.,

$${\sigma }_0^{\prime }+(P_b-p_s)\le 0$$

(35)

where ${\sigma }_0^{\prime }$ represents the initial effective stress, $P_b$ is the wave pressure acting on seabed, while $p_s$ is the wave-induced transient pore pressure. The value of $P_b-p_s$ is equal to the excess pore pressure generated by the wave loading ($|u_e|$).

Figure 13 shows the wave-induced transient liquefaction area around the immersed tunnel at three typical times ($t=12$ s, $t=$ 14 s and $t=$ 15.5 s) separately. As shown in the figure, the transient liquefaction area moves along the direction of the wave propagation. The previously liquefied area is able to recover as the wave trough go away. This process is repeated periodically under the cyclic wave loading. The maximum liquefaction depth in this case is 0.8 m below the seabed surface, as seen in Figure 13a, while the soil that covered the tunnel is fully liquefied during one wave period which illustrated in Figure 13b. Thus, the back filling soil above the tunnel can not protect the immersed tunnel any more in this circumstance. Moreover, the maximum liquefaction depth of the rightward seabed of the tunnel is 0.6 m, which is slightly shallow than that in leftwards of 0.8 m.

6. Parametric Study

In this section, parametric studies of the seabed liquefaction around the immersed tunnel are carried out, which include the wave characteristics, soil properties and the effects of current specifically.

6.1. Effects of Wave Characteristics

Basically, the wave-induced oscillatory liquefaction phenomenon of the porous seabed is highly relative to the water pressure propagating on it. As a consequence, the relation of the wave characteristics and the liquefaction depth in the vicinity of the tunnel is discussed in detail with respect to the wave period (T), wave height (H), and still water depth (d), separately.

Table 2. Parameters used for in parametric study of wave characteristics.

Wave Characteristics	Value	Unit
Wave period (T)	8 or various	s
Wave height (H)	2 or various	m
Water depth (d)	30 or various	m
Soil Properties	Value	Unit
Poisson’s ratio (${\nu }_s$)	0.20	-
Permeability ($K_s$)	$1.0\times {10}^{-5}$	m/s
Porosity ($n_s$)	0.42	-
Degree of saturation ($S_r$)	0.95	-
Shear modulus (G)	$5.0\times {10}^6$	Pa
Tunnel buried depth (b)	0.5	m

lists the input data involving the wave variables and other parameters. Figure 14 characterizes the maximum liquefaction depth of the vertical section at $x=$ 80 m which is on the left side of the tunnel in regard of the different wave conditions.

It is well known that the wave height is in positive correlation to the value of the water pressure acting on the seabed precisely. Furthermore, the water pressure will affect the wave-induced pore water pressure in seabed deposits. As illustrated in Figure 14a, the maximum liquefaction depth grows deeper with an increase of the wave height, which shows the positive relationship between the wave height and the potential of liquefaction.

The wave period is another key wave parameter that influences the wave length directly, which will further affect the seabed transient liquefy process. In this study, the wave periods are from 8 to 10 s. From Figure 14b, the maximum liquefaction depth is progressively increase with the rise of wave periods. It can be concluded that the liquefaction degree is more intense in the case of a long wave period.

Lastly, the influence of still water depth is discussed. The water depth ranges from 30 to 40 m. Unlike wave height and period, the maximum oscillatory liquefaction depth decreases with the raising of still water depth, as shown in Figure 14c, which indicates that the wave-induced transient liquefaction is more likely to occur in the shallow water areas.

6.2. Effects of Soil Properties

Besides the wave characteristics, the deposit conditions are found to be essential in a wave-seabed structure interactions. The effects of two parameters of soil properties are discussed individually, i.e., degree of saturation ($S_r$), and soil permeability ($K_s$). A series of topical values are selected in this section, which are from $1.0\times {10}^{-5}$ to $1.0\times {10}^{-3}$ for soil permeability, and 0.925–0.975 for degree of saturation, respectively. The wave condition and other parameters used in this section can be found in

Table 3. Parameters used for in parametric study of soil properties.

Wave Characteristics	Value	Unit
Wave period (T)	10	s
Wave height (H)	4	m
Water depth (d)	30	m
Soil Properties	Value	Unit
Poisson’s ratio (${\nu }_s$)	0.20	-
Permeability ($K_s$)	$1.0\times {10}^{-5}$ or various	m/s
Porosity ($n_s$)	0.42	-
Degree of saturation ($S_r$)	0.925 or various	-
Shear modulus (G)	$5.0\times {10}^6$	Pa
Tunnel buried depth (b)	0.5	m

.

To examine how the different soil parameters affect the seabed liquefaction around the immersed tunnel, Figure 15 and Figure 16 are depicted to show the vertical distribution of the excess pore pressure ($|u_e|/{\sigma }_0^{\prime }$) and the maximum liquefaction depth of the vertical section at $x=$ 80 m in seabed with various soil properties. Figure 15b illustrates that the maximum liquefaction depth is larger in the seabed with relative low degree of saturation. While the same conclusion can be draw for the various of soil permeability in Figure 16b, i.e., the smaller the soil permeability adopted, the deeper the liquefaction occurs in the vicinity of the tunnel.

6.3. Effects of Current

In reality, the fluid circumstance above the seabed are quite complex, and need to considered in the wave–current interaction. The motion of the current is able to influence the wave propagation, which further affects the seabed liquefaction process. This section aims to investigate the relationship between seabed liquefaction and the current around the immersed tunnel. The seabed response under a second-order Stokes wave ($T=$ 10 s, $H=$ 4 m and $d=$ 35 m) with a series of following currents ($U_0=$ 0.5, 1.0, 1.5 m/s) and opposing currents ($U_0$ = −0.5, −1.0, −1.5 m/s) are compared with the no current case. Other parameters used in this study are listed in

Table 4. Parameters used in study of effects of current.

Wave Characteristics	Value	Unit
Wave period (T)	10	s
Wave height (H)	4	m
Water depth (d)	35	m
Current velocity ($U_0$)	−1.5 or various	m/s
Soil Properties	Value	Unit
Poisson’s ratio (${\nu }_s$)	0.20	-
Permeability ($K_s$)	$1.0\times {10}^{-5}$	m/s
Porosity ($n_s$)	0.42	-
Degree of saturation ($S_r$)	0.95	-
Shear modulus (G)	$5.0\times {10}^6$	Pa
Tunnel buried depth (b)	0.5	m

.

As seen in Figure 17, the maximum liquefied depth in seabed around the tunnel ($x=$ 80 m) is 0.4 m, 0.6 m and 0.7 m when the current velocity $U_0$ takes on −1.5 m/s, 0 m/s and 1.5 m/s, which indicate the following current could deeper the liquefaction zone while the opposing current could decrease the liquefaction depth. Moreover, Figure 18 shows the liquefaction area around the immersed tunnel (when wave trough travels above the cross section $x=80$ m) triggered by the wave combined different currents velocities. It can be concluded that not only the liquefied depth, but the liquefaction zone changes around the tunnel are also positively relative to the velocity of the current. Thus, oscillatory liquefaction is more likely to occur under the following current and the opposing current is able to decrease the potential of liquefaction.

7. Conclusions

In this study, an integrated numerical model including a LRBFCM seabed model based on Biot “$u-p$” approximation is proposed to investigate the flow field dynamics and corresponding seabed behaviour around an immersed tunnel under second-order stokes wave combined current with various velocities. The oscillatory liquefaction under the different wave characteristics, soil properties, and current velocities are discussed in detail. From the numerical result, the following conclusions can be drawn:

(1)

The newly developed meshfree model is well validated by comparison with the analytical solution, the laboratory experiments and previous numerical results. The LRBFCM is examined to be reliable in simulation of wave-induced oscillatory liquefaction behaviour of a seabed. Moreover, on the consequence of mesh-independence, the present model is more efficient than the conventional model based on FEM (Finite element method) or FVM (Finite volume method).

(2)

The existence of the immersed tunnel affects surrounding seabed dynamic behaviours, which are able to weaken the displacement and dynamic pore pressure change nearby. Furthermore, the maximum liquefied depth on the right of the tunnel is smaller than that on the left.

(3)

The wave-induced transient liquefaction around the tunnel is highly relative to the wave characteristics. It is found that the seabed liquefaction is more likely to occur under a shallow water area with the waves of large wave height and long period.

(4)

The parametric studies show that the soil properties around tunnel have a significant impact on the liquefaction behaviour as well. The smaller the soil permeability and degree of saturation adopted, the deeper the liquefaction occurs in the vicinity of the tunnel.

(5)

Based on the numerical result, the occurrence of currents could obviously affect wave-induced liquefaction. The following current can aggravate the seabed liquefaction while the opposing current can decrease the liquefied risk around the tunnel.