Modeling of dynamic response of poroelastic soil layers under wave loading

Abstract

In this paper, the dynamic response of saturated and layered soils under harmonic waves is modeled using the finite element method. The numerical results are then verified by corresponding analytical solutions which are also developed by the author. The equations governing the dynamics of porous media are written in their fully dynamic form and possible simplifications are introduced based on the presence of inertial terms associated with solid and fluid phases. The response variations are presented in terms of pore water pressure and shear stress distributions within the layers. It is determined that a set of non-dimensional parameters and their respective ratios as a result of layering play a major role in the dynamic response.

Keywords

dynamic response of soils / coupled flow-deformation / finite elements / analytical solution / harmonic waves

Introduction

In practical applications of Geotechnical Engineering, often times engineers come across a series of soil layers deposited beneath the foundations regardless of where they are located. In marine environment, this has particular importance in the response of seabed and seabed-structure systems considering the constant effect of the incoming cyclic waves causing similar alternating shear and normal stresses acting on these systems. The studies conducted thus far highlight mostly the response and instability of a single soil layer under various loads generally extrapolating the variation of unknowns to deeper soils or sometimes simply calculating a weighted average of some key soil variables. However, as initially emphasized by some studies such as [1] and later presented by others that a single layer approach cannot fully represent the actual mechanics of a multi-layer case even if some averaging takes place [2,3]. Jeng and his colleagues [4-6] developed various analytical solutions including the inertial effects and wave nonlinearity on the soil response. Solutions considering multiple soil layers to the response have been limited [1,7,8] which also neglected the inertial terms. Mesgouez and Mesgouez [9] developed a semi-analytical approach to the wave propagation problem in a multi-layer poro-visco-elastic system due to transient loads. Therefore, the dynamic response of a multi-layer soil system needs to be evaluated in terms of the distributions of internal variables such as pore water pressure within the layers. For that generalized poroelasticity equations developed originally by [10] who later included inertial terms [11,12] are used.

In this study, two dimensional plane strain finite element models are developed for the dynamic response of multi-layer seabed soil under harmonic progressive waves. The governing equations are presented in their fully dynamic form (FD) and possible simplifications of partially dynamic (PD) and quasi-static (QS) formulations written in terms of the negligence of inertial terms are also introduced. The finite element models are numerically solved considering a set of boundary and initial conditions and the results are verified with their analytical counterparts which are also developed by the author. Subsequently, dynamic response is evaluated and presented in terms of pore water pressure variations in the medium. The effect of surface layer with relatively low permeability along with the underlying more granular and permeable layers on the response is presented for various soil and wave parameters. The study of the effect of inertial terms revealed that unless the hydraulic conductivities of adjacent soil layers or the applied wave frequency is considerably high (which is not typical for the sea environment) the inertial terms are insignificant.

Governing equations

The three basic laws of: 1) Constitutive relation for the stress-strain relationship, 2) Momentum balance of solid and liquid phases and 3) Mass balance equation for the fluid phase, are necessary to develop a system of equations defining the dynamic response of soil. Soils can be viewed as two phase porous materials considering the solid grains with different size and shape, and a viscous fluid (generally water) allowed to flow through its pores. The air phase is assumed to be fully dissolved in the fluid constituting a single compressible fluid. This section briefly summarizes the governing equations proposed first by Biot [10-12].

Constitutive relationship

The effective stress is defined as:

(1){\sigma }_{ij}\prime ={\sigma }_{ij}+{\delta }_{ij}p,

where {\sigma }_{ij} is total stress, {\sigma }_{ij}\prime is effective stress, {\delta }_{ij} is Kronecker delta and p is the pore water pressure. Tensile stress is taken positive. Average strain, {ϵ}_{ij} is,

(2){ϵ}_{ij}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}(u_{i,j}+u_{j,i}),

with u_{i,j} being the derivatives of the displacement of soil skeleton with respect to spatial coordinates. The stress-strain relationship is defined as,

(3){\sigma }_{ij}\prime =D_{ijkl}{ϵ}_{kl},

where D_{ijkl} is the elastic material matrix.

Momentum balance equations

Total equilibrium of the system and the equilibrium of the fluid phase are written respectively as;

(4){\sigma }_{ij,j}+\rho g_i-{\rho }_f{\ddot{\overline{w}}}_i-\rho {\ddot{u}}_i=\mathrm{0},

(5)-p_{,i}+{\rho }_f{g}_i-{\displaystyle \frac{{\dot{\overline{w}}}_i}{k_i}}{\rho }_f{g}_i-{\rho }_f{\ddot{u}}_i-{\displaystyle \frac{{\rho }_f}{n}}{\ddot{\overline{w}}}_i=\mathrm{0.}

where ρ and {\rho }_{\mathrm{f}} are the densities of the mixture and fluid; {\ddot{u}}_i={\displaystyle \frac{{\mathrm{d}}^{\mathrm{2}}{u}_i}{\mathrm{d}{t}^{\mathrm{2}}}} is the acceleration of the soil skeleton, {\overline{w}}_i is the average displacement of the fluid phase relative to the solid and {\ddot{\overline{w}}}_i is the average acceleration of the fluid relative to that of solid which is equal to n{\ddot{w}}_i with porosity being n. Subscript indices indicate spatial derivatives. The Darcy’s law with k_i being the coefficient of hydraulic conductivity is introduced in Eq. (5). g_i is the body force acceleration.

Mass balance equation

The law of the conservation of mass yields,

