Uncoupled state space solution to layered poroelastic medium with anisotropic permeability and compressible pore fluid

Abstract

This paper presents an uncoupled state space solution to three-dimensional consolidation of layered poroelastic medium with anisotropic permeability and compressible pore fluid. Starting from the basic equations of poroelastic medium, and introducing intermediate variables, the state space equation usually comprising eight coupled state vectors is uncoupled into two sets of equations of six and two state vectors in the Laplace-Fourier transform domain. Combined with the continuity conditions between adjacent layers and boundary conditions, the uncoupled state space solution of a layered poroelastic medium is obtained by using the transfer matrix method. Numerical results show that the anisotropy of permeability and the compressibility of pore fluid have remarkable influence on the consolidation behavior of poroelastic medium.

Keywords

uncoupled state space solution / layered poroelastic medium / three-dimensional consolidation / anisotropic permeability / compressible pore fluid

Introduction

Saturated soils and rocks are usually regarded as poroelastic medium, and the general theory of three-dimensional consolidation of poroelastic medium, which considers coupling between solid and fluid, was first introduced by Biot [1]. Many studies [2-12] have been conducted on Biot’s consolidation problem, however, in almost all of these studies, permeability was assumed to be isotropic. In the actual engineering situations, the soils and rocks are formed through a sedimentation process which produces the horizontal stratification. Consequently, the horizontal permeability often differs from that of the vertical [13-15]. Some experimental results have shown that horizontal permeability of the soil may be an order of magnitude, or even more greater than vertical permeability. On the other hand, the fluid phase of saturated soils and rocks is usually in the form of liquid-gas mixture, in which the gas is regarded as minute bubbles occupying in the liquid when the degree of saturation is high (approximately more than 70%), thereby the fluid phase behaves as homogeneous compressible pore fluid [16-27]. These studies further pointed out that the compressibility of pore fluid evidently influences the consolidation process and pore pressure distribution. Therefore, it is necessary to consider the anisotropy of permeability and the compressibility of the fluid phase of saturated soils and rocks in the consolidation calculation. Up to now, few studies considering the anisotropy of permeability and the compressibility of the fluid phase at the same time are conducted on the Biot’s consolidation of poroelastic medium. Booker and Carter [16] analyzed the rate of consolidation with a point sink in an elastic half space by considering the anisotropy of permeability and the compressibility of pore fluid. Chen [17] used the state vector method to analyze the axisymmetric consolidation of layered half space with anisotropic permeability and compressible pore fluid in the cylindrical coordinate system, which was suitable for circular loading. Recently, Ai and Wu [18] studied the plane-strain consolidation of a layered soil, and the anisotropy of permeability and the compressibility of pore fluid were taken into account. An examination of publications indicated that uncoupled state space solution to three-dimensional layered poroelastic medium with anisotropic permeability and compressible pore fluid has not been obtained. Therefore, the main purposes of this paper can be summarized as follows: (1) to present uncoupled state space solution to three-dimensional consolidation of layered poroelastic medium with anisotropic permeability and compressible pore fluid subjected to an arbitrary rectangle surface load; (2) to investigate the influence of the anisotropy of permeability and the compressibility of pore fluid on the consolidation behavior of poroelastic medium.

Uncoupled state space equations and their solutions

The governing equations of three-dimensional Biot’s consolidation of a saturated poroelastic medium with anisotropic permeability and compressible pore fluid can be expressed as follows [28,29]:

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

(1b){\sigma }_{ij,j}=\mathrm{0},

(1c){\sigma }_{ij}=\mathrm{2}G{ϵ}_{ij}+\lambda e{\delta }_{ij}-\sigma {\delta }_{ij},

(2){\displaystyle \frac{\partial e}{\partial t}}+n\beta {\displaystyle \frac{\partial \sigma }{\partial t}}={\displaystyle \frac{\mathrm{1}}{{\gamma }_{\mathrm{w}}}}\left(k_x{\displaystyle \frac{{\partial }^{\mathrm{2}}}{\partial x^{\mathrm{2}}}}+k_y{\displaystyle \frac{{\partial }^{\mathrm{2}}}{\partial y^{\mathrm{2}}}}+k_z{\displaystyle \frac{{\partial }^{\mathrm{2}}}{\partial z^{\mathrm{2}}}}\right)\sigma ,

where {ϵ}_{ij} is Cauchy’s strain tensor, u_i is the displacement in the x_i direction; e=u_{i,i}, is the volumetric strain; {\sigma }_{ij} is the total stress tensor, \sigma is the excess pore pressure (positive under compression), G and v are the shear modulus and Poisson’s ratio of poroelastic medium, respectively; \lambda is the Lame constant defined as \lambda ={\displaystyle \frac{\mathrm{2}Gv}{\mathrm{1}-\mathrm{2}v}}; {\sigma }_{ij} is the kronecker delta, n is porosity, \beta is the compressibility coefficient of pore fluid, {\gamma }_w is the unit weight of water, k_x, k_y and k_z are the coefficient of the permeability in the x, y and z coordinate directions, respectively; let k_x=k_y=k_h, k_z=k_v, and anisotropic ratio parameter \gamma =k_h/k_v.

According to Darcy’s law, the flux in the z direction can be defined as:

(3)Q={\displaystyle \frac{k_v}{{\gamma }_{\mathrm{w}}}}{\displaystyle \frac{\partial \sigma }{\partial z}}.

In order to solve the partial differential Eqs. (1)-(2), we take a Laplace transform with respect to the parameter t and a double Fourier transform with respect to the parameters x and y in Eqs. (1)-(2).

The Laplace transform and its inversion [30] are defined as follows:

\tilde{f}(x,y,z,s)={\displaystyle {\int }_{\mathrm{0}}^{\infty }f(x,y,z,t){\mathrm{e}}^{-st}\mathrm{d}t},f(x,y,z,t)={\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathrm{\pi }i}}{\displaystyle {\int }_{\gamma -i\infty }^{\gamma +i\infty }\tilde{f}(x,y,z,s){\mathrm{e}}^{st}\mathrm{d}s},

