Geotechnical FEM verification: the KATAI 2D benchmark suite

A geotechnical finite element result is a number someone may build on. Before such a number deserves any trust, the software that produced it has to demonstrate — not claim — that it solves the equations it says it solves. This article is the record of how KATAI 2D does that: which problems it is tested against, where their reference solutions come from, and what the measured deviations are. KATAI 2D is software I write and maintain, so read this as a documented verification record, not a review.

One of the cases below is worked through in full detail elsewhere — inputs, formula, mesh-convergence band and the residual difference from the published PLAXIS value: reproducing a PLAXIS 2D benchmark to within 0.03%.

Two kinds of number appear below, and they are checked in two different places. The analytical reference values and every percentage are recomputed from their stated inputs by a script that runs on each build of this site, the same way the Mohr circle and grading curve articles are audited. The measured KATAI 2D values are quoted from the verification record of the open repository, where a named automated test re-asserts each of them, inside a declared tolerance band, on every build.

Verification is the narrow claimPermalink to “Verification is the narrow claim”#

Verification and validation are different questions. Verification asks whether the code solves its governing equations correctly — against solutions known exactly. Validation asks whether those equations describe the physical ground. Everything in this article is the first question. A code that reproduces Terzaghi’s consolidation series has demonstrated correct time integration of the consolidation equation; whether that equation captures a particular estuarine clay is a separate matter that no benchmark suite can settle.

The suite admits three kinds of oracle, in this order of preference:

1. A closed-form solution, stated in full in the test that uses it.
2. An independent computation path that shares no code with the solver being tested — for example, every sparse linear-solver backend is checked against a dense factorisation with a different storage layout, agreeing to $4.4\times10^{-16}$.
3. A published benchmark, with the primary source cited.

Self-consistency — the code agreeing with itself — does not count as any of these. Two structural rules keep the record honest: the benchmark inputs are checked into the repository as plain .k2d files, re-solved from the file on every build (so nothing is hand-picked after the fact), and the verification matrix with its bibliography is generated from declarations inside the tests themselves, with a build gate that fails when the two drift apart. As of version 0.8.1 the record holds 58 declared cases and 26 checked-in benchmark files under 153 automated tests, and continuous integration reruns the full suite on every push on a configuration with no proprietary component.

That growth is not padding. Between 0.6.0 and 0.8.0 the work was to take capabilities that were implemented and tested at the element and run them along the path a user actually takes — from a file, through the mesher, through the driver. Six were taken that way for the first time; five of the six were wrong, and each one had been returning a converged result with a green suite behind it. An interface drawn along a fixed boundary was welded shut and returned 5,401,612 kN/m where the closed form gives 60. Deactivating the soil pinned the beams standing in it and reported max|u| = 0.000000e+00 with no diagnostic at all. That is the argument for this kind of suite stated as a measurement rather than as a principle: the tests that pass are not the ones that tell you something.

Collapse loads: the sharpest testPermalink to “Collapse loads: the sharpest test”#

Limit loads are where a displacement finite element code has the least room to hide. The elastic range forgives; a collapse mechanism does not. Three classical solutions anchor this class.

Prandtl’s bearing factor. For a weightless Tresca soil ($\varphi = 0$) under a strip footing, the exact bearing factor is

$$N_c = 2 + \pi = 5.14$$

The corpus file loads the footing past collapse and takes the honest non-convergence point as the oracle: the load fraction that equilibrates is the limit load. The measured factor is 5.175, +0.6% above exact. The same problem on a structured 15-node mesh measures +1.1%.

Cox’s circular footing. Cox (1962) solved the indentation of a ponderable cohesive-frictional soil by the slip-line method; for the parameters of the benchmark ($c = 1.6$ kPa, $\varphi = 30°$, $\gamma = 16$ kN/m³, footing radius 1 m) the limit pressure is

$$p_{\max} = 141\,c = 225.6\ \text{kPa}$$

This is an axisymmetric problem — the integration is radius-weighted, so it exercises a different kinematic path than any plane-strain case. The flow rule is associated ($\psi = \varphi$), because a slip-line solution is the associated limit load; running it non-associated would mix a modelling difference into a verification number. KATAI measures 233.9 kPa, +3.7% from the checked-in file’s own 0.25 m mesh. The mesh study is part of the record: a 0.5 m mesh puts only two elements across the footing radius and over-predicts by about +9% — the classic coarse-mesh bearing bias, recorded in the same place as the agreement.

Davis & Booker’s heterogeneous clay. Davis & Booker (1973) give the collapse pressure of a smooth rigid strip footing on clay whose undrained strength grows linearly with depth,

$$p_{\max} = \rho\left[(2+\pi)\,c_0 + \tfrac{1}{4}\,B\,c_{\text{inc}}\right]$$

For $c(z) = 1 + 2z$ kPa and $B = 2$ m the published value is 7.80 kPa. This case pins the depth-profile machinery — strength and stiffness varying with depth — against a theoretical limit load. KATAI measures 7.91 kPa, +1.4%. (The strength profile is a Tresca envelope, the $\varphi = 0$ special case of the Mohr–Coulomb geometry derived here.)