(6){\dot{ϵ}}_{ii}+{\dot{\overline{w}}}_{i,i}+{\displaystyle \frac{n}{K_f}}\dot{p}=\mathrm{0},

where K_f is the bulk modulus of the pore fluid modified to account for slight unsaturation through the degree of saturation (S) where the pore fluid can still be treated as of single phase and the fluid compressibility, \overline{\beta }, Ref. [1,13] is defined using,

(7)\overline{\beta }={\displaystyle \frac{\mathrm{1}}{K_f}}={\displaystyle \frac{\mathrm{1}}{K_p}}+{\displaystyle \frac{\mathrm{1}-S}{p_{\mathrm{0}}}},

with p₀ being the absolute pore pressure and K_p, the bulk modulus of the fluid.

Various formulations

As indicated before, depending on the motions of the pore fluid and the solids, it is possible to have a number of idealizations (also referred to as formulations). These formulations are briefly summarized here.

Fully dynamic formulation

Considering the inertial forces associated with both solid and fluid phases, the fully dynamic (complete) formulation is summarized as,

(8a){\sigma }_{ij,j}+\rho g_i-\rho {\ddot{u}}_i-{\rho }_f{\ddot{\overline{w}}}_i=\mathrm{0},

(8b)-p_{,i}+{\rho }_{\mathrm{f}}{g}_i-{\rho }_{\mathrm{f}}{\ddot{u}}_i-{\displaystyle \frac{{\rho }_{\mathrm{f}}{\ddot{\overline{w}}}_i}{n}}-{\displaystyle \frac{{\rho }_{\mathrm{f}}g}{k_i}}{\dot{\overline{w}}}_i=\mathrm{0},

(8c){\dot{u}}_{i,i}+{\dot{\overline{w}}}_{i,i}+{\displaystyle \frac{n}{K_{\mathrm{f}}}}\dot{p}=\mathrm{0.}

The above set has the field variables: relative fluid displacement, w_i, solid displacement, u_i, and pore fluid pressure, p which results in the so called u-w-p formulation. The field variable, U_i, the total fluid displacement, is introduced into the equations in terms of,

(9)U_i=u_i+{\overline{w}}_i/n,

leading to a final set with two field variables, u_i and U_i providing inertial decoupling. Here g is the gravitational acceleration. In this new form, pore pressure is evaluated using the volumetric strains of individual phases,

(10)p=-Q(u_{i,i}(\mathrm{1}-n)+n{U}_{i,i}),

where Q=K_{\mathrm{f}}/n. In this paper this formulation is called the fully dynamic (FD) formulation.

Partially dynamic formulation

Neglecting the inertial terms associated with the relative pore fluid motion, the governing equations are reduced down to two with the field variables now being u_i and p. They take the form,

(11a){\sigma }_{ij,j}+\rho g_i=\rho {\ddot{u}}_i,

(11b){\dot{u}}_{i,i}+{\left({\displaystyle \frac{-k_i}{{\rho }_{\mathrm{f}}g}}{p}_{,i}+{\displaystyle \frac{k_i{g}_i}{g}}-{\displaystyle \frac{k_i}{g}}{\ddot{u}}_i\right)}_{,i}+{\displaystyle \frac{n}{K_{\mathrm{f}}}}{\dot{p}}_i=\mathrm{0.}

Quasi-static formulation

If all the inertial terms are neglected, the equations are now free from accelerations yielding,

(12a){\sigma }_{ij,j}+\rho g_i=\mathrm{0},

(12b){\dot{u}}_{i,i}+{\left({\displaystyle \frac{-k_i}{{\rho }_{\mathrm{f}}g}}{p}_{,i}+{\displaystyle \frac{k_i{g}_i}{g}}\right)}_{,i}+{\displaystyle \frac{n}{K_{\mathrm{f}}}}{\dot{p}}_i=\mathrm{0.}

Finite element modeling of soil layers under wave loading

A mathematical model to the dynamic response of a saturated multi-layer soil under harmonic waves is illustrated in Fig. 1. The governing equations described above for all formulations are discretized using classical finite elements. In obtaining the weak form, these equations are multiplied by corresponding virtual displacements, δu, δw, δp and δU, and integrated within the domain of interest. Discretization of the FD formulation leads to approximating the two field variables, u and U, in terms of finite element interpolation using,

(13a)\mathit{u}={\mathit{N}}_{\mathit{u}}\mathit{U},\dot{\mathit{u}}={\mathit{N}}_{\mathit{u}}\dot{\mathit{U}},\ddot{\mathit{u}}={\mathit{N}}_{\mathit{u}}\ddot{\mathit{U}},

(13b)\mathit{U}={\mathit{N}}_{\mathit{U}}{\mathit{U}}_{\mathit{t}},\dot{\mathit{U}}={\mathit{N}}_{\mathit{U}}{\dot{\mathit{U}}}_{\mathit{t}},\ddot{\mathit{U}}={\mathit{N}}_{\mathit{U}}{\ddot{\mathit{U}}}_{\mathit{t}},

(13c)\delta \mathit{u}={\mathit{N}}_{\mathit{u}},\delta \mathit{w}={\mathit{N}}_{\mathit{w}},\delta \mathit{U}={\mathit{N}}_{\mathit{U}},\delta \mathit{p}={\mathit{N}}_{\mathit{p}},

where N_u, N_w, N_p and N_U are the appropriate shape functions. The weak form of the equations then becomes,

