CONTINUUM MECHANICS

V. Z. PARTON

AXISYMMETRIC PROBLEM IN THE THEORY OF CONSOLIDATION OF FLUID-SATURATED COMPRESSIBLE POROUS MEDIA

(Presented by Academician V. V. Shuleikin, 5 XI 1964)

To solve the axisymmetric problem of consolidation of a half-space deformed by normal pressure (Fig. 1), we shall use the equations of M. Biot \((^1)\). The Poisson coefficient \(\nu\) (as in M. Biot’s work \((^2)\)) is assumed equal to zero, which leads to a somewhat overestimated expression for the deflection. However, the presence under real conditions of entrapped air, which is not taken into account within the framework of the present theory \((^1)\), increases the settlement, which in a certain sense justifies the assumption introduced.

The equations of M. Biot \((^1)\), which describe consolidation, in the axisymmetric case have the form \((\nu = 0)\):

Fig. 1

\[ \left(\Delta-\frac{1}{\rho^2}\right)u_\rho-\frac{\partial e}{\partial \rho}-\frac{1}{G}\frac{\partial \sigma}{\partial \rho}=0,\qquad \Delta w-\frac{\partial e}{\partial z}-\frac{1}{G}\frac{\partial \sigma}{\partial z}=0,\qquad \Delta e=\frac{1}{c}\frac{\partial e}{\partial t}, \tag{1} \]

where

\[ \Delta=\partial^2/\partial\rho^2+\rho^{-1}\partial/\partial\rho+\partial^2/\partial z^2;\qquad e=-(\partial u_\rho/\partial\rho+u_\rho/\rho+\partial w/\partial z). \]

Here \(u_\rho, w\) are components of the displacement of particles of the solid phase; \(G\) is the shear modulus; \(\sigma\) is the pore pressure; \(c=2Gk\) is the coefficient of consolidation; \(k\) is the permeability coefficient.

We shall use the representations of the components in terms of two functions, introduced by McNamee and Gibson \((^3)\):

\[ \begin{aligned} u_\rho&=-\partial E/\partial\rho+z\,\partial S/\partial\rho,\\ \sigma_{\rho\rho}&=2G\left[(\partial^2/\partial\rho^2-\Delta)E-z\,\partial^2S/\partial\rho^2+\partial S/\partial z\right],\\ w&=-\partial E/\partial z+z\,\partial S/\partial z-S,\\ \sigma_{zz}&=2G\left[(\partial^2/\partial z^2-\Delta)E-z\,\partial^2S/\partial z^2+\partial S/\partial z\right],\\ \sigma&=2G(\partial S/\partial z-\Delta E),\\ \sigma_{\theta\theta}&=2G\left[(\rho^{-1}\partial/\partial\rho-\Delta)E-\rho^{-1}z\,\partial S/\partial\rho+\partial S/\partial z\right],\\ \sigma_{z\rho}&=2G(\partial^2E/\partial\rho\partial z-z\,\partial^2S/\partial\rho\partial z). \end{aligned} \tag{2} \]

