Theory & formulation

What is solved, written down.

An open-source program can be read; that is not the same as being documented. This page states the equations behind each analysis type of KATAI 2D and the numerical methods used to solve them, with the verification case that pins each one. The constitutive laws have their own page; the archived Scientific Manual is at its DOI.

1

Effective stress and equilibrium

Soil is treated as a porous skeleton saturated by water. Following Terzaghi, the skeleton deforms in response to the effective stress — the total stress less the pore pressure — and the constitutive law relates effective stress to strain:

σ=σ′−pw m,m=[ 1    1    1    0 ]T,pw=psteady+pexcess\boldsymbol{\sigma} = \boldsymbol{\sigma}' - p_w\,\mathbf{m}, \qquad \mathbf{m} = [\,1\;\;1\;\;1\;\;0\,]^{\mathsf T}, \qquad p_w = p_{\mathrm{steady}} + p_{\mathrm{excess}}Terzaghi’s principle, on this page with stresses tension-positive and the water pressure p_w positive in compression. The pore pressure is the sum of a steady part, from the phreatic level or a flow analysis, and an excess part generated by undrained loading.
∇ ⁣⋅σ+ρ g=0in Ω,u=uˉ on Γu,σ n=tˉ on Γt\nabla\!\cdot\boldsymbol{\sigma} + \rho\,\mathbf{g} = \mathbf{0} \quad\text{in } \Omega, \qquad \mathbf{u} = \bar{\mathbf{u}} \ \text{on } \Gamma_u, \qquad \boldsymbol{\sigma}\,\mathbf{n} = \bar{\mathbf{t}} \ \text{on } \Gamma_tStatic equilibrium with prescribed displacements and tractions.

The initial state is established either by the K0 procedure, which sets σ′h = K0 σ′v from the weight of the overburden and generates no displacement, or by gravity loading for ground whose surface or layering is not horizontal. For a submerged block the K0 state satisfies equilibrium to a residual of about 5×10−13 kPa (KV-NUM-003); in axisymmetry the identity holds to 8×10−13 on 15-node elements (KV-CST-013).

2

Finite element discretisation

The weak form of equilibrium, with displacements interpolated from nodal values, gives the discrete balance of internal and external forces:

u≈N a,ε=B a,r(a)=∫ΩBTσ dΩ  −  ∫ΩNTρ g dΩ  −  ∫ΓtNTtˉ dΓ=0\mathbf{u} \approx \mathbf{N}\,\mathbf{a}, \quad \boldsymbol{\varepsilon} = \mathbf{B}\,\mathbf{a}, \qquad \mathbf{r}(\mathbf{a}) = \int_{\Omega}\mathbf{B}^{\mathsf T}\boldsymbol{\sigma}\,\mathrm{d}\Omega \;-\; \int_{\Omega}\mathbf{N}^{\mathsf T}\rho\,\mathbf{g}\,\mathrm{d}\Omega \;-\; \int_{\Gamma_t}\mathbf{N}^{\mathsf T}\bar{\mathbf{t}}\,\mathrm{d}\Gamma = \mathbf{0}Displacement interpolation and the residual of equilibrium.

Elements are isoparametric triangles with six nodes (quadratic displacement) or fifteen (quartic). In plane strain the strain vector is [εxx, εyy, γxy] with εzz = 0. In axisymmetry the model is the r–z half-section, the hoop strain is carried and every integral is weighted by the radius:

ε=[ εr, εz, γrz, εθ ]T,εθ=urr,dΩ=2πr dr dz\boldsymbol{\varepsilon} = \bigl[\,\varepsilon_r,\ \varepsilon_z,\ \gamma_{rz},\ \varepsilon_\theta\,\bigr]^{\mathsf T}, \quad \varepsilon_\theta = \frac{u_r}{r}, \qquad \mathrm{d}\Omega = 2\pi r\,\mathrm{d}r\,\mathrm{d}zAxisymmetric kinematics.

Meshes are generated by Delaunay refinement (Ruppert 1995) with a minimum-angle bound of 20°, below the ≈20.7° limit for which the algorithm is guaranteed to terminate. Structural elements — plates, anchors, geogrids, embedded beams — are line elements whose nodes coincide with those of the continuum mesh, and interfaces are zero-thickness elements with a Coulomb strength.