(14a)$$\begin{gathered} \left({\displaystyle \underset{\mathrm{\Omega }}{\int }{(\nabla N_u)}^{\mathrm{T}}\left[{(\mathrm{1}-n)}^{\mathrm{2}}Q+\underset{˜}{\underset{˜}{D}}\right]\nabla N_u\mathrm{d}\mathrm{\Omega }}\right)U+\left({\displaystyle \underset{\mathrm{\Omega }}{\int }{(\nabla N_u)}^{\mathrm{T}}(n(\mathrm{1}-n)Q)\nabla N_U\mathrm{d}\mathrm{\Omega }}\right)U_t+\left({\displaystyle \underset{\mathrm{\Omega }}{\int }{(N_u)}^{\mathrm{T}}((\mathrm{1}-n){\rho }_s)N_u\mathrm{d}\mathrm{\Omega }}\right)\ddot{U} \\ -\left\{{\displaystyle \underset{\mathrm{\Omega }}{\int }{(N_u)}^{\mathrm{T}}[{\rho }_{\mathrm{f}}g{n}^{\mathrm{2}}{k}^{-\mathrm{1}}]N_U\mathrm{d}\mathrm{\Omega }}\right\}{\dot{U}}_t+\left\{{\displaystyle \underset{\mathrm{\Omega }}{\int }{(N_u)}^{\mathrm{T}}[{\rho }_{\mathrm{f}}g{n}^{\mathrm{2}}{k}^{-\mathrm{1}}]N_u\mathrm{d}\mathrm{\Omega }}\right\}\dot{U}={\displaystyle \underset{\mathrm{\Gamma }}{\int }{(N_u)}^{\mathrm{T}}(\sigma \prime )\mathrm{d}\mathrm{\Gamma }}-{\displaystyle \underset{\mathrm{\Gamma }}{\int }{(N_u)}^{\mathrm{T}}n(\mathrm{1}-n)p\mathrm{d}\mathrm{\Gamma }} \\ +{\displaystyle \underset{\mathrm{\Omega }}{\int }{(N_u)}^{\mathrm{T}}(\mathrm{1}-n){\rho }_{s}g\mathrm{d}\mathrm{\Omega }}, \end{gathered}$$

(14b)$$\begin{gathered} \left({\displaystyle \underset{\mathrm{\Omega }}{\int }{(\nabla N_U)}^{\mathrm{T}}}(n(\mathrm{1}-n)Q)\nabla N_u\mathrm{d}\mathrm{\Omega }\right)U+\left({\underset{\mathrm{\Omega }}{\int }(\nabla N_U)}^{\mathrm{T}}(n^{\mathrm{2}}Q)\nabla N_U\mathrm{d}\mathrm{\Omega }\right)U_t+\left({\displaystyle \underset{\mathrm{\Omega }}{\int }{(N_U)}^{\mathrm{T}}}(n{\rho }_{\mathrm{f}})N_U\mathrm{d}\mathrm{\Omega }\right){\ddot{U}}_t \\ +\left({\displaystyle \underset{\mathrm{\Omega }}{\int }{(N_U)}^{\mathrm{T}}}({\rho }_{\mathrm{f}}g{n}^{\mathrm{2}})k^{-\mathrm{1}}{N}_U\mathrm{d}\mathrm{\Omega }\right){\dot{U}}_t-\left({\displaystyle \underset{\mathrm{\Omega }}{\int }{(N_U)}^{\mathrm{T}}}({\rho }_{\mathrm{f}}g{n}^{\mathrm{2}})k^{-\mathrm{1}}{N}_u\mathrm{d}\mathrm{\Omega }\right)\dot{U}=-{\displaystyle \underset{\mathrm{\Gamma }}{\int }{(N_U)}^{\mathrm{T}}n(n)p}\mathrm{d}\mathrm{\Gamma } \\ +{\displaystyle \underset{\mathrm{\Omega }}{\int }{(N_U)}^{\mathrm{T}}n{\rho }_{\mathrm{f}}g\mathrm{d}\mathrm{\Omega }}, \end{gathered}$$

where \underset{\approx }{D} is the constitutive matrix and \nabla N_u is the gradient of the shape function matrix. This can be written in a matrix form as,

(15)\left[\begin{array}{cc}{\mathbf{M}}_{\mathbf{u}}& \mathbf{0}\\ \mathbf{0}& {\mathbf{M}}_{\mathbf{U}}\end{array}\right]\left\{\begin{array}{c}\ddot{\mathbf{U}}\\ {\ddot{\mathbf{U}}}_{\mathbf{t}}\end{array}\right\}+\left[\begin{array}{cc}{\mathbf{C}}_{\mathbf{u}}& -{\mathbf{C}}_{\mathbf{u}\mathbf{U}}\\ -{\mathbf{C}}_{\mathbf{u}\mathbf{U}}^{\mathbf{T}}& {\mathbf{C}}_{\mathbf{U}}\end{array}\right]\left\{\begin{array}{c}\dot{\mathbf{U}}\\ {\dot{\mathbf{U}}}_{\mathbf{t}}\end{array}\right\}+\left[\begin{array}{cc}{\mathbf{K}}_{\mathbf{u}}& {\mathbf{K}}_{\mathbf{u}\mathbf{U}}\\ {\mathbf{K}}_{\mathbf{u}\mathbf{U}}^{\mathbf{T}}& {\mathbf{K}}_{\mathbf{U}}\end{array}\right]\left\{\begin{array}{c}\mathbf{U}\\ {\mathbf{U}}_{\mathbf{t}}\end{array}\right\}=\left\{\begin{array}{c}{\mathbf{f}}_{\mathbf{u}}\\ {\mathbf{f}}_{\mathbf{U}}\end{array}\right\},

