Numerical and analytical solutions of poroelastic problems

Verruijt, Arnold

doi:10.1680/jgere.15.00006

For a number of problems of consolidation of soils (or poroelasticity), analytical solutions and numerical solutions are presented. The numerical solutions are primarily based on the method of finite elements in space and a single finite difference in time, so that the solution is obtained in the form of a stepwise process. Because of the mathematical complications in deriving analytical solutions, it is to be expected that engineers will use mainly numerical methods of solution. For this reason, it is important to consider the accuracy of these methods by comparing them with analytical solutions. For the problems considered in this paper, it appears that the numerical results are very close to the analytical results.

Notation

C_f: compressibility of fluid
C_s: compressibility of solids
G: shear modulus
K: compression modulus
k: permeability
m_v: compressibility of soil
n: porosity
p: fluid pressure
q: load
r: radial coordinate
S: storativity
t: time
U: finite-element functional
z: vertical coordinate
α: Biot’s coefficient
γ_f: volumetric weight of fluid
Δ: increment
δ: dimensionless time parameter
ϵ: volumetric strain
σ: total stress
σ^′: effective stress

Introduction

The theory of poroelasticity has been developed on the basis of Terzaghi’s analysis of the delay in the deformation of a loaded clay sample due to the outflow of water from the pores (Terzaghi, 1923) and Biot’s generalisation to three dimensions and soils consisting of compressible particles and a compressible fluid (Biot, 1941). The linear theory of poroelasticity (or consolidation) has reached a stage where there is practically general consensus on the basic equations; see, for example, Detournay and Cheng (1993), Verruijt (1995), Lewis and Schrefler (1998), De Boer (2000), Wang (2000), Rudnicki (2001), Coussy (2004), Gambolati et al. (2006) and Cheng (2016). Unfortunately, there is no general agreement on the definitions and the notations of the basic physical parameters in the theory of poroelasticity or consolidation. In this paper, the variables and the notations are taken in agreement with the usual presentations of soil mechanics, as in the books by Terzaghi (1943), Lambe and Whitman (1969), Craig (1997) and others. This means, for instance, that compressive stresses are considered as positive. It may be noted that Biot (1941) used his own set of notations, particularly for the coupling variables, and in his later publications (Biot and Willis, 1957), he modified these again. The notations in the presentation of the basic equations in this paper are based on the works of Gassmann (1951), Geertsma (1957), Skempton (1960) and Bishop (1973). The term ‘poroelasticity’ for the mechanics of a porous elastic material was introduced by Geertsma (1957), noting the analogy of the basic equations with the theory of thermoelasticity. He also expressed the parameters of the theory into basic engineering parameters, such as the compressibilities of the fluid and the solid particles.

Although many analytical solutions of problems of poroelasticity have been published – see, for example, Wang (2000), Schrefler (2002), Schanz (2009) and Cheng (2016) – these usually are restricted to simple geometries and simple boundary conditions. For this reason, engineers have put much effort in the development of numerical models. In this paper, some of these models are discussed and compared with the results of analytical solutions, particularly for determining the accuracy of the numerical models.

Terzaghi’s classical problem

In Terzaghi’s original theoretical work, which was inspired by the measurements of delayed deformations in a compression test on a clay sample (see Figure 1), the pore fluid and the soil particles were both assumed to be incompressible, so that the only mechanism of deformation was a rearrangement of the particles. In modern presentations, these assumptions are no longer made, and this generalised theory will be used here. A number of one-dimensional (1D) problems will be considered first, with two-dimensional (2D) problems to be considered later.

Figure 1

View large Download slide

Terzaghi’s classical problem

The analytical solution

For 1D deformation of a homogeneous porous material, the basic equation is the storage equation.

α \frac{\partial ϵ}{\partial t} + S \frac{\partial p}{\partial t} = \frac{k}{γ_{f}} \frac{\partial^{2} p}{\partial z^{2}}

1

which is based on the conservation of mass of both the particles and the fluid and Darcy’s law for the flow of the fluid with respect to the solid particles. In this equation, p is the pore pressure, ϵ is the volume change of the porous material, α is Biot’s coefficient, S is the storativity of the pore space, k is the permeability coefficient (sometimes known as the hydraulic conductivity) of the porous material and γ_f is the unit weight of the pore fluid. The storativity can be expressed as

S = n C_{f} + (α - n) C_{s}

