Linear Benchmark: Terzaghi’s Consolidation Problem

The first example is Terzaghi’s consolidation problem. It was introduced by Karl von Terzaghi in [234]. Even though the problem is one-dimensional, it will be solved within a two-dimensional domain. It is an instructive example and a good test for the Biot stabilization. The results pre-sented here are not new, but stated for the sake of completeness and comprehensiveness of the theory. Similar numerical analysis of Terzaghi’s problem can for instance be found in [95, 246].

The setup comprises a porous block (height 2h), in which deformation and flow are restricted to the vertical direction, see Figure 4.5. Plain strain conditions are assumed. The bottom of the porous block is fixed, while on the top a constant load q = 1 is applied. For the fluid, the up-per and lower boundaries are considered drained, i.e. a constant pressure p = 0 is enforced.

Inertia terms are neglected for both the skeleton and the fluid. However, the dynamically chang-ing porosity (with initial valueφ0 = 0.2) is considered via the instationary term in the balance of mass. Linear kinematics are used for this example. A likewise linear stress-strain relation (St.Venant-Kirchhoff law, see equation (2.40)) with a Young’s modulus E = 1.0·10⁶ and a Poisson’s ratio ν = 0.0 determines the constitutive behavior of the skeleton. An incompress-ible microscopic solid phase is assumed, which leads to the Biot relation (3.93) with b = 1.0 and1/N = 0.0. The isotropic permeability is chosen ask = 1.0·10⁻⁶. Under those conditions, the exact solution for the pressure field can be derived [242] as

p

See also Appendix A.3.1 for some explanations of this form of the exact solution. For discretiza-tion in space, 470 linear, 3-node elements (approximate element size h = 0.1) are used. The one-step-θscheme withθ = 1.0and a time step∆t = 1.0·10⁻⁴ is applied. The analytical

solu-Figure 4.5: Terzaghi consolidation problem: Schematic of geometry and boundary conditions.

tion and the results of the simulation at different time instances without the Biot stabilization are depicted in Figure 4.6. It can be seen, that the sudden application of loading induces a pressure

0 0.5 1 1.5 2

0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6

y coordinate

pressure

exact solution t=0.0001 t=0.01 t=0.1 t=0.3 t=0.5 t=1.0

(a) pressure

0 0.5 1 1.5 2

-1 -0.8 -0.6 -0.4 -0.2 0 0.2

y coordinate

normal stress exact solution

t=0.0001 t=0.01 t=0.1 t=0.3 t=0.5 t=1.0

(b) normal skeleton stress

Figure 4.6: Terzaghi consolidation problem: Evolution of fluid pressure and skeleton stresses.

Linear 3-node elements without Biot stabilization are used.

jump and unphysical oscillations in the numerical solution predominantly at early times (left Subfigure 4.6). At the beginning, the fluid carries the main part of the load. With time, the load-ing is redistributed to the skeleton, leadload-ing to completely drained conditions at zero pressure in the stationary state. The numerical complexity of the problem originates from the high pressure gradients. This is particularly critical at small values of the dimensionless timecvt/h² and small time steps: The smaller the time step is chosen, the higher the pressure gradient, that needs to be resolved, becomes. The skeleton stresses, i.e. σ +p^fI, exhibit similar oscillations at lower magnitude (right Subfigure 4.6). Note, that this example is particularly designed to amplify the

oscillations. For bigger time steps or moderate permeability, the problem is far less prone to instabilities. Further, it is worth mentioning, that at later times, the oscillations are almost com-pletely annihilated and the numerical solution fits the exact solution very well.

Now, a way to diminish the oscillations is presented. The oscillations can be smoothed by a large enough Biot stabilization, as depicted in Figure 4.7. It significantly improves the initial

0

Figure 4.7: Terzaghi consolidation problem: Evolution of fluid pressure and skeleton stresses.

Linear 3-node elements with Biot stabilization are used.

pressure and stress solution. However, the stabilization has major drawbacks for linear elements especially at later time instances. As the strong form of the residual is used, the stabilization is not consistent for linear interpolations. The stress divergence term∇₀ ·P in equation (4.134) cannot be represented (the smoothing mainly stems from the penalization of pressure gradients).

Due to this inconsistency, a high numerical damping is introduced into the system. Thus, the pressure and the stress tend much faster to a constant level than the exact solution predicts.

Actually, the unstabilized form is far more accurate at later times here. Yet, there are some reme-dies for this problem. The most straightforward way is to simply use elements, which are more consistent. Of course, this might lead to an increase of computational costs. For quadratic el-ements, the stress divergence within the strong residual is not vanishing and better results can be obtained, see Figure 4.8. There are still some oscillations at early times, this time most evi-dent in the skeleton stress solution. Note, however, that the stabilization parameter was probably not yet optimal, such that the results could be further improved. Due to the better consistency, also good results are achieved at later time instances. Another approach is the reconstruction of the stress divergence. Also for linear elements, information of the stress divergence can be obtained from the nodal displacements. For instance, a gradient field can be computed by solv-ing a least squares problem, see [127, sec. 4.4.1]. Alternatively, other projection methods can be applied, which even have proven superconvergent characteristics. Results obtained with such a method taken from [269, sec. 15.4] are illustrated in Figure 4.9. Even though the oscillations are still visible and also the dissipative contribution of the stabilization terms can be seen at later times, there is a significant improvement compared to the results without recovery of the stress divergence, cf. Figure 4.7. Other ways to handle the pressure oscillations are even more fundamental. Mathematically, the problem lies in the violation of the inf-sup-condition for

0

Figure 4.8: Terzaghi consolidation problem: Evolution of fluid pressure and skeleton stresses.

Quadratic 9-node elements with Biot stabilization are used.

0

Figure 4.9: Terzaghi consolidation problem: Evolution of fluid pressure and skeleton stresses.

Linear 3-node elements with Biot stabilization are used. The stress divergence was computed by a patch-wise recovery [269, sec. 15.4].

placement and pressure. Stabilization techniques as the Biot stabilization used here circumvent the inf-sup-condition. Clearly, another way is to fulfill it. Discontinuous Galerkin methods for porous media have been designed to have this ability [189]. Such methods are not considered in this thesis. Another way, which uses standard Lagrange finite elements, is to choose suitable function spaces to fulfill the inf-sub condition. The simplest example is a Taylor-Hood element for displacement and pressure, see [162] for a discussion. There, quadratic functions are used for the displacements and linear interpolation for the pressure. In this case, no Biot stabilization is necessary, but it implies some more complexity of the element formulation and additional computational costs. This is also not presented here.