Here \(\sigma_{\rho\rho}, \ldots\) are the total stresses, whereas \(\sigma'_{\rho\rho}=\sigma_{\rho\rho}-\sigma\) are the effective stresses, acting only on the elastic skeleton \((^4)\). With this representation the equations (1) take the form

\[ \Delta(\Delta-c^{-1}\partial/\partial t)E=0,\qquad \Delta S=0. \tag{3} \]

Let us consider the problem of compression of a saturated porous medium bounded by the plane \(z=0\). We shall assume that the porous medium under consideration is deformed under the action of a pressure \(p(\rho)\), distributed symmetrically with respect to the \(z\)-axis and acting normally to the boundary. It is assumed that the fluid contained in the soil freely seeps to the surface (both on the loading area and outside it), so that the water pressure at the surface is constant and equal to atmospheric pressure (settlement of a structure on a specially laid drainage layer).

The required functions must satisfy the following boundary and initial conditions:

\[ \begin{aligned} &1^\circ.\quad \Delta E=0,\quad t=0,\quad (e=0).\\ &2^\circ.\quad E,\ S\to0,\quad z\to\infty.\\ &3^\circ.\quad 2G\,\partial^2E/\partial z^2=p(\rho),\quad t>0,\quad z=0.\\ &4^\circ.\quad \partial E/\partial z=0,\quad \rho\geqslant0,\quad z=0.\\ &5^\circ.\quad \partial S/\partial z-\Delta E=0,\quad \rho\geqslant0,\quad t>0,\quad z=0. \end{aligned} \tag{4} \]

Relations \(3^\circ, 4^\circ\) express the conditions of equality of the normal effective stress and the load on the surface, as well as the equality to zero of the shear stress. Condition \(5^\circ\) corresponds to the requirement of free seepage of the fluid onto the surface. Applying to (3) the two-dimensional Laplace–Hankel transform, we obtain

Fig. 2

\[ (d^2/dz^2-\xi^2)\,[d^2/dz^2-(\xi^2+p/c)]\,\widetilde E=0, \]
\[ (d^2/dz^2-\xi^2)\widetilde S=0, \tag{5} \]

\[ \widetilde E(\xi,z,p)=\int_0^\infty\int_0^\infty \rho J_0(\xi\rho)e^{-pt}E(\rho,z,t)\,d\rho\,dt, \]

\[ \widetilde S(\xi,z,p)=\int_0^\infty\int_0^\infty \rho J_0(\xi\rho)e^{-pt}S(\rho,z,t)\,d\rho\,dt . \tag{6} \]

The solution of (5), bounded at infinity, has the form

\[ \widetilde E=c_1\exp\left[-\sqrt{\xi^2+p/c}\,z\right]+c_2e^{-\xi z}, \qquad \widetilde S=c_3e^{-\xi z}. \]

Conditions \(3^\circ\)—\(5^\circ\) (4), which determine \(c_1,c_2,c_3\), in the new variables \((\xi,z=0,p)\) take the form:

\[ 2G\,d^2\widetilde E/dz^2=\overline p(\xi)/p, \qquad d\widetilde E/dz=0, \qquad d\widetilde S/dz-(d^2/dz^2-\xi^2)\widetilde E=0; \tag{7} \]

The solution of (5) satisfying (7) can be represented in the form

\[ \widetilde E=\frac{\overline p(\xi)}{2Gp}\,\varphi(\xi,p) \left( \exp\left[-\sqrt{\xi^2+\frac pc}\,z\right] -\frac{\sqrt{\xi^2+p/c}}{\xi}e^{-\xi z} \right), \]

\[ \widetilde S=-\frac{\overline p(\xi)}{2Gc\xi}\,\varphi(\xi,p)e^{-\xi z}, \tag{8} \]

where

\[ \varphi(\xi,p)= \frac{1}{\xi^2+p/c-\xi\sqrt{\xi^2+p/c}} = \frac cp\left(1+\frac{\xi}{\sqrt{\xi^2+p/c}}\right). \]

Let us compute the vertical deflection of the surface \(z=0\):

\[ \widetilde w(\xi,0,p)=-\widetilde S(\xi,0,p) = \frac{\overline p(\xi)}{2G\xi p} \left(1+\frac{\xi}{\sqrt{\xi^2+p/c}}\right). \tag{9} \]

Applying the inverse Laplace and Hankel transforms, we find

\[ w=\frac{1}{2G}\int_0^\infty \overline p(\xi)\,[1+\operatorname{erf}(\xi\sqrt{ct})]\,J_0(\xi\rho)\,d\xi, \tag{10} \]

\[ \overline p(\xi)=\int_0^\infty \rho p(\rho)J_0(\xi\rho)\,d\rho; \qquad \operatorname{erf}(x)=\frac{2}{\sqrt\pi}\int_0^x e^{-\alpha^2}\,d\alpha . \tag{11} \]

1. Let us first consider the simplest case, namely the case of a concentrated force \(P\) acting at the point \(z=0,\rho=0\)*.

If \(p(\rho)=P\delta(\rho)/2\pi\rho\), where \(\delta(\rho)\) is the Dirac delta function, then, according to (11), \(\overline p(\xi)=P/2\pi\). Substituting into (10), we obtain:

\[ w=\frac{P}{4\pi G} \left[ \frac{2}{\sqrt\pi} \int_0^{\sqrt{ct}}\int_0^\infty \xi e^{-(\alpha\xi)^2}J_0(\xi\rho)\,d\xi\,d\alpha + \int_0^\infty J_0(\xi\rho)\,d\xi \right] = \]

\[ = \frac{P}{4\pi G} \left[ \frac{2}{\sqrt\pi}\int_0^{\sqrt{ct}} \frac{1}{2\alpha^2}\exp\left(-\frac{\rho^2}{4\alpha^2}\right)\,d\alpha + \frac{1}{\rho} \right] = \frac{P}{4\pi G} \left( \frac{1}{\rho} + \frac{2}{\rho\sqrt\pi}\int_{\rho/2\sqrt{ct}}^\infty e^{-\beta^2}\,d\beta \right) = \]

\[ = \frac{P}{4\pi G\rho} \left[ 1+\operatorname{erfc}\left(\frac{\rho}{2\sqrt{ct}}\right) \right], \qquad \operatorname{erfc}(x)=\frac{2}{\sqrt\pi}\int_x^\infty e^{-\alpha^2}\,d\alpha =1-\operatorname{erf}(x). \tag{12} \]

* The stress \(\sigma_{zz}\) arising along the line of action of the force \((\rho=0)\) was found by Paria (\(^{5}\)) for a soil model taking into account the compressibility of the particles of the solid phase.

As should be expected, as \(t \to \infty\), when all the liquid is squeezed out, we obtain the elastic solution for a concentrated force:

\[ w=\frac{P}{2\pi G\rho}. \]

At the initial instant the entire load is applied to the elastic skeleton, and the deflection is of a purely elastic character: \(w_0=w\big|_{t=0}=P/4\pi G\rho\). The settlement due to consolidation is \(w_s=w-w_0\), i.e.

\[ w_s=\frac{P}{4\pi G\rho}\operatorname{erfc}\left(\frac{\rho}{2\sqrt{ct}}\right). \tag{13} \]

Fig. 3

In Fig. 2 the dependence \(4\pi G w_s/P=w_s(\rho)\) is given \((\sqrt{ct}\) is the parameter). The state at the origin of coordinates is physically impossible.

2. Let the load now act uniformly over the area of a circle of radius \(a\), and let the pressure be absent outside this area. Then

\[ w_s=\frac{P}{4\pi G} \iint_{x_1^2+y_1^2\le a^2} \frac{1}{\sqrt{(x-x_1)^2+(y-y_1)^2}}\, \operatorname{erfc}\left[ \frac{\sqrt{(x-x_1)^2+(y-y_1)^2}}{2\sqrt{ct}} \right]\,dx_1\,dy_1 . \tag{14} \]

Replacing \(e^{-\alpha^2}\) by a series, we obtain in polar coordinates

\[ w_s=\frac{P}{4\pi G} \left[ \int_0^a r\,dr\int_0^{2\pi} \frac{d\theta}{\sqrt{\rho^2+r^2-2\rho r\cos\theta}} -\frac{2}{\sqrt{\pi}}\sum_{k=0}^{\infty} \frac{(-1)^k}{k!(2k+1)(2\sqrt{ct})^{2k+1}} \times \right. \]

\[ \left. \times \int_0^a r\,dr\int_0^{2\pi} (\rho^2+r^2-2\rho r\cos\theta)^k\,d\theta \right]. \tag{15} \]

After evaluating the integrals, we find:

\[ w_s=\frac{P}{4\pi G} \left[ 4aJ\left(\frac{\rho}{a}\right) -2\sqrt{\pi}\sum_{k=0}^{\infty}\sum_{m=0}^{k} \frac{(-1)^k\rho^{2k+2}}{(2k+1)k!(m+1)(2\sqrt{ct})^{2k+1}} \binom{k}{m}^{2} \left(\frac{a}{\rho}\right)^{2m+2} \right], \tag{16} \]

where

\[ J\left(\frac{\rho}{a}\right)= \begin{cases} \displaystyle \int_0^{\pi/2} \sqrt{1-(\rho/a)^2\sin^2\alpha}\,d\alpha =E(\pi/2,\rho/a), & 0\le \rho/a\le 1, \\[2.2ex] \displaystyle \int_0^{\pi/2} \frac{d\alpha}{\rho/a+\sqrt{(\rho/a)^2-\sin^2\alpha}}, & 1\le \rho/a . \end{cases} \tag{17} \]

Let us rewrite (16), introducing the dimensionless parameters \(\rho/a=\xi\), \(a/\sqrt{ct}=\eta\):

\[ w_s=\frac{P}{G}\left(\frac{ct}{\pi}\right)^{1/2} \left[ \frac{\eta}{\sqrt{\pi}}J(\xi) -\sum_{k=0}^{\infty}\sum_{m=0}^{k} \frac{(-1)^k}{(2k+1)k!(m+1)2^{2k+2}} \binom{k}{m}^{2} \xi^{2k-2m}\eta^{2k+2} \right]. \tag{18} \]

The settlement curves, calculated by formula (18) on the “Ural” electronic digital computer, are presented in Fig. 3 as functions of \(\xi\) at the instants of time corresponding to \(\eta = a/\sqrt{ct} = {}^8/_2, {}^8/_3, {}^8/_4, {}^8/_5, {}^8/_6\), and are compared directly with the settlement \(w_{s\infty}=Pa/G\sqrt{\pi}\eta\), which would occur at the same instants of time if the entire surface of the half-space were loaded. These settlement values are represented by horizontal lines with ordinates \(1Pa/8G\sqrt{\pi},\ 2Pa/8G\sqrt{\pi},\ldots\). The slope of the surface-deflection curve at the boundary of the loading area \((\xi=1)\) is infinite \((\partial w_s/\partial \xi|_{\xi=1}=\infty)\), which is connected with the infinitely large stress arising here. It is interesting to note that at first the settlement occurs mainly because of seepage of water to the surface directly on the loading area, and here \(w_s \cong w_{s\infty}\). At subsequent instants of time the unloaded region restricts the settlement of the loading area because of elastic stresses arising between the two regions (the restraining effect of the unloaded region).

Fig. 4

The magnitude of the deflection at the center of the loading area can be found from (15), if one sets \(\rho=0\), in the form

\[ \frac{w_s}{w_{s\infty}} = \frac{\sqrt{\pi}}{2}\eta - \sum_{k=0}^{\infty} \frac{(-1)^k}{(2k+1)(2k+2)2^{2k+1}k!}\eta^{2k+2}. \tag{19} \]

Let us simplify the expression of the sum entering into (19):

\[ \varphi(x) = 2\sum_{k=0}^{\infty} \frac{(-1)^k x^{2k+2}}{(2k+1)(2k+2)k!}, \qquad \text{where } x=\frac{\eta}{2}, \qquad \varphi''(x) = 2\sum_{k=0}^{\infty} \frac{(-1)^k x^{2k}}{k!} = 2e^{-x^2}. \]

The solution of the resulting equation \((\varphi(0)=\varphi'(0)=0)\) is

\[ \varphi(x) = 2\int_{0}^{x}(x-\xi)e^{-\xi^2}\,d\xi = x\sqrt{\pi}\operatorname{erf}(x)+e^{-x^2}-1. \]

Returning to the original variable, we obtain

\[ \frac{w_s}{w_{s\infty}} = \frac{\sqrt{\pi}}{2}\eta \left[ 1-\operatorname{erf}\left(\frac{\eta}{2}\right) \right] - e^{-\eta^2/4} + 1 = \frac{\sqrt{\pi}}{2}\eta\,\operatorname{erfc}\left(\frac{\eta}{2}\right) - e^{-\eta^2/4} + 1. \tag{20} \]

This ratio, characterizing the restraining effect of the unloaded regions, is presented in Fig. 4 as a function of \(1/\eta=\sqrt{ct}/a\). The larger the size of the loading area, the later this effect appears.

The author expresses his gratitude to L. S. Frank for his attention to the work.

Institute of Mechanics
Academy of Sciences of the USSR

Received
4 XI 1964

REFERENCES

1. M. Biot, J. Appl. Phys., 12, 155 (1941).
2. M. Biot, J. Appl. Phys., 12, 426 (1941).
3. J. McNamee, R. E. Gibson, Quarterly J. Mech. and Appl. Math., 13, part 1, 1 (1960).
4. K. Terzaghi, Theory of Soil Mechanics, 1961.
5. G. Paria, J. Math. and Phys., 36, No. 4 (1958).