where

K_u={\displaystyle \underset{\mathrm{\Omega }}{\int }{(\nabla N_u)}^{\mathrm{T}}}{(\mathrm{1}-n)}^{\mathrm{2}}Q\nabla N_u\mathrm{d}\mathrm{\Omega }+{\displaystyle \underset{\mathrm{\Omega }}{\int }{(\nabla N_u)}^{\mathrm{T}}}[D]\nabla N_u\mathrm{d}\mathrm{\Omega },K_{uU}={\displaystyle \underset{\mathrm{\Omega }}{\int }{(\nabla N_u)}^{\mathrm{T}}}n(\mathrm{1}-n)Q\nabla N_U\mathrm{d}\mathrm{\Omega },

K_U={\underset{\mathrm{\Omega }}{\int }(\nabla N_U)}^{\mathrm{T}}(n^{\mathrm{2}}Q)\nabla N_U\mathrm{d}\mathrm{\Omega },C_{uU}={\displaystyle \underset{\mathrm{\Omega }}{\int }{(N_u)}^{\mathrm{T}}}({\rho }_{\mathrm{f}}g{n}^{\mathrm{2}}){[k]}^{-\mathrm{1}}{N}_U\mathrm{d}\mathrm{\Omega },

C_u={\displaystyle \underset{\mathrm{\Omega }}{\int }{(N_u)}^{\mathrm{T}}({\rho }_{\mathrm{f}}g{n}^{\mathrm{2}}){[k]}^{-\mathrm{1}}{N}_u\mathrm{d}\mathrm{\Omega }},C_U={\displaystyle \underset{\mathrm{\Omega }}{\int }{(N_U)}^{\mathrm{T}}}({\rho }_{\mathrm{f}}g{n}^{\mathrm{2}}){[k]}^{-\mathrm{1}}{N}_U\mathrm{d}\mathrm{\Omega },

M_u={\displaystyle \underset{\mathrm{\Omega }}{\int }{(N_u)}^{\mathrm{T}}}(\mathrm{1}-n){\rho }_s{N}_u\mathrm{d}\mathrm{\Omega },C_U={\displaystyle \underset{\mathrm{\Omega }}{\int }{(N_U)}^{\mathrm{T}}}({\rho }_{\mathrm{f}}g{n}^{\mathrm{2}}){[k]}^{-\mathrm{1}}{N}_U\mathrm{d}\mathrm{\Omega },

\{f_u\}={\displaystyle \underset{\mathrm{\Gamma }}{\int }{(N_u)}^{\mathrm{T}}\sigma \prime \mathrm{d}\mathrm{\Gamma }}-{\displaystyle \underset{\mathrm{\Gamma }}{\int }{(N_u)}^{\mathrm{T}}n(\mathrm{1}-n)p\mathrm{d}\mathrm{\Gamma }}+{\displaystyle \underset{\mathrm{\Omega }}{\int }{(N_u)}^{\mathrm{T}}(\mathrm{1}-n){\rho }_{s}g\mathrm{d}\mathrm{\Omega }},

\{f_U\}=-{\displaystyle \underset{\mathrm{\Gamma }}{\int }{(N_U)}^{\mathrm{T}}n(n)p\mathrm{d}\mathrm{\Gamma }}+{\displaystyle \underset{\mathrm{\Omega }}{\int }{(N_U)}^{\mathrm{T}}n{\rho }_{\mathrm{f}}g\mathrm{d}\mathrm{\Omega }}\mathrm{.}

The dynamic Eq. (15) can further be written as:

(16)\mathbf{M}\ddot{X}+\mathrm{C}\dot{X}+\mathbf{K}X=F,

where X is the unknown vector of field variables, M is the total mass matrix, C is the total damping matrix, K is the total stiffness matrix and F is the force vector. Spatial integration is carried out using the Gauss-Quadrature method through reduced integration points (i.e., one point integration of four-noded quadrilaterals and four point integrations for eight-noded quadrilaterals). Integration in time is carried out by using the implicit Newmark Method. Thus, at time step n + 1 with an increment Δt, we have the following equations for the acceleration, velocity and displacement values respectively,

(17a){\ddot{X}}_{n+\mathrm{1}}={\ddot{X}}_n+\delta {\ddot{X}}_{n+\mathrm{1}},

(17b){\dot{X}}_{n+\mathrm{1}}={\dot{X}}_n+\mathrm{\Delta }t\left[(\mathrm{1}-\gamma ){\ddot{X}}_n+\gamma {\ddot{X}}_{n+\mathrm{1}}\right],

(17c)X_{n+\mathrm{1}}=X_n+\mathrm{\Delta }t{\dot{X}}_n+{\displaystyle \frac{\mathrm{\Delta }{t}^{\mathrm{2}}}{\mathrm{2}}}\left[(\mathrm{1}-\mathrm{2}\beta ){\ddot{X}}_n+\mathrm{2}\beta {\ddot{X}}_{n+\mathrm{1}}\right].

Here, γ and β are the Newmark parameters controlling the stability and convergence of the numerical solution. The finite element formulations of the simplified PD and QS forms follow the same steps described above. However, the field variables will then be u and p for both and the final matrix system for PD formulation becomes,