The linear systems are solved by sparse direct factorisation. The released binaries use the vendored Eigen library and contain no proprietary component; a build from source can use Intel oneMKL PARDISO instead. Every backend is checked against a dense LU factorisation on a code path that shares nothing with it: 4.4×10−16 (KV-NUM-002).

3

Nonlinear solution

The load of a phase is applied in increments, and each increment is iterated to equilibrium by Newton’s method:

KT(k) δa(k)=− r(a(k)),a(k+1)=a(k)+η δa(k),KT=∫ΩBT Dep B dΩ\mathbf{K}_T^{(k)}\,\delta\mathbf{a}^{(k)} = -\,\mathbf{r}\bigl(\mathbf{a}^{(k)}\bigr), \qquad \mathbf{a}^{(k+1)} = \mathbf{a}^{(k)} + \eta\,\delta\mathbf{a}^{(k)}, \qquad \mathbf{K}_T = \int_{\Omega}\mathbf{B}^{\mathsf T}\,\mathbf{D}_{ep}\,\mathbf{B}\,\mathrm{d}\OmegaNewton iteration with a line-search step length η and the tangent stiffness assembled from the constitutive tangent.
  • Stress integration. At every stress point the constitutive law is integrated over the strain increment by return mapping, with the increment subdivided under an error tolerance so that stress-dependent moduli are re-evaluated within the step. The tolerance is a per-phase input and is stored with the result. Before this was done, the Hardening Soil answer moved by 2.9 percentage points with the number of load steps; it now moves by 0.020% (KV-NUM-009).
  • Line search. A step that would increase the out-of-balance force is shortened. Where a monotone search stalls on a well-supported model, the increment is retried with a non-monotone window before it is cut back.
  • An honest ending. If no step along the Newton direction reduces the residual, the run stops and says stalled line search; it says mechanism only when the stiffness parameter confirms one. The distinction is tested (KV-NUM-012): a stalled search published as a bearing capacity is a wrong answer that looks like a result.

4

Convergence criteria

An increment is accepted when a family of errors is below the tolerated error, not one:

eglobal=∥r∥∥fref∥≤eˉ,CSP=∫ΔεT Δσ dΩ∫ΔεT De Δε dΩe_{\mathrm{global}} = \frac{\lVert \mathbf{r} \rVert}{\lVert \mathbf{f}_{\mathrm{ref}} \rVert} \le \bar{e}, \qquad \mathrm{CSP} = \frac{\displaystyle\int \Delta\boldsymbol{\varepsilon}^{\mathsf T}\,\Delta\boldsymbol{\sigma}\,\mathrm{d}\Omega}{\displaystyle\int \Delta\boldsymbol{\varepsilon}^{\mathsf T}\,\mathbf{D}^{e}\,\Delta\boldsymbol{\varepsilon}\,\mathrm{d}\Omega}The global force error — out-of-balance force over a reference load — and the current stiffness parameter: the work actually done in the increment over the work an elastic body would have done. CSP = 1 for an elastic step and tends to 0 as a mechanism forms; the reference load is weighted by it, which tightens the test where it matters.
elocal=∥σeq−σc∥Tmax⁡≤eˉat every stress pointe_{\mathrm{local}} = \frac{\lVert \boldsymbol{\sigma}^{\mathrm{eq}} - \boldsymbol{\sigma}^{\mathrm{c}} \rVert}{T_{\max}} \le \bar{e} \quad \text{at every stress point}The local error: the distance between the stress a point must carry for equilibrium and the stress its own constitutive law returns, normalised by a reference shear stress at that point.

Both are computed, both are applied, and a results file records which were met and how many stress points were inaccurate. With only the global criterion enforced, the two measures differ by a factor of 63 000 on the verification case; with both, settlements agree to 0.0048% (KV-NUM-013).

5

Staged construction

A phase changes the configuration — clusters and members switched on or off, loads, water, prescribed displacements — and the out-of-balance force between the old configuration and the new one is applied in steps:

fext(i)=fext(i−1)+ΣMstage (fextnew−fext(i−1)),ΣMstage:0→1\mathbf{f}^{(i)}_{\mathrm{ext}} = \mathbf{f}^{(i-1)}_{\mathrm{ext}} + \Sigma M_{\mathrm{stage}}\,\Bigl(\mathbf{f}^{\mathrm{new}}_{\mathrm{ext}} - \mathbf{f}^{(i-1)}_{\mathrm{ext}}\Bigr), \qquad \Sigma M_{\mathrm{stage}}: 0 \to 1The stage multiplier is the proportion of the out-of-balance force that has been applied. A phase may stop short of 1 by request.

For a linear material the response must be exactly proportional to the multiplier, and a stage taken in two halves must arrive where the single stage does; both hold to 10−12 (KV-EXC-002). Removing a soil cluster removes its weight and its stress; the heave of an excavated floor matches the closed form to +0.0% (KV-EXC-001).

6

Consolidation and coupled flow–deformation

Equilibrium of the mixture is coupled to continuity of the pore fluid (Biot 1941):

∇ ⁣⋅ ⁣(σ′−pw m)+ρ g=0,mTε˙+nKw p˙w+∇ ⁣⋅q=0,q=− k ∇ ⁣(y+pwγw)\nabla\!\cdot\!\bigl(\boldsymbol{\sigma}' - p_w\,\mathbf{m}\bigr) + \rho\,\mathbf{g} = \mathbf{0}, \qquad \mathbf{m}^{\mathsf T}\dot{\boldsymbol{\varepsilon}} + \frac{n}{K_w}\,\dot{p}_w + \nabla\!\cdot\mathbf{q} = 0, \qquad \mathbf{q} = -\,\mathbf{k}\,\nabla\!\left(y + \frac{p_w}{\gamma_w}\right)Biot’s equations: momentum balance of the mixture, mass balance of the pore water with porosity n and water bulk modulus K_w, and Darcy’s law in terms of the hydraulic head.

Displacement and pore pressure are discretised in space and the system is integrated in time by the backward Euler scheme, which is unconditionally stable. In semi-discrete form,

K a˙+L p˙=f˙,LTa˙+S p˙+H p=qext\mathbf{K}\,\dot{\mathbf{a}} + \mathbf{L}\,\dot{\mathbf{p}} = \dot{\mathbf{f}}, \qquad \mathbf{L}^{\mathsf T}\dot{\mathbf{a}} + \mathbf{S}\,\dot{\mathbf{p}} + \mathbf{H}\,\mathbf{p} = \mathbf{q}_{\mathrm{ext}}Stiffness K, coupling L, compressibility S and permeability H matrices. With an elastoplastic skeleton K is the tangent stiffness and each time step is iterated.

The scheme is first-order accurate, and the order is measured rather than assumed: time refinement of the Terzaghi column gives an observed order of 0.985 and 0.992 on two triplets of step sizes, and a grid-convergence band of ±0.25% on the file’s own 120 steps that contains the actual error of 0.20% (KV-NUM-010). A phase may also march until a target is reached; the stop time is then within +1.45% of the closed form, an error that halves with each doubling of the steps per logarithmic level (KV-CON-003).

The fully coupled phase type solves the same pair of equations together with the groundwater flow of the phase; in the saturated limit it reproduces the consolidation phase to +0.00000% (KV-STR-008).

7

Groundwater flow

∇ ⁣⋅ ⁣(kr(Se) k ∇h)=n ∂S∂t,h=y+pwγw\nabla\!\cdot\!\bigl(k_r(S_e)\,\mathbf{k}\,\nabla h\bigr) = n\,\frac{\partial S}{\partial t}, \qquad h = y + \frac{p_w}{\gamma_w}Flow in terms of the hydraulic head; the right-hand side vanishes in steady state. The relative permeability k_r falls from 1 below the phreatic surface according to the retention curve.

In unconfined flow the free surface is found as part of the solution. On a potential seepage face the boundary condition is itself unknown — zero pressure where water exits, no flow where it does not — and is resolved iteratively. The check that the result is a solution is the global mass balance: 8×10−14 on the Charny dam, where the discharge is within +1.02% of the exact value (KV-FLW-001).

Wells and drains are line conditions on discharge or head; a barrier across the flow is an interface that is impermeable or carries a hydraulic resistance. Each has its own closed-form check (KV-FLW-002, -003, -004).

8

Strength reduction

cr=c′F,φr=arctan⁡ ⁣(tan⁡φ′F),FoS=sup⁡ { F:r(a; cr,φr)=0 has a solution }c_r = \frac{c'}{F}, \qquad \varphi_r = \arctan\!\left(\frac{\tan\varphi'}{F}\right), \qquad \mathrm{FoS} = \sup\,\{\,F : \mathbf{r}(\mathbf{a};\,c_r, \varphi_r) = \mathbf{0} \ \text{has a solution}\,\}The factor of safety as the largest strength reduction factor for which equilibrium exists (Griffiths & Lane 1999).