where s denotes the Laplace transform parameter. The overbar ‘~’ is used to denote the Laplace transform of a given variable.

The double Fourier transform and its inversion [31] are defined as:

\overline{f}({\xi }_x,{\xi }_y,z,s)={\displaystyle \frac{\mathrm{1}}{\mathrm{4}{\mathrm{\pi }}^{\mathrm{2}}}}{\displaystyle {\int }_{-\infty }^{\infty }{\displaystyle {\int }_{-\infty }^{\infty }\tilde{f}(x,y,z,s)}}{\mathrm{e}}^{-i(x{\xi }_x+y{\xi }_y)}\mathrm{d}x\mathrm{d}y,

\tilde{f}(x,y,z,s)={\displaystyle {\int }_{-\infty }^{\infty }{\displaystyle {\int }_{-\infty }^{\infty }\overline{f}({\xi }_x,{\xi }_y,z,s)}{\mathrm{e}}^{i(x{\xi }_x+y{\xi }_y)}\mathrm{d}{\xi }_x\mathrm{d}{\xi }_y}.

where {\xi }_x and {\xi }_y denote the double Fourier transform parameters, respectively. The overbar ‘-’ is hereafter used to denote the Laplace-Fourier transform of a given variable.

The state space equation of three-dimensional consolidation usually comprises eight coupled state vectors [11,12]. In order to simplify derivation and calculation, intermediate variables are introduced as follows:

(4a)U={\displaystyle \frac{\partial u_x}{\partial x}}+{\displaystyle \frac{\partial u_y}{\partial y}},X={\displaystyle \frac{\partial {\sigma }_{xz}}{\partial x}}+{\displaystyle \frac{\partial {\sigma }_{yz}}{\partial y}},

(4b)V={\displaystyle \frac{\partial u_x}{\partial y}}-{\displaystyle \frac{\partial u_y}{\partial x}},Y={\displaystyle \frac{\partial {\sigma }_{xz}}{\partial y}}-{\displaystyle \frac{\partial {\sigma }_{yz}}{\partial x}}.

Deriving the variables in Eqs. (4a) and (4b) for z-derivative yields:

(5a){\displaystyle \frac{\partial U}{\partial z}}={\displaystyle \frac{{\partial }^{\mathrm{2}}{u}_x}{\partial x\partial z}}+{\displaystyle \frac{{\partial }^{\mathrm{2}}{u}_y}{\partial y\partial z}},

(5b){\displaystyle \frac{\partial X}{\partial z}}={\displaystyle \frac{{\partial }^{\mathrm{2}}{\sigma }_{xz}}{\partial x\partial z}}+{\displaystyle \frac{{\partial }^{\mathrm{2}}{\sigma }_{yz}}{\partial y\partial z}},

(5c){\displaystyle \frac{\partial V}{\partial z}}={\displaystyle \frac{{\partial }^{\mathrm{2}}{u}_x}{\partial y\partial z}}-{\displaystyle \frac{{\partial }^{\mathrm{2}}{u}_y}{\partial x\partial z}},

(5d){\displaystyle \frac{\partial Y}{\partial z}}={\displaystyle \frac{{\partial }^{\mathrm{2}}{\sigma }_{xz}}{\partial y\partial z}}-{\displaystyle \frac{{\partial }^{\mathrm{2}}{\sigma }_{yz}}{\partial x\partial z}}.

Combining with Eqs. (1a) and (1c), we find:

(6a){\displaystyle \frac{\partial u_x}{\partial z}}={\displaystyle \frac{{\sigma }_{xz}}{G}}-{\displaystyle \frac{\partial u_z}{\partial x}},

(6b){\displaystyle \frac{\partial u_y}{\partial z}}={\displaystyle \frac{{\sigma }_{yz}}{G}}-{\displaystyle \frac{\partial u_z}{\partial y}},

(6c){\displaystyle \frac{\partial u_z}{\partial z}}=\left({\displaystyle \frac{\mathrm{2}G}{M}}-\mathrm{1}\right)\left({\displaystyle \frac{\partial u_x}{\partial x}}+{\displaystyle \frac{\partial u_y}{\partial y}}\right)+{\displaystyle \frac{\mathrm{1}}{M}}(\sigma +{\sigma }_z).

Adopting the Laplace-Fourier transform to Eqs. (5a) and (6), we obtain:

(7){\displaystyle \frac{\mathrm{d}\overline{U}}{\mathrm{d}z}}=i{\xi }_x{\displaystyle \frac{\mathrm{d}{\overline{u}}_x}{\mathrm{d}z}}+i{\xi }_y{\displaystyle \frac{\mathrm{d}{\overline{u}}_y}{\mathrm{d}z}}=i{\xi }_x\left({\displaystyle \frac{\mathrm{1}}{G}}{\overline{\sigma }}_{xz}-i{\xi }_x{\overline{u}}_z\right)+i{\xi }_y\left({\displaystyle \frac{\mathrm{1}}{G}}{\overline{\sigma }}_{yz}-i{\xi }_y{\overline{u}}_z\right)={\xi }^{\mathrm{2}}{\overline{u}}_z+{\displaystyle \frac{\mathrm{1}}{G}}\overline{X},

(8){\displaystyle \frac{\mathrm{d}{\overline{u}}_z}{\mathrm{d}z}}=\left({\displaystyle \frac{\mathrm{2}G}{M}}-\mathrm{1}\right)\left({\xi }_{x}i{\overline{u}}_x+{\xi }_{y}i{\overline{u}}_y\right)+{\displaystyle \frac{\mathrm{1}}{M}}\left({\overline{\sigma }}_z+\overline{\sigma }\right)=\left({\displaystyle \frac{\mathrm{2}G}{M}}-\mathrm{1}\right)\overline{U}+{\displaystyle \frac{\mathrm{1}}{M}}\left({\overline{\sigma }}_z+\overline{\sigma }\right),