(18)\left[\begin{array}{cc}{\mathbf{M}}_{\mathrm{s}}& 0\\ {\mathbf{M}}_{\mathrm{s}\mathrm{p}}& 0\end{array}\right]\left\{\begin{array}{c}\ddot{\mathbf{U}}\\ \ddot{\mathbf{P}}\end{array}\right\}+\left[\begin{array}{cc}0& 0\\ {\mathbf{C}}^{\mathrm{T}}& {\mathbf{C}}_{\mathrm{p}}\end{array}\right]\left\{\begin{array}{c}\dot{\mathbf{U}}\\ \dot{\mathbf{P}}\end{array}\right\}+\left[\begin{array}{cc}{\mathbf{K}}_{\mathrm{s}}& -\mathbf{C}\\ 0& {\mathbf{K}}_{\mathrm{p}}\end{array}\right]\left\{\begin{array}{c}\mathbf{U}\\ \mathbf{P}\end{array}\right\}=\left\{\begin{array}{c}{f}_1\\ f_2\end{array}\right\}.

For QS form, only the mass matrix is neglected in the above equation and it becomes,

(19)\left[\begin{array}{cc}0& 0\\ {\mathbf{C}}^{\mathrm{T}}& {\mathbf{C}}_{\mathrm{p}}\end{array}\right]\left\{\begin{array}{c}\dot{\mathbf{U}}\\ \dot{\mathbf{P}}\end{array}\right\}+\left[\begin{array}{cc}{\mathbf{K}}_{\mathrm{s}}& -\mathbf{C}\\ 0& {\mathbf{K}}_{\mathrm{p}}\end{array}\right]\left\{\begin{array}{c}\mathbf{U}\\ \mathbf{P}\end{array}\right\}=\left\{\begin{array}{c}{f}_1\\ f_2\end{array}\right\}.

Here,

[K_s]={\displaystyle \underset{\mathrm{\Omega }}{\int }{(\nabla N_u)}^{\mathrm{T}}[D](\nabla N_u)\mathrm{d}\mathrm{\Omega }},[K_p]={\displaystyle \underset{\mathrm{\Omega }}{\int }{(\nabla N_p)}^{\mathrm{T}}{\displaystyle \frac{[k]}{{\rho }_{\mathrm{f}}g}}(\nabla N_p)\mathrm{d}\mathrm{\Omega }},[C]={\displaystyle \underset{\mathrm{\Omega }}{\int }{(\nabla N_u)}^{\mathrm{T}}\{m\}[N_p]\mathrm{d}\mathrm{\Omega }},

[C_p]={\displaystyle \underset{\mathrm{\Omega }}{\int }{(N_p)}^{\mathrm{T}}{\displaystyle \frac{\mathrm{1}}{Q}}(N_p)\mathrm{d}\mathrm{\Omega }},[M_s]={\displaystyle \underset{\mathrm{\Omega }}{\int }{(N_u)}^{\mathrm{T}}\rho (N_u)\mathrm{d}\mathrm{\Omega }},

[M_{sp}]={\displaystyle \underset{\mathrm{\Omega }}{\int }{(\nabla N_p)}^{\mathrm{T}}{\displaystyle \frac{[k]}{g}}(N_u)\mathrm{d}\mathrm{\Omega }},[f_{\mathrm{1}}]={\displaystyle \underset{\mathrm{\Gamma }}{\int }{(N_u)}^{\mathrm{T}}\{\sigma \}\mathrm{d}\mathrm{\Gamma }}+{\displaystyle \underset{\mathrm{\Omega }}{\int }{(N_u)}^{\mathrm{T}}\rho \{b\}\mathrm{d}\mathrm{\Omega }},

[f_{\mathrm{2}}]={\displaystyle \underset{\mathrm{\Omega }}{\int }{(\nabla N_p)}^{\mathrm{T}}{\displaystyle \frac{[k]}{{\rho }_{\mathrm{f}}g}}({\rho }_{\mathrm{f}}b)\mathrm{d}\mathrm{\Omega }}+{\displaystyle \underset{\mathrm{\Gamma }}{\int }{(N_p)}^{\mathrm{T}}\left({\displaystyle \frac{k}{{\rho }_{\mathrm{f}}g}}\right)(\nabla p)\mathrm{d}\mathrm{\Gamma }}\mathrm{.}

These three formulations are implemented in an all-purpose FE code and the FE formulations are used to analyze the response of a layered seabed in the free field condition. A harmonic wave load is applied on the seabed surface and the bottom is modeled impermeable and fixed allowing horizontal fluid displacement as can be seen in Fig. 1. First the numerical code is verified through analytical solutions for a single layer using all the formulations as presented in Fig. 2. The analytical solution is briefly discussed below.

Analytical solution

An analytical solution to the dynamics of multi-layer soil under waves is also developed and briefly discussed here. Equation (10) can be re-written as,

(20)Q{\left(u_{i,i}+{\overline{w}}_{i,i}\right)}_i=-p_{,i}.

If we substitute Eq. (14) into Eq. (8b) and use the effective stress relation of Eq. (1), neglecting the body forces yields the two final equations for layers numbered with subscript (l = 1,2…N) as,

(21a){(\sigma {\prime }_{ij,j})}_l+Q_l{\left[{\left(u_{i,i}+{\overline{w}}_{i,i}\right)}_i\right]}_l={\rho }_l{\left({\ddot{u}}_i\right)}_l+{({\rho }_{\mathrm{f}})}_l{\left({\ddot{\overline{w}}}_i\right)}_l,