Elastic settlementPermalink to “Elastic settlement”#

Giroud’s rigid strip. For a smooth rigid strip footing on an elastic layer, Giroud (1972) gives the force at a prescribed settlement $s$ as

$$F = \frac{2\,(1+\nu)\,G\,B\,s}{\rho}, \qquad \rho = 0.88$$

which for $G = 500$ kPa, $\nu = 1/3$, half-width $B = 1$ m and $s = 10$ mm evaluates to 15.15 kN/m. KATAI drives the footing by a prescribed displacement and reads the force from the support reactions — the discrete internal force $B^{T}\sigma$ at the fixed degrees of freedom, with no recovery smoothing — measuring 15.32 kN/m, +1.1%.

Gibson’s incompressible soil. Gibson (1967) derived the surprising result that a strip load on an incompressible half-space whose shear stiffness grows linearly from zero at the surface settles uniformly under the load, by

$$s = \frac{q}{2\,\mathrm{d}G/\mathrm{d}z} = 0.050\ \text{m}$$

for the benchmark’s $q = 10$ kPa and $\mathrm{d}G/\mathrm{d}z = 100$ kPa/m. The benchmark models a 4 m finite layer, which must settle less than the half-space; KATAI measures 0.0454 m on that layer. The case doubles as a discretisation stress test: at $\nu = 0.495$ a low-order displacement element locks volumetrically, and a stiffness vanishing at the surface amplifies any such error. The 15-node quartic triangles carry both without special treatment.

Two further elastic cases pin stress distribution rather than settlement: the integrated Boussinesq solution for a uniform strip load (+0.9% to +1.4% at depth) and Flamant’s line-load solution (+0.0% to +1.3%).

Strength reduction: the slope problemPermalink to “Strength reduction: the slope problem”#

A factor of safety by $\varphi$–$c$ reduction is defined as the number by which $c'$ and $\tan\varphi'$ must be divided to bring the slope to failure. Griffiths & Lane (1999) made this the standard finite-element treatment of slope stability, with non-convergence as the failure criterion, and their Example 1 — the homogeneous 2:1 slope with $c'/\gamma H = 0.05$, $\varphi' = 20°$ — is the reference here. Their paper reports a finite-element factor of 1.4 and cites 1.380 from the Bishop & Morgenstern (1960) charts for the same slope.

KATAI solves the checked-in file — whose initial procedure is the safety analysis, so loading the file and solving it is the reduction — and measures FoS = 1.384: +0.3% against the charts, −1.1% against their published finite-element value. The suite also asserts that the mechanism genuinely displaces (0.40 m at the reduced strength), so the factor comes from a real slip surface, not a numerical artefact. A second slope from the same paper, run the same way, measures 1.010 against a published multi-method consensus around 0.99 — Bishop 0.988, Spencer 0.987, and 0.997 by another finite-element code.

Time and flowPermalink to “Time and flow”#

Terzaghi’s consolidation series. One-dimensional consolidation has an exact series solution for the average degree of consolidation,

$$U(T_v) = 1 - \sum_{m=0}^{\infty} \frac{2}{M^{2}}\,e^{-M^{2} T_v}, \qquad M = \frac{\pi}{2}\,(2m+1)$$

The corpus column is sampled at four time factors and compared point by point:

Charny’s dam discharge. For unconfined flow through a rectangular dam with a seepage face, Charny’s exact theorem gives the discharge

$$q = \frac{k\,(h_1^{2} - h_2^{2})}{2L}$$

independently of the shape of the free surface. KATAI’s free-surface seepage solve reproduces it to +1.02%, with a global mass balance of $8\times10^{-14}$ — inflow and outflow agree to round-off, which is a statement about the discretisation’s internal consistency on top of the agreement with the theorem.

The undrained column. A confined column loaded undrained is a one-line closed form (Skempton’s limit: the pore water takes the total stress increment); KATAI reproduces settlement and mid-column effective stress to +0.0%. And beneath all the time-dependent cases sits the geostatic baseline: the $K_0$ procedure on a submerged block leaves a residual of order $5\times10^{-13}$ kPa — the initial state is equilibrium to machine precision, not approximately.

Seismic dynamicsPermalink to “Seismic dynamics”#

The dynamic solver is checked at two levels. A resonant column driven at its fundamental frequency $f_1 = V_s/4H$ is compared with the damped shear-beam closed form (the standard derivation is in Kramer 1996, ch. 7): peak surface displacement −0.6%, peak acceleration −0.2%.

Above the closed form sits a real record. The suite ships a digitisation of the El Centro 1940 NS accelerogram inside the benchmark file and pins its identity first — PGA 0.31882 g at 2.02 s, against the published ≈0.319 g at ≈2 s — then its 5%-damped response spectrum: peak spectral amplification 2.87× at 0.19 s, inside the published 2.0–3.5 band for stiff-site recordings, with the 3 s ordinate at 0.118 g. Identity before response: if the input record were wrong, every site-response number computed from it would be wrong in a way no closed form would catch.