where {\xi }_x^{\mathrm{2}}+{\xi }_y^{\mathrm{2}}={\xi }^{\mathrm{2}} and M=\lambda +\mathrm{2}G.

From Eq. (1b), we have:

(9a){\displaystyle \frac{\partial {\sigma }_{xz}}{\partial z}}=-{\displaystyle \frac{\partial {\sigma }_x}{\partial x}}-{\displaystyle \frac{\partial {\sigma }_{xy}}{\partial y}},

(9b){\displaystyle \frac{\partial {\sigma }_{yz}}{\partial z}}=-{\displaystyle \frac{\partial {\sigma }_{yx}}{\partial x}}-{\displaystyle \frac{\partial {\sigma }_y}{\partial y}},

(9c){\displaystyle \frac{\partial {\sigma }_z}{\partial z}}=-{\displaystyle \frac{\partial {\sigma }_{zx}}{\partial x}}-{\displaystyle \frac{\partial {\sigma }_{zy}}{\partial y}}.

Applying the Laplace-Fourier transform to Eqs. (5b) and (9), we have:

(10){\displaystyle \frac{\mathrm{d}\overline{X}}{\mathrm{d}z}}=i{\xi }_x{\displaystyle \frac{\mathrm{d}{\overline{\sigma }}_{xz}}{\mathrm{d}z}}+i{\xi }_y{\displaystyle \frac{\mathrm{d}{\overline{\sigma }}_{yz}}{\mathrm{d}z}}={\xi }_x^{\mathrm{2}}{\overline{\sigma }}_x+{\xi }_y^{\mathrm{2}}{\overline{\sigma }}_y+\mathrm{2}{\xi }_x{\xi }_y{\overline{\sigma }}_{xy},

(11){\displaystyle \frac{\mathrm{d}{\overline{\sigma }}_z}{\mathrm{d}z}}=-{\xi }_{x}i{\overline{\sigma }}_{zx}-{\xi }_{y}i{\overline{\sigma }}_{zy}=-\overline{X},

where the expressions of {\overline{\sigma }}_x, {\overline{\sigma }}_y and {\overline{\sigma }}_{xy} can be derived from Eq. (1c), therefore Eq. (10) can be further re-structured as follows:

(12){\displaystyle \frac{\mathrm{d}\overline{X}}{\mathrm{d}z}}=\mathrm{4}G\left(\mathrm{1}-{\displaystyle \frac{G}{M}}\right){\xi }^{\mathrm{2}}\overline{U}+\left(\mathrm{1}-{\displaystyle \frac{\mathrm{2}G}{M}}\right){\xi }^{\mathrm{2}}{\overline{\sigma }}_z-{\displaystyle \frac{\mathrm{2}G}{M}}{\xi }^{\mathrm{2}}\overline{\sigma }.

Applying the Laplace-transform to Eqs. (2) and (3), and then combining with Eq. (8), we can obtain:

(13){\displaystyle \frac{\mathrm{d}\overline{\sigma }}{\mathrm{d}z}}={\displaystyle \frac{{\gamma }_{\mathrm{w}}}{k_h}}\overline{Q},

(14){\displaystyle \frac{\mathrm{d}\overline{Q}}{\mathrm{d}z}}={\displaystyle \frac{k_z}{{\gamma }_{\mathrm{w}}}}{\displaystyle \frac{{\mathrm{d}}^{\mathrm{2}}\overline{\sigma }}{\mathrm{d}{z}^{\mathrm{2}}}}={\displaystyle \frac{\mathrm{2}Gs}{M}}\overline{U}+{\displaystyle \frac{s}{M}}{\overline{\sigma }}_z+({\displaystyle \frac{s}{M}}+sn\beta +{\displaystyle \frac{k_h}{{\gamma }_{\mathrm{w}}}}{\xi }^{\mathrm{2}})\overline{\sigma }.

Similarly, adopting the Laplace-Fourier transform to Eqs. (4b), (5c) and (6), we obtain:

(15){\displaystyle \frac{\mathrm{d}\overline{V}}{\mathrm{d}z}}=i{\xi }_y{\displaystyle \frac{\mathrm{d}{\overline{u}}_x}{\mathrm{d}z}}-i{\xi }_x{\displaystyle \frac{\mathrm{d}{\overline{u}}_y}{\mathrm{d}z}}=i{\xi }_y\left({\displaystyle \frac{\mathrm{1}}{G}}{\overline{\sigma }}_{xz}-i{\xi }_x{\overline{u}}_z\right)-i{\xi }_x\left({\displaystyle \frac{\mathrm{1}}{G}}{\overline{\sigma }}_{yz}-i{\xi }_y{\overline{u}}_z\right)={\displaystyle \frac{\mathrm{1}}{G}}\overline{Y}.

Applying the Laplace-Fourier transform to Eqs. (4b), (5d) and (9), it results in:

(16){\displaystyle \frac{\mathrm{d}\overline{Y}}{\mathrm{d}z}}=i{\xi }_y{\displaystyle \frac{\mathrm{d}{\overline{\sigma }}_{xz}}{\mathrm{d}z}}-i{\xi }_x{\displaystyle \frac{\mathrm{d}{\overline{\sigma }}_{yz}}{\mathrm{d}z}}={\xi }_x{\xi }_y\left({\overline{\sigma }}_x-{\overline{\sigma }}_y\right)+({\xi }_y^{\mathrm{2}}-{\xi }_x^{\mathrm{2}}){\overline{\sigma }}_{xy}=G{\xi }^{\mathrm{2}}\overline{V},

where the expressions of {\overline{\sigma }}_x, {\overline{\sigma }}_y and {\overline{\sigma }}_{xy} can be derived from Eq. (1c).

Assembling Eqs. (7), (8), (11), (12), (13), (14), (15), and (16), we can obtain two sets of equations of six and two state vectors in the Laplace-Fourier transform domain, which can be expressed in matrix form as follows:

(17a){\displaystyle \frac{\mathrm{d}}{\mathrm{d}z}}\left[{\overline{B}}_{\mathrm{1}}({\xi }_x,{\xi }_y,z,s)\right]=A_{\mathrm{1}}({\xi }_x,{\xi }_y,s){\overline{B}}_{\mathrm{1}}({\xi }_x,{\xi }_y,z,s),

(17b){\displaystyle \frac{\mathrm{d}}{\mathrm{d}z}}\left[{\overline{B}}_{\mathrm{2}}({\xi }_x,{\xi }_y,z,s)\right]=A_{\mathrm{2}}({\xi }_x,{\xi }_y,s){\overline{B}}_{\mathrm{2}}({\xi }_x,{\xi }_y,z,s).

where,

A_{\mathrm{1}}=\left[\begin{array}{llllll}\mathrm{0}& {\xi }^{\mathrm{2}}& {\displaystyle \frac{\mathrm{1}}{G}}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ ({\displaystyle \frac{\mathrm{2}G}{M}}-\mathrm{1})& \mathrm{0}& \mathrm{0}& {\displaystyle \frac{\mathrm{1}}{M}}& {\displaystyle \frac{\mathrm{1}}{M}}& \mathrm{0}\\ (\mathrm{4}G-{\displaystyle \frac{\mathrm{4}{G}^{\mathrm{2}}}{M}}){\xi }^{\mathrm{2}}& \mathrm{0}& \mathrm{0}& (\mathrm{1}-{\displaystyle \frac{\mathrm{2}G}{M}}){\xi }^{\mathrm{2}}& -{\displaystyle \frac{\mathrm{2}G}{M}}{\xi }^{\mathrm{2}}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& -\mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& {\displaystyle \frac{{\gamma }_w}{k_v}}\\ s{\displaystyle \frac{\mathrm{2}G}{M}}& \mathrm{0}& \mathrm{0}& {\displaystyle \frac{s}{M}}& s\delta +{\displaystyle \frac{k_h}{{\gamma }_w}}{\xi }^{\mathrm{2}}& \mathrm{0}\end{array}\right],A_{\mathrm{2}}=\left[\begin{array}{ll}\mathrm{0}& {\displaystyle \frac{\mathrm{1}}{G}}\\ G{\xi }^{\mathrm{2}}& \mathrm{0}\end{array}\right],

{\overline{\mathit{B}}}_{\mathrm{1}}({\xi }_x,{\xi }_y,z,s)={\left[\overline{U},{\overline{u}}_z,\overline{X},{\overline{\sigma }}_z,\overline{\sigma },\overline{Q}\right]}^{\mathrm{T}},{\overline{\mathit{B}}}_{\mathrm{2}}({\xi }_x,{\xi }_y,z,s)={\left[\overline{V},\overline{Y}\right]}^{\mathrm{T}},\delta ={\displaystyle \frac{\mathrm{1}}{M}}+n\beta .

Let \alpha ={\displaystyle \frac{\mathrm{1}+n\beta (\lambda +\mathrm{2}G)}{\mathrm{1}+n\beta (\lambda +G)}}, which is a factor related to the compressibility of pore fluid (obviously ), then \delta ={\displaystyle \frac{\alpha G}{M[(\mathrm{1}-\alpha )M+\alpha G]}}.

Equations (17a) and (17b) are ordinary differential equations in the form of matrix, and their solutions are:

(18a){\overline{\mathit{B}}}_1\left({\xi }_x,{\xi }_y,z,s\right)=\mathit{\Phi }\left({\xi }_x,{\xi }_y,z,s\right){\overline{\mathit{B}}}_1\left({\xi }_x,{\xi }_y,0,s\right),

(18b){\overline{\mathit{B}}}_{\mathrm{2}}({\xi }_x,{\xi }_y,z,s)=\mathit{\Psi }({\xi }_x,{\xi }_y,z,s){\overline{\mathit{B}}}_{\mathrm{2}}({\xi }_x,{\xi }_y,\mathrm{0},s),

where \mathit{\Phi }({\xi }_x,{\xi }_y,z,s)={\mathrm{e}}^{z{\mathit{A}}_{\mathrm{1}}({\xi }_x,{\xi }_y,s)}, \mathit{\Psi }({\xi }_x,{\xi }_y,z,s)={\mathrm{e}}^{z{\mathit{A}}_{\mathrm{2}}({\xi }_x,{\xi }_y,s)}. Based on the theorem of Cayley-Hamilton [32], we have derived \mathit{\Phi }({\xi }_x,{\xi }_y,z,s) and \mathit{\Psi }({\xi }_x,{\xi }_y,z,s), whose elements are listed in the Appendix.

Solutions of layered poroelastic medium

As shown in Fig. 1, an arbitrary load is applied to the surface of a saturated poroelastic medium with n layers and the components of the load in the x, y and z directions are p\left(x,y,0,t\right),r\left(x,y,0,t\right) and f(x,y,\mathrm{0},t), respectively. Let the thickness of ith layer be \mathrm{\Delta }{H}_i=H_i-H_{i-\mathrm{1}}, where H_i and H_{i-\mathrm{1}} are the depths from the ground surface to the bottom and top of the ith layer, respectively.

It is assumed that the surface of the medium is completely permeable, then:

(19a){\overline{\sigma }}_z({\xi }_x,{\xi }_y,\mathrm{0},s)=-\overline{f}\overline{\sigma }({\xi }_x,{\xi }_y,\mathrm{0},s)=\mathrm{0},

(19b)\overline{X}({\xi }_x,{\xi }_y,\mathrm{0},s)=i{\xi }_x\overline{p}+i{\xi }_y\overline{r}\overline{Y}({\xi }_x,{\xi }_y,\mathrm{0},s)=i{\xi }_y\overline{p}-i{\xi }_x\overline{r}.

The bottom of poroelastic medium is supposed to be fixed and completely impermeable, then:

(20)\overline{\mathit{U}}({\xi }_x,{\xi }_y,H_n,s)=\overline{\mathit{V}}({\xi }_x,{\xi }_y,H_n,s)={\overline{\mathit{u}}}_z({\xi }_x,{\xi }_y,H_n,s)=\overline{\mathit{Q}}({\xi }_x,{\xi }_y,H_n,s)=\mathrm{0.}