(21b)Q_l{\left[{\left(u_{i,i}+{\overline{w}}_{i,i}\right)}_i\right]}_l={({\rho }_{\mathrm{f}})}_l{\left({\ddot{u}}_i\right)}_l+{\left({\displaystyle \frac{{\rho }_{\mathrm{f}}}{n}}\right)}_l{\left({\ddot{\overline{w}}}_i\right)}_l+{\left({\displaystyle \frac{{\rho }_{\mathrm{f}}{g}_i}{k_i}}\right)}_l{\left({\dot{\overline{w}}}_i\right)}_l.

Considering the elastic stress-strain relationship using the Lame’s parameters \lambda and G, we have,

(22){({\sigma }_{ij}\prime )}_l={\lambda }_l{({ϵ}_{kk})}_l{\delta }_{ij}+\mathrm{2}{G}_l{({ϵ}_{ij})}_l,

where {\mathbf{ϵ}}_{\mathbf{kk}} is volumetric strain. Writing the scalar form of Eq. (21) incorporating Eq. (22) will require for each layer to have each of these physical parameters {(l,k_x,k_z,n,\lambda ,G,K,K_{\mathrm{f}},\rho ,{\rho }_{\mathrm{f}})}_l. Considering all the response variables to be of the form, R(x,z,t)=R(z)e^{i(kx-\omega t)}, we can obtain a harmonic complex form of the equations where R(z) represents the amplitudes of the response variables which are the solid and fluid displacements, U_x,U_z,{\overline{W}}_x,{\overline{W}}_z, pore pressure, P and the stresses, {{\sigma }^{\prime }}_{xx}, {{\sigma }^{\prime }}_{zz}, {\sigma }_{xz}. However, now that there is N number of layers, a modified normalized depth ratio, Z={\displaystyle \frac{z}{lh}} is introduced into the linear set (with “l” being the thickness ratio of each layer and h the total depth) which is now correspondingly multiplied by l or l² terms. It is necessary to define a set of non-dimensional parameters as was also stated in Ref. [14,15]. So, \kappa ={\displaystyle \frac{Q}{K+Q}}, {\kappa }_{\mathrm{1}}={\displaystyle \frac{\lambda }{K+Q}}, {\kappa }_{\mathrm{2}}={\displaystyle \frac{G}{K+Q}} and \beta ={\displaystyle \frac{{\rho }_{\mathrm{f}}}{\rho }} stay the same but m=\psi lh is now modified to account for the effect of layering on the spatial variation of loading with \psi ={\displaystyle \raisebox{1ex}{$\mathrm{2}\mathrm{\pi }$}\!\left/ \!\raisebox{-1ex}{$L$}\right.} being the wave number and L the wavelength. The main parameters; {\mathrm{\Pi }}_{\mathrm{1}x}={\displaystyle \frac{k_x{V}_{\mathrm{c}}^{\mathrm{2}}}{g\beta \omega l^{\mathrm{2}}{h}^{\mathrm{2}}}} as the ratio of time for pore fluid flow in x-direction to time for compression wave to travel (k_x is the coefficient of hydraulic conductivity in the x-direction); {\mathrm{\Pi }}_{\mathrm{1}z}={\displaystyle \frac{k_z{V}_{\mathrm{c}}^{\mathrm{2}}}{g\beta \omega l^{\mathrm{2}}{h}^{\mathrm{2}}}} to account for the flow in z-direction with k_z being the coefficient of hydraulic conductivity in the z-direction and {\mathrm{\Pi }}_{\mathrm{2}}={\left({\displaystyle \frac{\omega lh}{V_{\mathrm{c}}}}\right)}^{\mathrm{2}} as the ratio of the rate of dynamic loading to the speed of compression wave, V_{\mathrm{c}}, are now also modified due to layering. In addition, it is now necessary to account for the ratios of the non-dimensional parameters of each layers as, R_{{\mathrm{\Pi }}_{\mathrm{1}x}}={\displaystyle \frac{{({\mathrm{\Pi }}_{\mathrm{1}x})}_{l-\mathrm{1}}}{{({\mathrm{\Pi }}_{\mathrm{1}x})}_l}}, R_{{\mathrm{\Pi }}_{\mathrm{1}\mathrm{z}}}={\displaystyle \frac{{({\mathrm{\Pi }}_{\mathrm{1}\mathrm{z}})}_{l-\mathrm{1}}}{{({\mathrm{\Pi }}_{\mathrm{1}\mathrm{z}})}_l}}, R_{{\mathrm{\Pi }}_{\mathrm{2}}}={\displaystyle \frac{{({\mathrm{\Pi }}_{\mathrm{2}})}_{l-\mathrm{1}}}{{({\mathrm{\Pi }}_{\mathrm{2}})}_l}}, R_m={\displaystyle \frac{m_{l-\mathrm{1}}}{m_l}}. It should be noted here that these ratios are used mainly for two distinct layers in a few layer system. The governing equations, in terms of these non-dimensional parameters are written in a matrix form (provided in the Appendix), where we can neglect the wave number, k, and take m = 0 for a single layer, it becomes the 1-D solution provided that {\mathrm{\Pi }}_{\mathrm{1}x}=\mathrm{0} and {\mathrm{\Pi }}_{\mathrm{1}z}={\mathrm{\Pi }}_{\mathrm{1}}. Therefore, for each layer there is {(V_{\mathrm{c}},\kappa ,{\kappa }_{\mathrm{1}},{\kappa }_{\mathrm{2}},\beta ,{\mathrm{\Pi }}_{\mathrm{1}x},{\mathrm{\Pi }}_{\mathrm{1}z},{\mathrm{\Pi }}_{\mathrm{2}})}_l. The boundary conditions are: at z = - h, {(U_x)}_N={(U_z)}_N={(W_z)}_N=\mathrm{0} and at z = 0, vertical effective stress and shear stress of the first layer (l = 1) vanish, {({{\sigma }^{\prime }}_{zz})}_{\mathrm{1}}={({\sigma }_{xz})}_{\mathrm{1}}=\mathrm{0} while pore pressure takes the form p_{\mathrm{1}}=q_{\mathrm{0}}{e}^{i(\psi x-\omega t)} with q_{\mathrm{0}}={\displaystyle \frac{{\gamma }_{w}H}{\mathrm{2}\cosh (kd)}} where H is the wave height. Using these boundary conditions, a 6 × 6 matrix with entries {({\xi }_{ij})}_l (i, j = 1,2,...,6) is obtained in terms of eigenvalues {({\eta }_i)}_l that are calculated from the sixth order characteristic equations of layers in the form,

