Short note on the Biot soil consolidation

Jos van Kan

March 3, 2005

Abstract

We introduce briefly a computational scheme for the time dependent Biot soil consolidation model. The scheme is equivalent to solving 6 Laplacian equations per time step. We argue that the scheme is O(Δt^{2}) in the pore pressure p and the dilatation div u in which u is the soil displacement vector.

1 Introduction

The Biot soil consolidation model is given by

−μΔ̃u − (λ + μ)grad div u + grad p = 0  (1)
\frac{∂}{∂t}(div u) − \frac{κ}{η}Δp = f(x, t)  (2)

with boundary conditions p = 0 on Γ and u = 0 on Γ. Initial conditions are given for the dilatation div u. To facilitate the calculations we introduce a special symbol φ for the dilatation div u which transforms the system to:

\frac{∂φ}{∂t} − \frac{κ}{η}Δp = f(x, t)  (3)
−μΔ̃u − (λ + μ)grad φ + grad p = 0  (4)
div u = φ  (5)

This also introduces an extra boundary condition

φ = \frac{∂u_{n}}{∂n} on Γ.  (6)

2 Numerical model

We assume that we have a natural staggered space discretization. (φ- and ppoints must coincide). We also assume that we have the quantities p, u and φ on time level n with sufficient accuracy and the straightforward way to calculate an update δp, δu and δφ would be the solution of the following system of simultaneous linear equations:

δφ − \frac{1}{2}Δt\frac{κ}{η}Δ δp = Δt[\frac{κ}{η}Δ p^{n} + f^{n+1/2}]  (7)
−μΔ̃ δu − (λ + μ)grad δφ + grad δp = 0  (8)
δφ − div δu = 0  (9)

The accuracy of such an update would be O(Δt^{3}) at best (that is the error of the mid point integration) and for that to happen the error in p^{n} should be no worse than O(Δt^{2}). Fortunately these two characteristics are compatible.

This simultaneous system is a bit unwieldy and we therefore try a simplification that won’t lose us the overall accuracy. First of all we note, by taking the divergence of the second equation that

−(λ + 2μ)Δδφ + Δδp = 0.  (10)

Therefore we may write the first equation of the system as:

δφ − \frac{1}{2}Δt\frac{κ(λ + 2μ)}{η}Δ δφ = Δt[\frac{κ}{η}Δ p^{n} + f^{n+1/2}]  (11)

with boundary condition

δφ = \frac{∂δu_{n}}{∂n} on Γ.  (12)

Unfortunately, we do not know δu_{n} yet, but we can be sure, that it is of order O(Δt). So if we take δφ^{∗} = 0 on Γ instead, we have an error of O(Δt) in the boundary conditions, but only an error O(Δt^{2}) in the solution δφ^{∗} of system (11) on the interior points, since those boundary condtions are multiplied by Δt in the equations.

We define δψ = δp − (λ + 2μ)δφ. From equation (10) we now could calulate δψ^{∗} with an accuracy of O(Δt^{2}), provided we had boundary values to δφ^{∗} with that accuracy. We know δφ^{∗} in the interior with an accuracy of O(Δt^{2}), hence we can extrapolate to the boundary without loss of accuracy. From this solution δψ∗ we can obtain δp^{∗} with O(Δt^{2}) accuracy.

We plug δp^{∗} and δφ^{∗} into equation (8) to obtain δu^{∗} with O(Δt^{2}) accuracy. Note that this system is uncoupled in the u components, hence is equivalent to two Laplacian equations. This yields δu^{∗} with O(Δt^{2}) accuracy. Using this in the boundary condition of equation (11) we obtain an O(Δt^{3}) accurate solution δφ. From this we extrapolate the boundary contions for δφ once more to obtain O(Δt^{3}) accurate boundary conditions and an O(Δt^{3}) accurate solution δψ to (10). From δψ we get an O(Δt^{3}) accurate δp.

This suffices to obtain the pressure p and the dilatation φ with an global O(Δt^{2}) accuracy for the cost of 6 Laplacian like equations per time step. To get the displacement u to the same accuracy would require an additional two Laplacian equations.

There may be a stability problem.