.

According to the compatibility of displacement and the continuity of stress at the interface between any two adjacent layers, the following continuity conditions can be obtained:

(21a){\overline{\mathit{B}}}_{\mathrm{1}}({\xi }_x,{\xi }_y,H_i^{-},s)={\overline{\mathit{B}}}_{\mathrm{1}}({\xi }_x,{\xi }_y,H_i^{+},s),

(21b){\overline{\mathit{B}}}_{\mathrm{2}}({\xi }_x,{\xi }_y,H_i^{-},s)={\overline{\mathit{B}}}_{\mathrm{2}}({\xi }_x,{\xi }_y,H_i^{+},s),

where {\overline{\mathit{B}}}_{\mathrm{1}}({\xi }_x,{\xi }_y,H_i^{+},s) holds the value of variables at the surface of the (i+1)th layer when z=H_i and {\overline{\mathit{B}}}_{\mathrm{1}}({\xi }_x,{\xi }_y,H_i^{-},s) holds those at the bottom of the ith layer when z=H_i, so is {\overline{\mathit{B}}}_{\mathrm{2}}.

Applying equations (18) to all the poroelastic layers, the following equations can be obtained:

(22a){\overline{\mathit{B}}}_{\mathrm{1}}({\xi }_x,{\xi }_y,H_n,s)={\mathit{\prod }}_{\mathrm{1}}·{\overline{\mathit{B}}}_{\mathrm{1}}({\xi }_x,{\xi }_y,\mathrm{0},s),

(22b){\overline{\mathit{B}}}_{\mathrm{2}}({\xi }_x,{\xi }_y,H_n,s)={\mathit{\prod }}_{\mathrm{2}}·{\overline{\mathit{B}}}_{\mathrm{2}}({\xi }_x,{\xi }_y,\mathrm{0},s),

where {\mathit{\Pi }}_{\mathrm{1}}=\mathit{\Phi }({\xi }_x,{\xi }_y,\mathrm{\Delta }{H}_n,s)\mathit{\Phi }({\xi }_x,{\xi }_y,\mathrm{\Delta }{H}_{n-\mathrm{1}},s)\cdots \mathit{\Phi }({\xi }_x,{\xi }_y,\mathrm{\Delta }{H}_{\mathrm{1}},s) and {\mathit{\Pi }}_{\mathrm{2}}=\mathit{\Psi }({\xi }_x,{\xi }_y,\mathrm{\Delta }{H}_n,s)\mathit{\Psi }({\xi }_x,{\xi }_y,\mathrm{\Delta }{H}_{n-\mathrm{1}},s)\cdots \mathit{\Psi }({\xi }_x,{\xi }_y,\mathrm{\Delta }{H}_{\mathrm{1}},s).

For the known boundary condition, {\overline{\mathit{B}}}_{\mathrm{1}}({\xi }_x,{\xi }_y,\mathrm{0},s), {\overline{\mathit{B}}}_{\mathrm{1}}({\xi }_x,{\xi }_y,H_n,s), {\overline{\mathit{B}}}_{\mathrm{2}}({\xi }_x,{\xi }_y,\mathrm{0},s) and {\overline{\mathit{B}}}_{\mathrm{2}}({\xi }_x,{\xi }_y,H_n,s) can be derived from equations (22). Therefore, for a given depth z in the ith layer, the variables in the transformed domain can be obtained through the following relationships:

(23a){\overline{\mathit{B}}}_{\mathrm{1}}({\xi }_x,{\xi }_y,z,s)={\mathit{\Pi }}_{z\mathrm{1}}·{\overline{\mathit{B}}}_{\mathrm{1}}({\xi }_x,{\xi }_y,H_n,s),

(23b){\overline{\mathit{B}}}_{\mathrm{2}}({\xi }_x,{\xi }_y,z,s)={\mathit{\Pi }}_{z\mathrm{2}}·{\overline{\mathit{B}}}_{\mathrm{2}}({\xi }_x,{\xi }_y,H_n,s),

where {\mathit{\Pi }}_{z\mathrm{1}}=\mathit{\Phi }({\xi }_x,{\xi }_y,z-H_i,s)\mathit{\Phi }({\xi }_x,{\xi }_y,-\mathrm{\Delta }{H}_{i+\mathrm{1}},s)\cdots \mathit{\Phi }({\xi }_x,{\xi }_y,-\mathrm{\Delta }{H}_n,s) and {\mathit{\Pi }}_{z\mathrm{2}}=\mathit{\Psi }({\xi }_x,{\xi }_y,z-H_i,s)\mathit{\Psi }({\xi }_x,{\xi }_y,-\mathrm{\Delta }{H}_{i+\mathrm{1}},s)\cdots \mathit{\Psi }({\xi }_x,{\xi }_y,-\mathrm{\Delta }{H}_n,s).

Equations (23) are the solutions for the consolidation problem of the layered poroelastic medium in the transformed domain. By adopting the inversion of the Laplace-Fourier transform to {\overline{\mathit{B}}}_{\mathrm{1}}({\xi }_x,{\xi }_y,z,s) and {\overline{\mathit{B}}}_{\mathrm{2}}({\xi }_x,{\xi }_y,z,s), the real solutions for the stresses, displacements, excess pore pressure and flux in the layered poroelastic medium with anisotropic permeability and compressible constituents can be obtained.

Numerical results

In this paper, a numerical solution for the Laplace inverse transform has been presented by adopting Talbot’s method [30], i.e.:

(24)f(t)={\displaystyle \frac{\delta e^{\eta t}}{n}}{{\displaystyle \sum _{k=\mathrm{0}}^{n-\mathrm{1}}[e^{\phi \delta t}\{(\upsilon R-\beta I)\cos (\upsilon \theta \delta t)-(\upsilon I+\beta R)\sin (\upsilon \theta \delta t)\}]}}_{\theta ={\theta }_k},