The factor is located by a bracketing search on F, each trial being a full elastoplastic equilibrium solve from the same stress state. The search’s own resolution on the verification slope is 6.3×10−4, and within the useful range of tolerated errors the factor is bit-identical: 1.421021 at 10−3, 10−4, 10−6 and 10−8 (KV-NUM-007).

Mesh dependence with non-associated flow. When ψ < φ′ the limit load of the discrete problem is not unique and the shear band narrows with the elements. On a nested family of 142, 568 and 2272 elements the factor falls 1.4210 → 1.3867 → 1.3645 (KV-SLP-003). The engine reports this on every such run (K2D-A005).

9

Dynamics

M a¨+C a˙+fint(a)=fext(t),C=α M+β K\mathbf{M}\,\ddot{\mathbf{a}} + \mathbf{C}\,\dot{\mathbf{a}} + \mathbf{f}_{\mathrm{int}}(\mathbf{a}) = \mathbf{f}_{\mathrm{ext}}(t), \qquad \mathbf{C} = \alpha\,\mathbf{M} + \beta\,\mathbf{K}Equations of motion with Rayleigh damping, integrated in time by Newmark’s method.
α=2 ξ ω1 ω2ω1+ω2,β=2 ξω1+ω2,ξ(ω)=α2ω+β ω2\alpha = \frac{2\,\xi\,\omega_1\,\omega_2}{\omega_1 + \omega_2}, \qquad \beta = \frac{2\,\xi}{\omega_1 + \omega_2}, \qquad \xi(\omega) = \frac{\alpha}{2\omega} + \frac{\beta\,\omega}{2}Rayleigh coefficients for a target damping ratio ξ at the two circular frequencies ω₁ and ω₂. Between them the damping is lower than the target, outside them higher.
tn=ρ Vp u˙n,ts=ρ Vs u˙s,Vs=G/ρ,Vp=Eoed/ρt_n = \rho\,V_p\,\dot{u}_n, \quad t_s = \rho\,V_s\,\dot{u}_s, \qquad V_s = \sqrt{G/\rho}, \quad V_p = \sqrt{E_{\mathrm{oed}}/\rho}Viscous boundary tractions (Lysmer & Kuhlemeyer 1969). At a compliant base they absorb the downgoing wave while the input motion is applied as the upgoing one; at free-field lateral boundaries they act on the motion relative to a one-dimensional column.

A uniform layer on rigid rock driven at its fundamental frequency f1 = Vs/4H reaches the closed-form surface amplitude to −0.6% in displacement and −0.2% in acceleration (KV-DYN-002). The response spectrum of a computed or recorded motion is the peak response of a damped single-degree-of-freedom oscillator as a function of its period; for the El Centro 1940 record shipped with the suite the peak spectral amplification is 2.87 at 0.19 s, inside the published band (KV-DYN-003).

10

How far is a computed number from the exact solution of its own equations?

“+1.4% against Davis & Booker on a 0.5 m mesh” hides a question a reviewer is entitled to ask: how much of that 1.4% is the model, and how much is the mesh? A deviation that shrinks under refinement is discretisation; one that does not is the formulation, the domain size or the reference. KATAI 2D separates them with the procedure the verification literature has standardised (Roache 1994; Celik et al. 2008; ASME V&V 20), built into the engine. From three systematically refined solutions φ1 (finest), φ2, φ3 with cell sizes h1 < h2 < h3:

h=AtotalNelements,r21=h2h1,r32=h3h2,ε21=φ2−φ1,ε32=φ3−φ2h = \sqrt{\frac{A_{\mathrm{total}}}{N_{\mathrm{elements}}}}, \qquad r_{21} = \frac{h_2}{h_1}, \quad r_{32} = \frac{h_3}{h_2}, \qquad \varepsilon_{21} = \varphi_2 - \varphi_1, \quad \varepsilon_{32} = \varphi_3 - \varphi_2Representative cell size, from the mesh that was actually built, and the refinement ratios.
p=∣ln⁡∣ε32/ε21∣+q(p)∣ln⁡r21,q(p)=ln⁡ ⁣(r21 p−sr32 p−s),s=sign⁡ ⁣(ε32ε21)p = \frac{\bigl|\ln\lvert \varepsilon_{32}/\varepsilon_{21}\rvert + q(p)\bigr|}{\ln r_{21}}, \qquad q(p) = \ln\!\left(\frac{r_{21}^{\,p} - s}{r_{32}^{\,p} - s}\right), \quad s = \operatorname{sign}\!\left(\frac{\varepsilon_{32}}{\varepsilon_{21}}\right)Observed order of convergence, by fixed-point iteration; valid for unequal refinement ratios.
φext=r21 p φ1−φ2r21 p−1,GCIfine=Fs ∣(φ1−φ2)/φ1∣r21 p−1,Fs=1.25\varphi_{\mathrm{ext}} = \frac{r_{21}^{\,p}\,\varphi_1 - \varphi_2}{r_{21}^{\,p} - 1}, \qquad \mathrm{GCI}_{\mathrm{fine}} = F_s\,\frac{\lvert(\varphi_1 - \varphi_2)/\varphi_1\rvert}{r_{21}^{\,p} - 1}, \qquad F_s = 1.25Richardson extrapolation and the Grid Convergence Index: a band on the fine-grid value, with a factor of safety of 1.25 when the order is established by three grids.

The implementation is verified on triplets manufactured from φ = φexact + C hp: the order and the extrapolated value come back to 10−14 (KV-NUM-004). Applied to the engine, the procedure produced three findings worth stating, because none is flattering by construction:

  • On an elastic strip load the triplet is in the asymptotic range and the extrapolated value is closer to the closed form than the finest mesh (KV-NUM-005).
  • For consolidation the refinement axis is the time step, the observed order is 0.99, and the band contains the true error (KV-NUM-010).
  • For the Hardening Soil oedometer the mesh is not the axis the error is on — refining it changes the answer by round-off — so no band is published for that case, and the record says why (KV-NUM-009).

The full account is numerical-uncertainty.md in the repository.

11

The verification method

Code verification asks whether a program solves the equations it claims to solve. Three kinds of evidence are admitted in KATAI 2D, and a fourth is explicitly excluded:

  • A closed-form solution, stated in full inside the test — not called from the code under test.
  • An independent computation path that shares no code with the solver: a dense factorisation against the sparse one, a hand-built equivalent model against the feature.
  • A published benchmark, with the primary source and its locator cited.
  • Not admitted: self-consistency. The code agreeing with an earlier version of itself is a regression test, not verification.

Each of the 69 cases is declared in a machine-readable block inside the test that asserts it — quantity, source, expected value, tolerance band, oracle class — and the verification matrix and its bibliography are generated from those declarations. A gate in the suite fails the build when the document and the tests disagree, so the published record cannot drift from what is actually checked. The tolerance bands are stated together with the measured deviation: “2%, measured +0.6%” says both how tight the test is and how much margin it has.

The record is reproduced, with sources, on the verification page.

The formulation is public. So is the evidence for it.

KATAI 2D is free under Apache‑2.0, installs in one line and ships the input files of the benchmarks it is verified against. The fastest way to judge it is to rerun one.

Install KATAI 2D See the verification record Source on GitHub