Against an independent implementationPermalink to “Against an independent implementation”#

Closed forms verify the equations; an independent implementation locates the result among its peers. The PLAXIS 2D Validation Manual, Version 8 (Bentley Systems) publishes analytic references alongside PLAXIS 2D’s own results for a set of canonical problems, and four of those cases are rebuilt in the KATAI corpus exactly as the manual defines them. The protocol is stated so the numbers can be read fairly: the PLAXIS values are published values quoted unchanged — PLAXIS was not re-run — and the two programs do not share a mesh; what is identical is the problem definition and the closed-form reference both programs answer to.

Two readings are worth making explicit. On the elastic problems the programs agree with each other to well under one percent where the geometry is like-for-like — 0.5% apart on Giroud — and both sit below the half-space closed form on the finite Gibson layer, as a finite layer must. On the Cox collapse problem the two programs bracket the exact slip-line value from opposite sides, PLAXIS 2.5% below and KATAI 3.7% above. That pattern is the honest lesson of the whole exercise: against an exact limit load, a few percent of discretisation-dependent bias, of either sign, is the expected result class for any displacement finite-element formulation — the sign and size of the bias are properties of the discretisation, not of the mechanics. The full walkthrough, with every parameter and command-line transcript, is on the KATAI verification record.

What the record does not claimPermalink to “What the record does not claim”#

Twenty-seven cases verify twenty-seven specific behaviours. They do not verify the program in general, and nothing here validates a constitutive model against a real soil — that is a different kind of evidence with a different owner: the engineer holding the site data. The suite’s job is narrower and checkable: that the equations the manual states are the equations the code solves, within bands that are declared rather than implied, on inputs anyone can re-run:

katai validate tests/corpus/kv-fnd-014-davis-booker-strip-footing.k2d
katai solve    tests/corpus/kv-fnd-014-davis-booker-strip-footing.k2d

Every case on this page works that way — a plain input file, the public command line, and a test that fails the build if the number drifts. The suite also pins the result contract itself: the command line, the Python package and the in-process engine produce bit-identical results files from the same input. A verification record is not a property a program acquires once; it is re-earned on every build, which is the same reason this site audits its own arithmetic.

SourcesPermalink to “Sources”#

The verification matrix, its bibliography and every benchmark input file are in the open repository under docs/validation/ and tests/corpus/; the summary record is published at kataigeo.com/benchmarks. The scientific manual has a citable DOI: 10.5281/zenodo.20692282.

References:

• Prandtl, L. (1921). Über die Eindringungsfestigkeit (Härte) plastischer Baustoffe und die Festigkeit von Schneiden. ZAMM 1(1):15–20. doi:10.1002/zamm.19210010102 — the bearing-capacity wedge solution for a weightless Tresca half-space.
• Cox, A. D. (1962). Axially-symmetric plastic deformation in soils — II. Indentation of ponderable soils. International Journal of Mechanical Sciences 4(5):371–380. doi:10.1016/S0020-7403(62)80024-1 — the tabulated slip-line yield-point pressures for a circular punch on cohesive-frictional soil.
• Davis, E. H. & Booker, J. R. (1973). The effect of increasing strength with depth on the bearing capacity of clays. Géotechnique 23(4):551–563. doi:10.1680/geot.1973.23.4.551
• Giroud, J.-P. (1972). Tables pour le calcul des fondations, Vol. 1: Tassement. Dunod, Paris — the settlement tables the rigid-strip reference value comes from.
• Gibson, R. E. (1967). Some results concerning displacements and stresses in a non-homogeneous elastic half-space. Géotechnique 17(1):58–67. doi:10.1680/geot.1967.17.1.58
• Griffiths, D. V. & Lane, P. A. (1999). Slope stability analysis by finite elements. Géotechnique 49(3):387–403. doi:10.1680/geot.1999.49.3.387
• Bishop, A. W. & Morgenstern, N. (1960). Stability coefficients for earth slopes. Géotechnique 10(4):129–153. doi:10.1680/geot.1960.10.4.129
• Terzaghi, K. (1943). Theoretical Soil Mechanics. John Wiley & Sons, New York — the one-dimensional consolidation series.
• Charnyi, I. A. (1951). A rigorous derivation of Dupuit’s formula for unconfined seepage with seepage surface. Doklady Akademii Nauk SSSR 79:937–940 (in Russian) — the proof that the Dupuit discharge is exact for the rectangular dam.
• Kramer, S. L. (1996). Geotechnical Earthquake Engineering. Prentice Hall — ch. 7, the damped shear-beam solution.
• Imperial Valley earthquake, 1940-05-18, El Centro Terminal Substation, S00E (NS) component — the digitisation shipped with the suite records its provenance and download identity in the repository (tests/data/elcentro-1940-ns.md).
• PLAXIS 2D Validation Manual, Version 8. Bentley Systems — Sections 2.1, 2.2, 3.1 and 3.2, the source of the published PLAXIS 2D values quoted above; the four analytic reference values in the comparison table are as stated in that manual.