where {\theta }_k=k\mathrm{\pi }/n,k=\mathrm{0},\mathrm{1},\mathrm{2},\cdots ,n-\mathrm{1} and n is determined by the accuracy and efficiency; \phi =\theta \cot \theta and \beta =\theta +\phi (\phi -\mathrm{1})/\theta are the functions of \theta; the values of R and I are defined by the equation of \tilde{f}(\delta s_v+\eta )=R+iI, where s_v=\phi +i\upsilon \theta. Consequently, the realization of inverse Laplace transform is based on the three parameters of \delta ,\eta ,\upsilon.

Let the singularities of \tilde{f}(s) be at s_j=p_j+i{q}_j, we note:

(25){\eta }_{\mathrm{0}}=\max \left(\mathrm{0},\widehat{p}\right),

where \widehat{p}=\underset{j}{\max }{p}_j, then let:

(26)s_j^{\prime }=s_j-{\eta }_{\mathrm{0}}=p_j^{\prime }+i{q}_j,{\eta }^{\prime }=\eta -{\eta }_{\mathrm{0}},

where p_j^{\prime }=p_j-{\eta }_{\mathrm{0}}\le \mathrm{0}.

the most significant or dominant singularity s_{\mathrm{d}} is defined as:

(27){\displaystyle \frac{q_{\mathrm{d}}}{{\theta }_{\mathrm{d}}}}=\underset{q_j>\mathrm{0}}{\max }{\displaystyle \frac{q_j}{{\theta }_j}},

where {\theta }_j=\arg s_j^{\prime }.

Further let \omega =(\delta +{\eta }^{\prime })t. If q_{\mathrm{d}}t\le \omega {\theta }_{\mathrm{d}}/\mathrm{1.8}, we have:

(28)\delta =\omega /t,\eta ={\eta }_{\mathrm{0}},\upsilon =\mathrm{1},

else