(23){({\alpha }_{\mathrm{1}}D{D}^{\mathrm{6}}+{\alpha }_{\mathrm{2}}D{D}^{\mathrm{4}}+{\alpha }_{\mathrm{3}}D{D}^{\mathrm{2}}+{\alpha }_{\mathrm{4}})}_l=\mathrm{0},

where α are the coefficients defined as a function of non-dimensional parameters. Solutions of these characteristic equations lead to displacements evaluated first and being valid for -l_{j-\mathrm{1}}h\le z\le -l_{j}h and coefficients a_i obtained from {({\xi }_{ij})}_l and {(b_i,c_i,d_i)}_l evaluated through solving the matrix system. Here in place of the displacement vector, the eigenvector with entries {(b_i,c_i,d_i)}_l is used. Interested readers can refer to Ref. [15] for more details. The semi-analytical solution for the dynamic response of layered porous media requires the use of continuity conditions at layer interfaces. For Z=-l_j, {(U_x)}_j={(U_x)}_{j+\mathrm{1}} and similarly for Z=-(l_j+l_{j+\mathrm{1}}), {(U_x)}_{j+\mathrm{1}}={(U_x)}_{j+\mathrm{2}}. We can write the rest of the continuity conditions as; {(U_z)}_j={(U_z)}_{j+\mathrm{1}} , {(f_z)}_j={(f_z)}_{j+\mathrm{1}}, {({\sigma }_{zz})}_j={({\sigma }_{zz})}_{j+\mathrm{1}}, {({\tau }_{xz})}_j={({\tau }_{xz})}_{j+\mathrm{1}}, P_j=P_{j+\mathrm{1}}. These equations result in a linear system of

(24){[\xi ]}_{MxM}{\{a\}}_{Mx\mathrm{1}}={\{F\}}_{Mx\mathrm{1}},

where M= 6N, [\xi ] is the matrix of equations obtained from both boundary and continuity equations (or {({\xi }_{ij})}_l equations) that can be referred to as layer stiffness matrix, \{a\} is the vector of coefficients a_i and \{F\} is the forcing vector of the form, \{\mathrm{0},\mathrm{0},F,\cdots ,\mathrm{0}\} with F being the only nonzero term coming from the pressure boundary condition at the surface written as F={\displaystyle \frac{q_{\mathrm{0}}}{Q_{\mathrm{1}}}} with Q_{\mathrm{1}}={\left({\displaystyle \frac{K_f}{n}}\right)}_{\mathrm{1}} for the first layer. The linear system of Eq. (24) along with the continuity conditions are presented in the Appendix. The equation system of Eq. (24) should separately be written for all formulations. Then the wave-induced response variables for a particular layer, l, in terms of shear stress and pore water pressure is respectively calculated as,

(25){({\tau }_{xz})}_l=G_l{\left[{\displaystyle \sum _{j=\mathrm{1}}^{\mathrm{6}}{\left({\displaystyle \frac{{\eta }_j}{lh}}+i\psi b_j\right)}_l{a}_j{e}^{{\eta }_j\frac{z}{lh}}}\right]}_l{e}^{i(\psi x-\omega t)},

(26)p_l=-{\left({\displaystyle \frac{K_f}{n}}\right)}_l{\left[{\displaystyle \sum _{j=\mathrm{1}}^{\mathrm{6}}{\left(i\psi (\mathrm{1}+c_j)+{\displaystyle \frac{{\eta }_j}{lh}}(b_j+d_j)\right)}_l{a}_j{e}^{{\eta }_j\frac{z}{lh}}}\right]}_l{e}^{i(\psi x-\omega t)}.

Dynamic response of multi-layer systems

The dynamic response of a multi-layer soil is presented in terms of absolute values of pore pressure and shear stress variations in the medium. Here non-dimensional parameters for each layer and their ratios (R_{{\mathrm{\Pi }}_{\mathrm{1}x}}, R_{{\mathrm{\Pi }}_{\mathrm{1}z}}, R_{{\mathrm{\Pi }}_{\mathrm{2}}}, R_m) play a key role in defining the continuous behavior of multiple layers under dynamic loading. It is determined that in case there are only a few layers with distinct physical properties, depending on their respective thicknesses, the first distinct consecutive layers closest to the surface are important in identifying the necessary formulation governing the response. If there are more than two layers with thicknesses being relatively much larger or that there is larger difference between some of the typical physical properties (i.e., hydraulic conductivity, degree of saturation), then a weighted average of those properties using the relation below will be sufficient,

