Scheme Accuracy

The theoretical analysis of a time integration scheme accuracy and stability is generally based on a simplified problem (Zienkiewicz et al., 1988). Let us con­sider diffusion phenomena restricted to the linear case. Introducing the discre — tised field (Eq. 11.15) into the constitutive equations gives Darcy law (Eq. 11.7) (neglecting here the gravity term for the sake of simplicity and using the more gen­eral pressure p in place of u) in the following form:

fi = — K dtp = — K (dtNL )pL (11.32)

p p

Similarly the storage law (linear case) can be re-written:

S = rp = rNLpL (11.33)

where r is a storage parameter (cf. Eq. 11.8). Neglecting source terms, the weak form of the balance equation, Eq. 11.21, then produces:

j [SSp — ft dt (8p)]dV

V

/

к

— Nl pl diNK 8 pKdV = 0

p

VV

Considering that nodal values are not affected by the integration, this becomes:

(11.35)

which is valid for any arbitrary perturbation 8p. Thus:

RklP l + KklPl = 0 (11.36)

which is a simple system of linear equations with a time derivative, a storage matrix R (of which Rkl is an element) and a permeability matrix K (of which KKL is an element). One can extract the eigenvalues of this system and so arrive at a series of scalar independent equations of similar form:

pL + a2LpL = 0 (no summation) (11.37)

where L represents now the number of the eigenmode with the eigenvalue aL al­though it will be omitted in the following. The exact solution for Eq. 11.35 is a decreasing exponential function of time, t:

p(t) = p(t0)e-2t (11.38)

This problem represents the damping of a perturbation for a given eigenmode. Numerically, the modelling is approximated and numerical errors always appear. If Eq. 11.36 is well modelled, any numerical error will be rapidly damped, if the error source is not maintained. Following this analysis, the whole accuracy and stability discussion may be based on these last scalar Eqs. 11.37 and 11.38.

Introducing the time discretisation (Eqs. 11.26 and 11.27) into Eq. 11.37 gives:

PB-PA + a2 [(1 — r)Pa + трв] = 0 (11.39)

which allows evaluation of the end-of-step pressure as a function of the pressure at the beginning of the step:

Pb = UpA

with the amplification factor, U:

1 — (1 — t) a2 At
1 + та2 At

To ensure the damping process of the numerical algorithm, which is the sta­bility condition, it is strictly necessary that the amplification factor remains lower than unity:

-1 < U < 1

This condition is always verified if т > У2, and conditionally satisfied otherwise:

This last equation is not easy to verify, as it depends on the eigenvalues, which are generally not computed. Therefore, for classical diffusion processes considered in geomaterials, the condition т > /2 is generally used.

It should be noted that the amplification factor becomes negative for large time steps, except for the fully implicit scheme. Then, the perturbed pressure de­creases monotonically in amplitude but with changes of sign. This may be ques­tionable for some coupled phenomena, as it could induce oscillation of the coupled problem.

Let us now consider the accuracy of the numerical schemes. Developing in Tay­lor’s series the exact and numerical solution allows a comparison:

1 2 1 3

Uexact = 1 — X + X — X + …

2 6

^numerical = 1 — X + вX2 — в2X3 + …X = a2At (11.44)

It appears that only the Crank-Nicholson scheme т = /2 has second order accu­racy properties. However this conclusion is limited to infinitesimal time steps. For larger time steps, as in most numerical models, the Galerkin’s scheme т = 2 /3 gives the optimal compromise and should generally be used.

The whole discussion related to the stability and accuracy of the proposed time numerical schemes was based on eigenmodes of a linear problem. Can we extrapo­late them to general problems? The eigenvalue passage is only a mathematical tool allowing consideration of scalar problems, and has no influence on our conclusions. Conversely, the non-linear aspects could sometimes modify our conclusions. How­ever, it is impossible to develop the analysis for a general non-linear problem, and the preceding conclusions should be adopted as guidelines, as they appear to be fruitful in most cases.

Updated: 22 ноября, 2015 — 5:00 пп