(29)\{\begin{array}{l}\delta =k\mu /\varphi \hfill \\ \eta =\widehat{p}-\mu \cot \varphi \hfill \\ \upsilon =q_{\mathrm{d}}/\mu \hfill \end{array},

where \{\begin{array}{l}k=\mathrm{1.6}+\mathrm{12}/(q_{\mathrm{d}}t+\mathrm{25})\hfill \\ \varphi =\mathrm{1.05}+\mathrm{1050}/\max (\mathrm{553}\mathrm{,}\mathrm{800}-q_{\mathrm{d}}t)\hfill \\ \mu =\left(\omega /t+{\eta }_{\mathrm{0}}-\widehat{p}\right)/(k/\varphi -\cot \varphi )\hfill \end{array} and more details about the theory and range of application can be consulted in the reference [30].

On the other hand, the compound Gauss integral is introduced to realize the Fourier inverse transform [6], i.e.:

(30)f\left(x,y\right)={\displaystyle \sum _{j=-m}^m{\displaystyle \sum _{k=-n}^n{w}_j{w}_k\overline{f}\left({\xi }_{xj},{\xi }_{yk}\right){\mathrm{e}}^{i\left(x{\xi }_{xj}+y{\xi }_{yk}\right)}}},

where {\xi }_{xj},{\xi }_{yk} are sample points and w_j, w_k are the associated weights.

Booker and Small [5-7] had demonstrated the accuracy and efficiency of these numerical methods mentioned above for solving Biot’s consolidation problems. When a thick layer is analyzed, numerical calculation may overflow, and this problem can be solved by using the technique proposed by Ai et al. [33]. Based on the final solutions in the transformed domain and the application of the inverse transform methods mentioned above, we have compiled the corresponding computational procedure, thus the actual solutions at an arbitrary point (whether the time or the space) can be obtained. All the results presented below are obtained under the assumption of permeable surface and fixed and impermeable base.

Suppose that an uniform rectangle vertical load f is applied on the ground surface, occupying the area, -a\le x\le a, -b\le y\le b, at time t={\mathrm{0}}^{+}, and then held constant.

The Laplace and the double Fourier transform of the load f yields:

(31)\overline{\mathit{f}}({\xi }_x,{\xi }_y,\mathrm{0},s)={\displaystyle \frac{\mathrm{1}}{{\mathrm{\pi }}^{\mathrm{2}}}}{\displaystyle \frac{f}{s}}{\displaystyle \frac{\sin {\xi }_{x}a}{{\xi }_x}}{\displaystyle \frac{\sin {\xi }_{y}b}{{\xi }_y}}.

A special case with \alpha =\mathrm{1.0} and γ = 1 is considered for a two-layered poroelasric medium, in which the pore fluid is incompressible and the permeability is isotropic. The computed results for a two-layered poroelastic medium by the proposed method are compared with those by Booker and Small [6]. The dimensionless surface settlement {\displaystyle \frac{G{u}_z}{bf}} at the center of the load is plotted against the time factor \tau ={\displaystyle \frac{\mathrm{2}G{k}_{v}t}{{\gamma }_{\mathrm{w}}{b}^{\mathrm{2}}}} in Fig.2. Obviously, the results from the solutions in this study match well with those of Booker and Small [6] using the finite layer method.

Consolidation of a single poroelastic layer

The influence of anisotropy of permeability on three-dimensional consolidation problems is investigated herein, in which the pore fluid is incompressible (\alpha =\mathrm{1.0}). Figures 3 and 4 present results and parameters of the three-dimensional consolidation problem of a single poroelastic layer. Different anisotropic ratios \gamma (0.05, 1, 5, 25) are selected to investigate the influence of the anisotropy of permeability on the consolidation behavior. Figure 3 shows that the surface settlements at different \gamma ratios are significantly different during the process of consolidation but the same at the initial and final times. Figure 4 also shows that for a given time factor \tau =\mathrm{0.1}, the dissipation of excess pore water pressure is significantly different at different \gamma ratios.

The influence of the compressibility of the pore fluid on the three-dimensional consolidation is investigated in this section, in which the permeability is isotropic (\gamma =\mathrm{1}). The compressibility parameters of the pore fluid a are selected as 1.0, 1.1, 1.2 and 1.3 to investigate the influence on the three-dimensional consolidation. The results of the dimensionless surface settlement {\displaystyle \frac{G{u}_z}{bq}} along with the dimensionless time factor \tau ={\displaystyle \frac{\mathrm{2}G{k}_{v}t}{{\gamma }_{\mathrm{w}}{b}^{\mathrm{2}}}} are presented in Fig. 5. Figure 5 shows that the increase of the compressibility factor \alpha promotes the initial surface settlement but does not affect the final surface settlement. The distribution of the dimensionless pore pressure \sigma /f along with the dimensionless depth z/h at a given time factor \tau =\mathrm{0.1} is presented in Fig. 6, which shows that the increase of the pore fluid compressibility factor \alpha accelerates the dissipation of the pore water pressure.

Consolidation of a layered poroelastic medium

A five-layered poroelastic medium with anisotropic permeability and compressible pore fluid is investigated to study the influence of anisotropic permeability and compressible pore fluid on the surface settlement. Three cases are examined.

Case 1: γ₁ ∶ γ₂ ∶ γ₃ ∶ γ₄ ∶ γ₅ = 1 ∶ 1 ∶ 1 ∶ 1 ∶ 1, α₁ ∶ α₂ ∶ α₃ ∶ α₄ ∶ α₅ = 1 ∶ 1 ∶ 1 ∶ 1 ∶ 1.

Case 2: γ₁ ∶ γ₂ ∶ γ₃ ∶ γ₄ ∶ γ₅ = 1 ∶ 1 ∶ 1 ∶ 1 ∶ 1, α₁ ∶ α₂ ∶ α₃ ∶ α₄ ∶ α₅ = 1.0 ∶ 1.3 ∶ 1.2 ∶ 1.1 ∶ 1.4.

Case 3: γ₁ ∶ γ₂ ∶ γ₃ ∶ γ₄ ∶ γ₅ = 5 ∶ 1 ∶ 0.05 ∶ 25 ∶ 1, α₁ ∶ α₂ ∶ α₃ ∶ α₄ ∶ α₅ = 1 ∶ 1 ∶ 1 ∶ 1 ∶ 1.

The other parameters of the five-layered poroelastic medium satisfy the following relationships: G₁ ∶ G₂ ∶ G₃ ∶ G₄ ∶ G₅ = 1 ∶ 4 ∶ 2 ∶ 2 ∶ 8; k_v1 ∶ k_v2 ∶ k _v3 ∶ k_v4 ∶ k_v5 = 4 ∶ 1 ∶ 8 ∶ 4 ∶ 2; h₁ ∶ h₂ ∶ h₃ ∶ h₄ ∶ h₅ = 2 ∶ 4 ∶ 1 ∶ 1 ∶ 2; and all the v are fixed to 0.25. The parameters and results are presented in Figs. 7–9, where h = h₁+h₂+h₃+h₄+h₅ is the total depth of the five-layered poroelstic medium. It can be seen from Fig. 7 that the anisotropy of permeability and the compressibility of pore fluid have great influence on the surface settlement of the poroelastic medium, and the early settlement of case 2 which has greater compressibility parameters is larger than those of other two cases, but at the later period, we get an opposite result. Figures 8 and 9 give the distribution of excess pore water pressure along depth at time factor \tau =\mathrm{0.01} and \tau =\mathrm{0.1}, respectively. The comparison of the two figures shows that the increase of parameters \alpha and \gamma is in favor of the dissipation of the pore water pressure, while their effects are different with time growth. Specifically speaking, the influence of the compressibility of pore fluid is more remarkable at the early stage of the consolidation (\tau =\mathrm{0.01}), however, the influence of the anisotropy of permeability becomes notable at the later period of the consolidation (\tau =\mathrm{0.1}). On the other hand, since the given time factor \tau =\mathrm{0.1} is already a large time parameter for the examined examples, the peaks of excess pore water pressure which may exist due to the inability of the water to be drained out in time at the beginning of the consolidation can not been found in Fig.9. Thus, it can be seen that the curves in Figs. 4, 6, and 9 are not so fluctuant as those in Fig.8.

Conclusions

This paper presents an uncoupled state space solution to three-dimensional consolidation of a layered poroelstic medium with anisotropic permeability and compressible pore fluid. Starting from the basic equations of poroelastic medium, and introducing intermediate variables, the state space equation usually comprising of eight coupled state vectors is uncoupled into two sets of equations of six and two state vectors in the Laplace-Fourier transform domain, therefore, further derivation and calculation is reduced. Based on the theorem of Cayley-Hamilton, the elements of the transfer matrix are readily got. Then the uncoupled state space solution of a layered poroelastic medium in transformed domain is obtained by using the transfer matrix method combined with the continuity conditions between adjacent layers and the boundary conditions. Numerical inversion is employed to get the actual solution. A two-layered poroelastic medium with isotropic permeability and incompressible pore fluid is investigated to verify the proposed solution and the results in this paper match well with those by Booker and Small [6]. Numerical analysis is carried out to investigate the effect of permeability anisotropy on the consolidation behavior of a single poroelastic layer, showing that anisotropy of permeability has a great influence on the consolidation behavior of the poroelastic medium. Numerical results also show that with increasing compressibility parameter, the poroelastic medium has greater initial surface settlement but the same final surface settlement, and the dissipation of the excess pore water pressure is accelerated. A five-layered soil with different anisotropic ratios and compressibility parameters is considered to study the influence of anisotropic permeability and compressible pore fluid on consolidation behavior of the layered soil, and it shows that, when other conditions are the same, the soil owning greater compressibility parameters would have larger settlement and degree of the dissipation of the pore water pressure at the early stage of the consolidation, but an opposite result is observed for the soil that have greater anisotropic ratios.

Compared to classical numerical methods, such as finite element or boundary element methods, the semi-analytical technique presented in this study has a higher efficiency of calculation, and can easily analyze the consolidation of a layered poroelastic medium with anisotropic permeability and compressible pore fluid.

References

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

[2] McNamee J, Gibson R E. Displacement functions and linear transforms applied to diffusion through porous elastic media. Quarterly Journal of Mechanics and Applied Mathematics, 1960, 13(1): 98-111

[3] McNamee J, Gibson R E. Plane strain and axially symmetric problem of the consolidation of a semi-infinite clay stratum. Quarterly Journal of Mechanics and Applied Mathematics, 1960, 13(2): 210-227

[4] Schiffman R L, Fungaroli A A. Consolidation due to tangential loads. In: Proceedings of the 6th International Conference on Soil Mechanics and Foundation Engineering. Montreal, 1965, 1: 188-192

[5] Booker J R, Small J C. Finite layer analysis of consolidation I. International Journal for Numerical and Analytical Methods in Geomechanics, 1982, 6(2): 151-171

[6] Booker J R, Small J C. Finite layer analysis of consolidation II. International Journal for Numerical and Analytical Methods in Geomechanics, 1982, 6(2): 173-194

[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] Wang J G, Fang S S. The state vector solution of axisymmetric Biot’s consolidation problems for multilayered poroelastic media. Mechanics Research Communications, 2001, 28(6): 671-677

[9] Wang J G, Fang S S. State space solution of non-axisymmetric Biot consolidation problems for multilayered poroelastic media. International Journal of Engineering Science, 2003, 41(15): 1799-1813

[10] Ai Z Y, Han J. A solution to plane strain consolidation of multi-layered soils. Soil and Rock Behavior and Modeling. ASCE: Geotechnical Special Publication, 2006, 150: 276-283

[11] Ai Z Y, Cheng Z Y, Han J. State space solution to three-dimensional consolidation of multi-layered soils. International Journal of Engineering Science, 2008, 46(5): 486-498

[12] Ai Z Y, Cheng Z Y. Transfer matrix solutions to plane-strain and three-dimensional Biot’s consolidation of multi-layered soils. Mechanics of Materials, 2009, 41(3): 244-251

[13] Booker J R, Carter J P. Elastic consolidation around a point sink embedded in a half-space with anisotropic permeability. International Journal for Numerical and Analytical Methods in Geomechanics, 1987, 11(1): 61-77

[14] Mei G X, Yin J H, Zai J M, Yin Z Z, Ding X L, Zhu G F, Chu L M. Consolidation analysis of a cross-anisotropic homogeneous elastic soil using a finite layer numerical method. International Journal for Numerical and Analytical Methods in Geomechanics, 2004, 28(2): 111-129

[15] Singh J S, Rani S, Kumar R. Quasi-static deformation of poroelastic half-space with anisotropic permeability by two-dimensional surface load. Geophysical Journal International, 2007, 170(3): 1311-1327

[16] Booker J R, Carter J P. Withdrawal of a compressible pore fluid from a point sink in an isotropic elastic half space with anisotropic permeability. International Journal of Solids and Structures, 1987, 23(3): 369-385

[17] Chen G J. Consolidation of multilayered half space with anisotropic permeability and compressible constituents. International Journal of Solids and Structures, 2004, 41(16-17): 4567-4586

[18] Ai Z Y, Wu C. Analysis of plane strain consolidation for a multi-layered soil with anisotropic permeability and compressible constituents. Chineses Journal of Theoretical and Applied Mechanics, 2009, 41: 801-807

[19] Skempton A W. The pore-pressure coefficients A and B. Geotechnique, 1954, 4(4): 143-147

[20] Cheng A H D, Liggett J A. Boundary integral equation method for linear porous-elasticity with applications to soil consolidation. International Journal for Numerical Methods in Engineering, 1984, 20(2): 255-278

[21] Yue Z Q, Selvadurai A P S, Law K T. Excess pore water pressure in a poroelastic seabed saturated with a compressible fluid. Canadian Geotechnical Journal, 1994, 31(6): 989-1003

[22] Senjuntichai T, Rajapakse R K N D. Exact stiffness method for quasi-statics of a multi-layered poroelastic medium. International Journal of Solids and Structures, 1995, 32(11): 1535-1553

[23] Pan E. Green’s functions in layered poroelastic half-space. International Journal for Numerical and Analytical Methods in Geomechanics, 1999, 23(13): 1631-1653

[24] Wang H F. Theory of Linear Poroelasticity with Applications to Geomechanics and Hydrogeology. Princeton: Princeton University Press, 2000

[25] Ai Z Y, Cheng Z Y. Plane strain Biot’s consolidation of multi-layered soils with compressible constituents. In: Characterization, Monitoring, and Modeling of GeoSystems. ASCE: Geotechnical Special Publication, 2008, 179: 702-709

[26] Ai Z Y, Wang Q S. Axisymmetric Biot’s consolidation of multi-layered soils with compressible constituents. In: Characterization, Monitoring, and Modeling of GeoSystems. ASCE: Geotechnical Special Publication, 2008, 179: 678-685

[27] Ai Z Y, Wu C, Han J. Transfer matrix solutions for three dimensional consolidation of a multi-layered soil with compressible constituents. International Journal of Engineering Science, 2008, 46(11): 1111-1119

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

[29] Verruijt A.Displacement functions in the theory of consolidation of thermoelasticity. Zeitschrift für angewandte Mathematik und Physik ZAMP, 1971, 22: 891-898

[30] Talbot A. The accurate numerical inversion of Laplace transforms. Journal of the Institute of Mathematics and Its Applications, 1979, 23(1): 97-120

[31] Sneddon I N. The Use of Integral Transform. New York: McGraw-Hill, 1972

[32] Pastel E C, Leckie F A. Matrix Methods in Elasto-Mechanics. New York: McGraw-Hill, 1963

[33] Ai Z Y, Yue Z Q, Tham L G, Yang M. Extended Sneddon and Muki solutions for multilayered elastic materials. International Journal of Engineering Science, 2002, 40(13): 1453-1483

RIGHTS & PERMISSIONS

Higher Education Press and Springer-Verlag Berlin Heidelberg