(27)A_{\mathrm{e}\mathrm{q}}={\displaystyle \frac{{\displaystyle \sum _{l=\mathrm{1}}^N{h}_l{A}_l}}{{\displaystyle \sum _{l=\mathrm{1}}^N{h}_l}}}.

Here A_{\mathrm{e}\mathrm{q}} is the equivalent response parameter to be weighted. Figures 3-10 present the dynamic response of a layered system for the three formulations in terms of pore pressure and shear stress variations. The effect of wave period (T) in terms of Π₂, the effect of hydraulic conductivity ratio in terms of R_{{\mathrm{\Pi }}_{\mathrm{1}}} and the effect of layer thickness ratio in terms of R_m on the response are illustrated. While there is no significant discrepancy between the PD and FD solutions in pore pressure distributions, Fig. 5 yields some difference due to the boundary layer effect. Figures 3 and 4 present the variation of pore pressure with depth for three R_{{\mathrm{\Pi }}_{\mathrm{1}}}, two m values and two wave periods. As obvious, the response in terms of PD and QS formulations differs for m and T values. Difference between the formulations increases as the m value increases irrespective of Π₂. It can be seen from some of the figures that formulations converging to the same response reduce the necessity of including the inertial terms in the system. In this case even the QS formulation is sufficient to capture the dynamic behavior. However, as the discrepancy increases, beyond maximum absolute value of 3%, inertial terms associated with the motion of the solid part and the relative motion of the fluid part become important.

Figures 6-8 present the normalized shear stress response within the depth. There is significant boundary layer effect apparent in the variation of shear stresses when the inertial terms are included (Fig. 7) for low wave periods. Also dynamic formulations become significant for higher values of m (Fig. 8).

The effect of layering is presented in Figs. 9 and 10 where varying values of respective ratios of non-dimensional parameters both illustrate the discrepancy in the formulations and their effect on the pore pressure and shear stress response.

Closure

In this study, solutions to the dynamic response of two-dimensional saturated and layered soils are derived. Finite element numerical formulations to the governing poroelasticity equations are used to solve the dynamic response of soil layers and the results are first validated with corresponding analytical solutions which are also developed by the author. For this, the three formulations of quasi-static, partially dynamic and fully dynamic are used. The dynamic response is evaluated in terms of absolute value distributions of pore pressure and shear stresses within the soil layers. Results indicate that the non-dimensional parameter ratios for consecutive layers using physical quantities play a key role in obtaining the response. Their variations determine what formulation to use to capture the behavior in a layered system. The solutions developed in this study can be of use to researchers and practitioners in obtaining the dynamic response of soils in a given layered system.

References

[1] Rahman M S, El-Zahaby K, Booker J R. A semi-analytical method for the wave induced seabed response. International Journal for Numerical and Analytical Methods in Geomechanics, 1994, 18(4): 213-236

[2] Ai Z Y, Zeng W Z. Analytical layer-element method for non-axisymmetric consolidation of multilayered soils. International Journal for Numerical and Analytical Methods in Geomechanics, 2012, 36(5): 533-545

[3] Ulker M B C. Pore pressure, stress distributions and instantaneous liquefaction of two-layer soil under waves. Journal of Waterway, Port, Coastal, and Ocean Engineering, 2012, 138(6): 435-450

[4] Jeng D S, Rahman M S, Lee T L. Effects of inertia forces on wave-induced seabed response. International Journal of Offshore and Polar Engineering, 1999, 9: 307-313

[5] Tsai C P, Lee T L, Hsu J R C. Effects of wave non-linearity on the standing wave-induced seabed response. International Journal for Numerical and Analytical Methods in Geomechanics, 2000, 24(11): 869-892

[6] Jeng D S, Cha D H. Effects of dynamic soil behavior and wave non-linearity on the wave induced pore pressure and effective stresses in porous seabed. Ocean Engineering, 2003, 30(16): 2065-2089

[7] Booker J R, Small J C. A method of computing the consolidation behavior of layered soils using direct numerical inversion of Laplace transforms. International Journal for Numerical and Analytical Methods in Geomechanics, 1987, 11(4): 363-380

[8] Hsu J R C, Jeng D S, Lee C P. Oscillatory soil response and liquefaction in an unsaturated layered seabed. International Journal for Numerical and Analytical Methods in Geomechanics, 1995, 19(12): 825-849

[9] Mesgouez A, Lefeuve-Mesgouez G. Transient solution for multilayered poroviscoelastic media obtained by an exact stiffness matrix formulation. International Journal for Numerical and Analytical Methods in Geomechanics, 2009, 33(18): 1911-1931

[10] Biot M A. General theory of three dimensional consolidation. Journal of Applied Physics, 1941, 12(2): 155-164

[11] Biot M A. Theory of elasticity and consolidation for a porous anisotropic solid. Journal of Applied Physics, 1955, 26(2): 182-185

[12] Biot M A. Mechanics of deformation and acoustic propagation in porous media. Journal of Applied Physics, 1962, 33(4): 1482-1498

[13] Okusa S. Wave induced stresses in unsaturated submarine sediments. Geotechnique, 1985, 35(4): 517-532

[14] Zienkiewicz O C, Chang C T, Bettess P. Drained, undrained, consolidating and dynamic behavior assumptions in soils. Geotechnique, 1980, 30(4): 385-395

[15] Ulker M B C, Rahman M S. Response of saturated porous media: Different formulations and their validity. International Journal for Numerical and Analytical Methods in Geomechanics, 2009, 33(5): 633-664

RIGHTS & PERMISSIONS

Higher Education Press and Springer-Verlag Berlin Heidelberg