2

where C_f is the compressibility of the fluid and C_s is the compressibility of the particle material. Biot’s coefficient can also be related to C_s by the equation

α = 1 - C_{s} / C_{m}

3

where C_m is the compressibility of the porous medium, the inverse of its compression modulus K.

C_{m} = 1 / K

4

In the case of 1D deformation, the volume change equals the vertical strain. Assuming linear elastic behaviour, this can be related to the vertical effective stress by the equation

\frac{\partial ϵ}{\partial t} = {- m}_{v} \frac{\partial σ_{zz}^{'}}{\partial t} = {- m}_{v} (\frac{\partial σ_{zz}}{\partial t} - α \frac{\partial p}{\partial t})

5

where m_v is the confined compressibility of the porous medium

m_{v} = \frac{1}{K + (4 / 3) G}

6

where K and G are the compression modulus and the shear modulus of the porous medium.

Elimination of the volume strain rate ∂ϵ/∂t from Equations 1 and 5 gives

\frac{\partial p}{\partial t} = \frac{α m_{v}}{S + α^{2} m_{v}} \frac{\partial σ_{zz}}{\partial t} + \frac{k}{γ_{f} (S + α^{2} m_{v})} \frac{\partial^{2} p}{\partial z^{2}}

7

This is the general differential equation for 1D consolidation.

The classical problem is that at time t = 0, a vertical load of magnitude q is applied, and this load is maintained for t > 0. This means that for t > 0, the total stress σ_zz is constant, so that Equation 7 reduces to

t > 0 : \frac{\partial p}{\partial t} = c_{v} \frac{\partial^{2} p}{\partial z^{2}}

8

where c_v is the consolidation coefficient.

c_{v} = \frac{k}{γ_{f} (S + α^{2} m_{v})}

9

Equation 8 is Terzaghi’s basic equation of 1D consolidation under constant total stress.

The initial condition can be established by noting that at the moment of loading, there cannot yet have been any fluid loss from the soil, so that it follows from Equation 7 that

t = 0 : p = p_{0} = \frac{α m_{v}}{S + α^{2} m_{v}} q

10

which is the initial condition of the problem. If the fluid and the solid particles are incompressible, α = 1 and S = 0, then p₀ = q, as was assumed in the original work by Terzaghi.

The boundary condition at the bottom of the sample is

t > 0, z = 0 : \frac{\partial p}{\partial z} = 0

11

and the boundary condition at the top of the sample is

t > 0, z = h : p = 0

12

The solution of the problem can best be obtained using the Laplace transformation and its complex inversion integral (Churchill, 1972). The final result is

\frac{p}{p_{0}} = \frac{4}{π} \sum_{k = 1}^{\infty} \frac{{(- 1)}^{k - 1}}{2 k - 1} \cos [(2 k - 1) \frac{π}{2} \frac{z}{h}] \exp [- {(2 k - 1)}^{2} \frac{π^{2}}{4} \frac{c_{v} t}{h^{2}}]

13

This is a well-known result that can be found, with its derivation, in many basic textbooks on soil mechanics. It should be noted that in some presentations, the thickness of the sample is denoted by 2h and both the upper and lower boundaries are considered fully drained. The centre of the sample, then, is located at z = h, and because of the symmetry, there is no flow at this location. A graphical presentation of Equation 13 will be shown later, comparing it with some simple numerical solutions.

The finite-difference method

As an alternative for the analytical solution presented earlier, a numerical method may be considered. The simplest one is the finite-difference method, which was first applied to consolidation problems by Gibson and Lumb (1953) and Abbott (1960). This method is particularly convenient for 1D problems.

The basic differential equation is Equation 8

t > 0 : \frac{\partial p}{\partial t} = c_{v} \frac{\partial^{2} p}{\partial z^{2}}

14

with the boundary conditions

z = 0 : \frac{\partial p}{\partial z} = 0

15

z = h : p = 0

16

and the initial condition

t = 0 : p = p_{0}

17

In the finite-difference method, the partial differentiations in Equation 14 are approximated by finite differences. This gives

\frac{p (z, t + Δ t) - p (z, t)}{Δ t} = c_{v} \frac{p (z + Δ z, t) - 2 p (z, t) + p (z - Δ z, t)}{{(Δ z)}^{2}}

18

When writing p (z, t) = p_i (t), p (z + Δz, t) = p_i+1(t) and p (z − Δz, t) = p_i−1(t), this equation can be written as

p_{i} (t + Δ t) = p_{i} (t) + δ [p_{i + 1} (t) - 2 p_{i} (t) + p_{i - 1} (t)]

19

where

δ = \frac{c_{v} Δ t}{{(Δ z)}^{2}}

20

If the values of the pore pressure at time t = 0 are known, as they are from the boundary condition in Equation 17, for i = 0,…,n, the values at time t = Δt can be calculated for i = 1,…,n − 1. The value for i = n follows from the boundary condition at the top (Equation 15), and the value for i = 0 can be calculated from the boundary condition at the bottom (Equation 16). The simplest way to do that is to introduce an extra value just below this boundary, for i = −1, assuming that p₋₁ = p₁, and then calculating p₀ from the algorithm in Equation 19 for i = 0. The symmetry condition then ensures that the first boundary condition is indeed satisfied. After calculating the values at time t = Δt, the next step can be made to calculate the values at time t = 2Δt and so on.

It may be noted that the simplicity of the procedure is a consequence of taking the values of the pressure in the right-hand side of Equation 18 or 19 all at the initial moment of time t. If these are taken at time t + Δt, or at some intermediate time, a more complex algorithm will be obtained, in which the new values of the pressure have to be calculated from a system of linear equations. Such an approach is called an implicit procedure. The explicit procedure used here is the simplest, and usually sufficiently accurate, provided that the process is stable.

Stability

The numerical algorithm in Equation 19 is stable only if the value of δ is small enough. If this value is taken too large, the process will become unstable, leading to fluctuating values of ever-increasing magnitude; see, for instance, Fox (1962). The stability criterion can be derived by assuming a state of zero pore pressures and then requiring that this is maintained. If it is assumed that there is a small error, with p _i−1(t) = ϵ, p_i (t) = −ϵ and p _i+1(t) = ϵ, it follows from Equation 19 that

p_{i} (t + Δ t) = (1 - 4 δ) ϵ

21

This is always smaller than ϵ provided that δ > 0, which means that the time step must be positive, which is a trivial condition. The value of p_i (t + Δt) must also be larger than −ϵ, however, if the process is not to lead to ever-increasing fluctuations. This leads to the condition δ < 0·5, which means that the time steps must satisfy the condition

Δ t < \frac{{(Δ z)}^{2}}{2 c_{v}}

22

This condition should always be satisfied when using the algorithm in Equation 19.

Figure 2 shows a comparison of the analytical solution and the numerical solution using finite differences. The analytical solution is shown by the fully drawn curves, and the numerical solution is shown by small dots. It appears that the agreement is very good. Because of the relative simplicity of the numerical method, it is widely used. Unfortunately, the finite-difference method is not so well suited for a non-homogeneous material – for instance, a layered soil. For such problems, the finite-element method is more suitable.

Figure 2

View large Download slide

Terzaghi’s classical problem: analytical and finite-difference solutions. Solid lines indicate analytical solutions, while dots indicate numerical solutions

The finite-element method

The finite-element method is particularly convenient for problems involving variable soil properties as will be demonstrated in this paper. For the general principles of the finite-element method, see, for instance, Zienkiewicz (1977) and Bathe (1982).

The problem to be considered is that of a soil column, which is subdivided into a number of elements, which may have different properties. The height of the elements may also be variable; see Figure 3.

Figure 3

View large Download slide

Column of soil and distribution of pore pressure

Basic equations

The basic differential equation of 1D poroelasticity of a non-homogeneous soil with a constant vertical load is (see Equation 14)

m \frac{\partial p}{\partial t} = \frac{\partial}{\partial z} (\frac{k}{γ_{f}} \frac{\partial p}{\partial z})

23

where m is the effective compressibility of the soil, k is the permeability and γ _f is the volumetric weight of the fluid (water). In general, for a soil with a compressible fluid and compressible particles, the effective compressibility of the soil is m = S + α ²/[K + (4/3)G] = S + α ² m _v. For a soil with an incompressible fluid and incompressible particles, the effective compressibility is m = m _v = 1/[K + (4/3)G]. It should be emphasised that the basic differential equation in the finite-element method, Equation 23, applies to a soil with variable permeability and compressibility. This variability is automatically taken into account in the development of the numerical equations, without any need for a separate treatment of the continuity of the flow at an interface between two regions with different properties. This is a distinct advantage of the finite-element method over other numerical methods.

The most common boundary conditions are that at the two boundaries, either the pore pressure or the specific discharge is prescribed. It will be assumed here that these boundary conditions are

z = 0 : \frac{k}{γ_{f}} \frac{\partial p}{\partial z} = 0

24

z = h : p = 0

25

The initial condition is assumed to be

t = 0 : p = p_{0}

26

This initial condition applies in the case of a load q that is applied at time t = 0 and then remains constant, so that

p_{0} = \frac{α m_{v}}{S + α^{2} m_{v}} q

27

The most common method of solution is to calculate the values of the pore pressure p in successive time steps. The appropriate equation may be obtained by integrating the differential equation in Equation 23 over a time step Δt. For a certain point in the soil column, this gives

Δ t \frac{d}{d z} (\frac{k}{γ_{f}} \frac{d \bar{p}}{d z}) - m [p (t + Δ t) - p (t)] = 0

28

where $\bar{p}$ is the average value of the pore pressure during the time step. It is assumed that this can be written as

\bar{p} = (1 - ϵ) p (t) + ϵ p (t + Δ t)

29

where ϵ is an interpolation parameter, 0 < ϵ < 1. This parameter can be used to simulate various interpolation schemes. A value ϵ = 1 indicates a backward interpolation, because in that case, the value during the entire time step is determined by the value at the end of the time step. A value ϵ = 0 indicates a forward interpolation, because then the value during the entire time step is determined by the value at the beginning of the time step. Linear interpolation, which may seem the most realistic, corresponds to the value ϵ = 1/2. Actually, it has been shown by Booker and Small (1975) that the process is always stable if 1/2 ≤ ϵ ≤ 1. It is therefore suggested to always use ϵ = 1/2 or, even safer (and simpler), ϵ = 1.

It follows from Equation 29 that

p (t + Δ t) - p (t) = \frac{\bar{p} - p (t)}{ϵ}

30

and writing p(t) = p, Equation 28 then gives

\frac{d}{d z} (\frac{k}{γ_{f}} \frac{d \bar{p}}{d z}) - \frac{m}{ϵ Δ t} (\bar{p} - p) = 0

31

If the average pore pressure during a time step $\bar{p}$ can be calculated using Equation 31, the value at the end of the time step can be determined from Equation 30: $p (t + Δ t) = p (t) + \bar{p} - p (t) / ϵ$ ⁠, and then a next time step can be taken.

Variational principle

There are several methods for establishing the system of numerical equations for the finite-element method. Here the equations will be derived using a variational principle, which is rather simple for the case of 1D problems. The variational principle states that the solution of the differential equation in Equation 31 can be obtained as the minimum of the functional

U = \frac{1}{2} \int_{0}^{h} [\frac{k}{γ_{f}} {(\frac{d f}{d z})}^{2} + \frac{m}{ϵ Δ t} (f^{2} - 2 f p)] d z

32

where f is an arbitrary function of z, except that it is required that for z = h, f = 0, indicating that the boundary condition of a given pore pressure (see Equation 25) should be satisfied. The other quantities in Equation 32, such as the soil properties, are the same as in Equation 31.

An example of a possible function f (z) is shown in Figure 3, assuming linear interpolation in each element. The proof that the function $f = \bar{p}$ indeed leads to a minimum of the functional U is given in many publications. It will be omitted here.

Finite elements

The soil column of height h is subdivided into n finite elements, denoted as R_ℓ, ℓ = 1,…,n. The functional U can now be expressed as

U = \sum_{ℓ = 1}^{n} U_{ℓ}

33

where

U_{ℓ} = \frac{1}{2} \int_{ℓ} [\frac{k}{γ_{f}} {(\frac{d f}{d z})}^{2} + \frac{m}{ϵ Δ t} (f^{2} - 2 f p)] d z

34

It is assumed that element ℓ is located between the two points at levels z ₁ and z ₂. In this element, the function f (the approximation of the average pore pressure $\bar{p}$ ⁠) is assumed to be given by a linear interpolation.

z_{1} < z < z_{2} : f = a z + b

35

where

a = \frac{c_{1} f_{1} + c_{2} f_{2}}{z_{2} - z_{1}} b = \frac{d_{1} f_{1} + d_{2} f_{2}}{z_{2} - z_{1}}

36

with

c_{1} = - 1 c_{2} = 1 d_{1} = z_{2} d_{2} = {- z}_{1}

37

It follows that

z_{1} < z < z_{2} : \frac{d f}{d z} = a = \frac{c_{1} f_{1} + c_{2} f_{2}}{z_{2} - z_{1}} = \frac{1}{Δ z} \sum_{i = 1}^{2} c_{i} f_{i}

38

where

Δ z = z_{2} - z_{1}

39

Noting that the length of the integration interval from z ₁ to z ₂ is Δz, it now follows that the first part of the integral in Equation 34 is

U_{ℓ 1} = \frac{1}{2} \int_{ℓ} \frac{k}{γ_{f}} {(\frac{d f}{d z})}^{2} = \frac{1}{2} \frac{k}{γ_{f}} \frac{1}{Δ z} \sum_{i = 1}^{2} \sum_{j = 1}^{2} c_{i} c_{j} f_{i} f_{j}

40

For the evaluation of the second part of the integral in Equation 34, it may be noted that

f^{2} = a^{2} z^{2} + 2 a b z + b^{2}

41

Furthermore, the initial pore pressure p can be expressed by a similar linear interpolation as f.

p = \frac{c_{1} p_{1}}{z_{2} - z_{1}} z + \frac{d_{1} p_{1} + d_{2} p_{2}}{z_{2} - z_{1}}

42

where the coefficients c_i and d_i have the same meaning as before; see Equation 37.

The evaluation of the second part of the integral now requires to integrate the expressions for f ² and fp over the interval from z ₁ to z ₂. Without loss of generality, it may be assumed that the origin z = 0 coincides with the average value of z ₁ and z ₂, so that z ₁ + z ₂ = 0 The necessary integrals are

\int_{ℓ} d z = Δ z

43

\int_{ℓ} z d z = 0

44

\int_{ℓ} z^{2} d z = \frac{{(Δ z)}^{3}}{12}

45

Using these integrals, the second part of the integral in Equation 34 is found to be

U_{ℓ 2} = \frac{1}{2} \frac{m}{ϵ Δ t Δ z} \sum_{i = 1}^{2} \sum_{j = 1}^{2} [d_{i} d_{j} + \frac{{(Δ z)}^{2}}{12} c_{i} c_{j}] f_{i} (f_{j} - 2 p_{j})

46

From Equations 40 and 46, it follows that the total contribution of the element to the functional is

U_{ℓ} = \frac{1}{2} \sum_{i = 1}^{2} \sum_{j = 1}^{2} (P_{i j} + R_{i j}) f_{i} f_{j} - \sum_{i = 1}^{2} \sum_{j = 1}^{2} R_{i j} f_{i} p_{j}

47

where

P_{i j} = \frac{1}{2} \frac{k}{γ_{f}} \frac{1}{Δ z} c_{i} c_{j}

48

R_{i j} = \frac{1}{2} \frac{m}{ϵ Δ t Δ z} [d_{i} d_{j} + \frac{{(Δ z)}^{2}}{12} c_{i} c_{j}]

49

After summation over all elements, the final expression for the functional U will be of the same general form as Equation 47, with the summation now being over all nodes – that is, from i = 0 to i = n. The matrices P and R can be determined by a summation over all elements. In this (numerical) process, the soil properties (k and m) and the element thickness (Δz) may be different for each element. This shows that this finite-element method is particularly well suited for non-homogeneous materials and layered soils with variable layer thicknesses. The time step Δt may also be different in the course of the analysis, but it may also be given a constant value, so that the matrix R _ij is constant in time, which may be convenient for the numerical computations.

Finally, the minimum value of the functional must be determined. This minimum can be found by requiring that

\frac{\partial U}{\partial f_{k}} = 0

50

for all f _k values where the pore pressure is not prescribed by a boundary condition. Because the matrices P _ij and R _ij are symmetric, this leads to the system of equations

\sum_{j = 0}^{n} (P_{i j} + R_{i j}) f_{j} = \sum_{j = 0}^{n} R_{i j} p_{j}

51

From this system of linear algebraic equations, the values of f_i may be determined, which gives an approximation of the average pore pressure in the first time step. The pore pressure at the end of the time step can then be determined from Equation 30, $p (t + Δ t) = p (t) + \bar{p} - p (t) / ϵ$ ⁠, and then a next time step can be taken. By stepping forward in time, the values of the pore pressure can be obtained as functions of time.

The system of linear equations can be solved by one of the many available algorithms – for instance, by Gauss elimination. For large systems, involving many points, it is worthwhile to take into account the banded structure of the matrices P and R in this process. This avoids storing many variables that are always equal to zero, and it also avoids arithmetic operations with such variables. Actually, in the present case of 1D consolidation, there are only three unknowns in each equation (f _j−1, f_j and f _j+1). The system matrix then is said to be a tridiagonal matrix, for which effective simple methods of solution exist; see, for instance, Press et al. (1988).

Terzaghi’s classical problem

Some results obtained using a finite-element program are shown in Figure 4. The problem refers to a layer of thickness h, subdivided into 40 sublayers. The soil data have been chosen such that the consolidation coefficient is c = 1·0 m²/d, so that consolidation is practically finished if ct/h ² = 2 or t = 2 d. The results can be compared with the results obtained by an analytical solution and a finite-difference solution, as shown in Figure 2. Again, the agreement is very good.

Figure 4

View large Download slide

Terzaghi’s classical problem: analytical and finite-element solutions. Solid lines indicate analytical solutions, while dots indicate numerical solutions

Two-layered soil

In this section, the analytical solution of the problem of 1D consolidation of a soil consisting of two layers with different permeabilities is presented and compared with results of a finite-element analysis. For reasons of simplicity, the compressibilities of the soil layers, the fluid and the particles are assumed to be the same in the two layers. The analytical method of solution of the problem is similar to the solution of a problem of heat conduction in a composite slab, as discussed by Carslaw and Jaeger (1959: p. 323). Two cases will be considered: one in which the lower boundary is impermeable and the upper boundary is fully drained, and another one in which both boundaries are fully drained. The finite-element model is a simple 1D model, using linear variations of the pore pressure in each element, as presented in the previous section, but now with a variability in the permeability.

The problem

The problem to be considered is the consolidation of a soil consisting of two layers; see Figure 5. The properties of the two layers are their thicknesses h ₁ and h ₂ and permeabilities k ₁ and k ₂ and a uniform compressibility m _v for both layers. The boundary conditions for the first problem are that the lower boundary is rigid and impermeable and that the upper boundary is fully drained and free to deform. The load is a step load of magnitude q, applied at time t = 0, resulting in an initial uniform pore pressure in the soil, p = p ₀ = (αm _v q/S + α ² m _v).

Figure 5

View large Download slide

Two layers of soil

The pore pressures will start to dissipate at time t = 0, but at a different rate in the two layers, depending on their consolidation coefficients, defined as

c_{1} = \frac{k_{1}}{γ_{f} (S + α^{2} m_{v})} c_{2} = \frac{k_{2}}{γ_{f} (S + α^{2} m_{v})}

52

where γ _f is the unit weight of the pore fluid, which is supposed to be γ _f = 10 kN/m³; S is the storativity, defined as S = nC _f + (α − n)C _s; n is the porosity; C _f is the compressibility of the fluid; C _s is the compressibility of the solid particles; C _m is the compressibility of the porous medium, the inverse of its compression modulus (C _m = 1/K); and m _v is the confined compressibility of the porous medium, m _v = 1/[K + (4/3) G]. It is to be noted that all properties are the same for the two layers, except for the permeability.

It is the purpose of this section to determine the pore pressure distribution as a function of time t and the vertical coordinate z.

Basic equations

Following Carslaw and Jaeger (1959: p. 323), it seems the simplest to assume that the vertical coordinate z = 0 is located at the interface of the two layers; see Figure 5.

In the upper part of the soil, the differential equation is

z > 0 : \frac{\partial p_{1}}{\partial t} = c_{1} \frac{\partial^{2} p_{1}}{\partial z^{2}}

53

In the lower part, the differential equation is

z < 0 : \frac{\partial p_{2}}{\partial t} = c_{2} \frac{\partial^{2} p_{2}}{\partial z^{2}}

54

The boundary conditions at the interface are that the pressure and the specific discharge are continuous.

z = 0 : p_{1} = p_{2}

55

z = 0 : k_{1} \frac{\partial p_{1}}{\partial z} = k_{2} \frac{\partial p_{2}}{\partial z}

56

The bottom of the soil is assumed to be impermeable, so that

z = {- h}_{2} : \frac{\partial p_{2}}{\partial z} = 0

57

At the top of the soil, full drainage is assumed.

z = h_{1} : p_{1} = 0

58

The initial condition is that at time t = 0, the pressure throughout the soil is given, determined by the load.

t = 0 : p = p_{0} = α m_{v} q / (S + α^{2} m_{v})

59

Analytical solution

The analytical solution of the problem can be obtained using the standard procedures of the Laplace transform theory. The only complication is the inverse transformation, for which the zeroes of the Laplace transform must be determined, along the negative axis of the complex solution. This can be done by first investigating what the likely distance between these zeroes is and then determining their precise location by searching for the zero values in succession. The number of zeroes to be determined should be chosen as large as n = 100 to achieve sufficient accuracy.

Finite-element solution

The computations for the finite-element solution are simple, but care must be taken that the permeabilities are different in the two layers. Some results are shown in Figure 6. The thickness of the upper layer is 1 m, and its permeability is 10 m/d. The thickness of the lower layer is 2 m, and its permeability is 0·1 m/d. The compressibility of both layers is 0·1 m²/kN. The figure shows the pore pressures after 0·01, 0·1, 1 and 10 d. The results of the analytical solution are shown by a continuous line, and the results of the finite-element computations are shown by small dots. The contrast of a factor of 100 in the permeability causes the pore pressures in the upper layer to be reduced rapidly, whereas in the lower layer, the small permeability results in a considerable time delay. The agreement between analytical and numerical results is very good. This may be considered as confirming their respective accuracies.

Figure 6

View large Download slide

Two layers: analytical and finite-element solutions. Permeability contrast 100:1. Solid lines indicate analytical solutions, while dots indicate numerical solutions

Another example is shown in Figure 7 for the case where the permeability in the upper layer is a factor of 100 smaller than the permeability in the lower layer. In this case, the greater permeability in the lower layer causes the pore pressure in this layer to be practically uniform, and the time delay is mostly caused by the small permeability in the upper layer. Again, the agreement between analytical and numerical results is very good.

Figure 7

View large Download slide

Two layers: analytical and finite-element solutions. Permeability contrast 1:100. Solid lines indicate analytical solutions, while dots indicate numerical solutions

The practically uniform pore pressure in the lower layer for a high permeability in that layer has the interesting consequence that the actual value of this permeability is not so important as long as it is considerably higher than the permeability in the upper layer, the drained layer. This can be illustrated by repeating the computations for a permeability contrast between lower and upper layers equal to 10 000. The results of such computations are shown in Figure 8. These results cannot be practically distinguished from the results given in Figure 7.

Figure 8

View large Download slide

Two layers: analytical and finite-element solutions. Permeability contrast 1:10000. Solid lines indicate analytical solutions, while dots indicate numerical solutions

Double drainage

To illustrate further the accuracy of the finite-element method, a variant of the problem may be considered in which both the upper boundary and the lower boundary are drained. The results of the analytical and numerical solutions are shown in Figure 9. As before, the results of the analytical and the numerical solutions are very close.

Figure 9

View large Download slide

Two layers: analytical and finite-element solutions. Permeability contrast 1:100. Solid lines indicate analytical solutions, while dots indicate numerical solutions

Radial flow problems

In this section, a number of 2D problems of poroelasticity are considered, in particular problems in which the deformation is axially symmetric, using a radial coordinate r and a vertical coordinate z. The finite-element method was developed for the solution of problems of 2D poroelasticity around 1970, with contributions from many authors; see, for instance, Sandhu and Wilson (1969) and Christian and Boehmer (1970). The procedures have gradually been improved and generalised by many authors. The most common procedures now use quadrangular elements with isoparametric shape functions and a derivation of the numerical equations using linear time steps and a Galerkin-type method in which the shape functions are used as weight functions in the approximation of the differential equations that should be satisfied. For problems involving a large number of elements, the solution requires sophisticated mathematical techniques, such as iterative procedures with advanced pre-conditioners; see Ferronato et al. (2010). The problems considered in this paper are simple enough to use standard-solution methods.

Radially loaded cylinder

A simple example of a cylindrical sample loaded along its radial boundary by a uniform load is shown in Figure 10. The upper and lower boundaries are impermeable and rigid. The basic geometry consists of four rectangular elements, with nine nodes. In the finite-element program, the elements can be refined by subdividing each element in four parts, and this has been repeated until the final number of elements is 4096 and the number of nodes is 4255. The program calculates the distribution of the pore pressure p for ten values of time. The results are shown in Figure 11, together with the analytical solution (De Leeuw, 1965). The agreement appears imperfect, but the largest error is 1·2%, which can be considered as satisfactory. It may be observed that the Mandel–Cryer effect (Cryer, 1963; Mandel, 1953), with an overshoot of the pore pressure in the centre of the sample for small values of time, is correctly obtained in both solutions.

Figure 10

View large Download slide

Cylindrical sample with radial load

Figure 11

View large Download slide

Radially loaded cylinder: pore pressures. Solid lines indicate analytical solutions, while dots indicate numerical solutions

Horizontal displacements in pumped aquifers

In groundwater hydraulics, problems of non-steady flow in aquifers have been studied extensively. In the first publications (Jacob, 1940; Theis, 1935), it was assumed, perhaps inspired by the success of Terzaghi’s work, that the soil in the aquifers deformed only in the vertical direction and that the horizontal displacements are assumed to be zero. This resulted in considerably simpler mathematics. In many subsequent studies for more complicated geometries, consisting of layered systems of aquifers and aquicludes (layers of low permeability), the same schematisation of vertical displacements only was used. Actually, this assumption seems to be in disagreement with the flow regime in the aquifers, which is that the flow of the groundwater is mainly horizontal, so that the Darcian friction forces on the soil are also horizontal. Then, it is to be expected that horizontal displacements occur in the aquifers. A more realistic approach is to use Biot’s theory of poroelasticity (Biot, 1941), which allows for displacements in all directions. Solutions for this more refined theory have been developed (Bear and Corapcioglu, 1981; Verruijt, 1969), and field measurements have demonstrated that in pumped aquifers, horizontal displacements are indeed observed (Wolff, 1970). If the vertical total stress and the shear stress on horizontal planes are assumed to vanish, a relatively simple set of equations is obtained, with the coefficient of consolidation defined as $c_{v}^{'} = k [K + (1 / 3) G] / γ_{f}$ ⁠, which is considerably smaller than the classical coefficient c _v = k [K + (4/3)G]/γ _f. As a result, the consolidation process takes longer, by about a factor of 2. Actually, the results from this refined analysis are very similar to the ones obtained by the classical theory of Theis (1935) and Jacob (1940), but the timescales are different by about a factor of 2.

As an example, the flow in a confined aquifer of radius R = 300 m and thickness of H = 10 m is considered; see Figure 12. The properties of the soil in the aquifer are such that the permeability is k = 1 m/d and $c_{v}^{'} = 100 m^{2} / d$ ⁠. The discharge of the well is Q = 10π = 31·4159 m³/d, and the dimensionless variable q is q = Qγ _f /πkH [K + (1/3)G] = 0·01, so that qR = 3 m. The results of an analytical solution and a numerical solution for the horizontal displacements with finite elements for axial symmetry are shown in Figure 13. The agreement appears to be good.

Figure 12

View large Download slide

A well in a circular confined aquifer

Figure 13

View large Download slide

Well: horizontal displacements for t = 10, 100 and 1000 d. Solid lines indicate analytical solutions, while dots indicate numerical solutions

Conclusion

It can be concluded that in all the problems considered, the results of the finite-element solutions, using simple algorithms, are in good agreement with the analytical solutions. It is particularly satisfying that the solutions for a column consisting of two soils with different permeabilities are so close to the analytical solution. The finite-element method well deserves its popularity in engineering practice.

Graphic. Refer to the image caption for details.

REFERENCES

Abbott

MB

1960

One-dimensional consolidation of multi-layered soils

Géotechnique

10

4

151

-

165

Numerical and analytical solutions of poroelastic problems

Notation

Introduction

Terzaghi’s classical problem

The analytical solution

The finite-difference method

Stability

The finite-element method

Basic equations

Variational principle

Finite elements

Terzaghi’s classical problem

Two-layered soil

The problem

Basic equations

Analytical solution

Finite-element solution

Double drainage

Radial flow problems

Radially loaded cylinder

Horizontal displacements in pumped aquifers

Conclusion

REFERENCES

Email Alerts

Suggested Reading

Related Chapters

Recommended for you

Cited By

Sharing Unavailable