Skip to content

Hardening Soil (MN)

The implementation of the Hardening Soil (HS) model uses an implicit (backward Euler) return-mapping framework with stress-dependent stiffness, a hyperbolic shear hardening law, and an elliptical volumetric cap.

Distinctively, this implementation utilises the Matsuoka–Nakai (MN) failure criterion to account for the intermediate principal stress effect (Lode angle dependency) via a shape factor \(\chi(\theta)\), and applies this factor consistently to both the shear (cone) and compression (cap) yield surfaces.

Contents

Overview

The model assumes non-linear elasticity and two distinct plasticity mechanisms:

  1. Shear mechanism: a non-associative hyperbolic hardening cone, modified by the Matsuoka–Nakai shape factor \(\chi\).
  2. Volumetric mechanism: an associative elliptical cap controlling plastic compression, also modified by \(\chi\) for consistency.

The integration is performed via a backward Euler return mapping scheme. If the trial state overshoots the frictional strength (based on the mobilised friction angle), the increment is integrated with sub-stepping. A tensile-apex (vertex) cut-off is enforced explicitly.

Stress Measures, Invariants and Lode Angle

Sign convention and mean stress

The solver stress components in stress follow an Abaqus-like sign convention (compression is negative). However, the mean stress measure used internally is positive in compression and computed as

\[ p = -\frac{\sigma_{11}+\sigma_{22}+\sigma_{33}}{3}. \]

Shifted mean stress (cohesive apex)

To account for cohesion \(c\), the stress origin is shifted to the tensile apex of the failure envelope:

\[ p^\prime = p + c \cot \varphi. \]

Matsuoka–Nakai shape factor

The MN shape factor \(\chi(\theta)\) is evaluated from the deviatoric invariants \(J_2\) and \(J_3\) of the current (trial) stress state. In the code, the Lode angle \(\theta\) is obtained from

\[ \sin(3\theta) = -\frac{3\sqrt{3}}{2}\,\frac{J_3}{J_2^{3/2}}, \]

with clamping to avoid numerical issues when \(J_2 \to 0\). The friction-angle dependent parameter is

\[ \delta = \frac{3-\sin\varphi}{3+\sin\varphi}, \]

and the shape factor is computed as

\[ \chi = \frac{\sqrt{3}\,\delta}{2 \sqrt{\delta^2-\delta+1}\,\cos\upsilon}, \]

where \(\upsilon\) is a function of \(\theta\). In the local return-mapping subroutines, \(\chi\) is passed in and treated as constant for that local solve.

Elasticity

The model uses isotropic linear elasticity, but with stress-dependent moduli following a power law. The stiffness update is based on a positive magnitude of the least compressive principal stress (minor principal stress magnitude)

\[ \sigma_3^{\ast} = \max\!\Big(-\max(\sigma_1,\sigma_2,\sigma_3),\ \sigma_{\min}\Big), \]

where \(\sigma_i\) are principal stresses in the solver sign convention (compression negative). The lower bound \(\sigma_{\min}\) prevents singular stiffness updates.

The unloading/reloading modulus is

\[ E_{ur} = E_{ur}^{\text{ref}} \left( \frac{\sigma_3^{\ast} + c \cot \varphi}{p_{\text{ref}} + c \cot \varphi} \right)^m, \]

and the initial tangent stiffness used in the hyperbolic hardening law scales identically:

\[ E_i = E_{i}^{\text{ref}} \left( \frac{\sigma_3^{\ast} + c \cot \varphi}{p_{\text{ref}} + c \cot \varphi} \right)^m. \]

The Lamé constants are computed from \(E_{ur}\) and \(\nu_{ur}\) and assembled into the isotropic elastic Jacobian.

Yield Surfaces

Shear Yield Function (Cone)

The shear yield function \(f_s\) is defined as the residual of the hyperbolic hardening law (Eq. 7.54 in Benz (2007)1), expressed in terms of the deviatoric stress measure \(q\) and the accumulated plastic shear strain \(\gamma^p\) (stored as statev(2)):

\[ f_s = \frac{3}{2} \frac{q}{E_i} - \left( \frac{3}{2} \frac{q}{E_{ur}} + \gamma^p \right)\Gamma_{\text{mob}}. \]

The mobilisation factor is

\[ \Gamma_{\text{mob}} = \frac{(\eta_{\text{mob}} - R_f \eta_{\varphi}) \sin \varphi_{\text{mob}}}{1 - \sin \varphi_{\text{mob}}}, \qquad \eta_{\varphi}=\frac{1-\sin\varphi}{\sin\varphi}, \qquad \eta_{\text{mob}}=\frac{1-\sin\varphi_{\text{mob}}}{\sin\varphi_{\text{mob}}}. \]

The mobilised friction angle is computed from the stress ratio including the MN factor:

\[ \varphi_{\text{mob}} = \arcsin\!\left(\frac{3q}{6\chi\,p^\prime + q}\right), \qquad p^\prime = p + c\cot\varphi, \]

with clamping in the code to avoid invalid \(\arcsin(\cdot)\) arguments and to avoid non-positive mobilisation.

Cap Yield Function

The cap yield function \(f_c\) closes the elastic domain along the isotropic axis. It is an ellipse in the \(p^\prime-q\) plane with Lode-dependent scaling of the \(q\)-axis via \(\chi\):

\[ f_c = \frac{q^2}{(\chi \alpha)^2} + (p^\prime)^2 - (p_p^\prime)^2 = 0, \]

where

  • \(\alpha\) is the cap aspect ratio parameter,
  • \(p_p\) is the (unshifted) preconsolidation stress stored as statev(3),
  • \(p_p^\prime = p_p + c\cot\varphi\) is the shifted preconsolidation stress.

Hardening Laws

  1. Shear hardening: driven by the accumulated plastic shear strain \(\gamma^p\) through the implicit definition of \(f_s\).

  2. Cap hardening: the preconsolidation stress \(p_p\) evolves with the cap plastic multiplier \(\lambda_c\) (the cap return mapping uses \(\Delta p_p = H_c\,\Delta\lambda_c\)). In the implementation, the stress-dependent hardening modulus is

\[ H_c = 2 H_{pp} \left( \frac{\sigma_3^{\ast} + c \cot \varphi}{p_{\text{ref}} + c \cot \varphi} \right)^m p^\prime, \qquad p^\prime=p+c\cot\varphi. \]

Flow Rules

  • Cap mechanism: associative (plastic potential equals \(f_c\)).
  • Shear mechanism: non-associative, using a Rowe-type stress–dilatancy approach via a mobilised dilatancy angle \(\psi_{\text{mob}}\).

First, a critical-state angle is computed:

\[ \varphi_{\text{cs}} = \arcsin\!\left(\frac{\sin\varphi-\sin\psi}{1-\sin\varphi\,\sin\psi}\right). \]

Then:

  • If \(\sin\varphi_{\text{mob}} - \sin\varphi_{\text{cs}} \ge 0\) (dilative branch),
\[ \psi_{\text{mob}}=\arcsin\!\left(\frac{\sin\varphi_{\text{mob}}-\sin\varphi_{\text{cs}}}{1-\sin\varphi_{\text{mob}}\sin\varphi_{\text{cs}}}\right). \]
  • Otherwise, a contractive branch is used in the code, based on \(\eta = |q/(p+c\cot\varphi)|\) and an exponential interpolation.

For the tensile-apex check, the code forms

\[ \tan\beta = \frac{6\sin\psi^\ast}{3-\sin\psi^\ast}, \qquad \psi^\ast = \min(\psi_{\text{mob}},\psi), \]

and enforces the apex if

\[ p_{\text{trial}} \le -0.97\,(q_{\text{trial}}\tan\beta + c\cot\varphi). \]

Integration Algorithm Notes

  • Elastic predictor: \(\sigma_{n+1}^{\text{trial}} = \sigma_n + D^e:\Delta\varepsilon\) using the current stress-dependent \(D^e\).

  • Sub-stepping trigger: if the trial mobilisation exceeds the friction angle, i.e. \(\varphi_{\text{mob}}^{\text{trial}} \ge \varphi\), the increment is sub-stepped. The number of subincrements is

    \[ n_{\text{sub}} = \max\!\left(\left\lfloor\frac{\|\Delta\varepsilon\|_2}{5\times10^{-6}}\right\rfloor,\ 10\right), \]

    and \(\Delta\varepsilon\) is split uniformly.

  • Return mapping sequence: depending on which yield functions are violated, the algorithm returns to the cone, the cap, or their intersection (“mixed”).

  • Tangent operator: the current implementation returns an elastic tangent (not a consistent algorithmic tangent) from all local return mappings; in sub-stepping the elastic tangent is accumulated over subincrements.

Internal Optimisation

The parameters \(\alpha\) and \(H_{pp}\) are internal cap parameters that are typically calibrated from:

  1. \(E_{oed}^{ref}\) (reference oedometer stiffness),
  2. \(K_0^{nc}\) (at-rest earth pressure coefficient for normally consolidated state).

An internal routine optimize_hs_internal_constants(...) is provided to compute \((\alpha, H_{pp})\) by matching a virtual oedometer response.

Objective

The routine runs a virtual oedometer loading path (plane strain-like, uniaxial strain condition) at the reference stress level and forms

\[ \mathbf{R}= \begin{bmatrix} E_{oed}^{\text{target}} - E_{oed}^{\text{sim}}(\alpha,H_{pp})\\ K_0^{\text{target}} - K_0^{\text{sim}}(\alpha,H_{pp}) \end{bmatrix}, \]

with

  • \(E_{oed}^{\text{sim}} = \Delta\sigma_{v}/\Delta\varepsilon_{v}\),
  • \(K_0^{\text{sim}} = \Delta\sigma_{h}/\Delta\sigma_{v}\).

Algorithm (as implemented)

  • Initial guesses: \(\alpha_0=1.20\), \(H_{pp,0}=2E_{oed}^{ref}\).
  • Jacobian: computed every iteration via central differences with a relative perturbation of \(1\%\).
  • Update: solved as a \(2\times2\) system using Cramer’s rule.

Convergence criterion (as coded)

The loop exits when either residual satisfies the relative tolerance (\(10^{-6}\)):

\[ |R_1|\le 10^{-6}\,E_{oed}^{ref}\quad \text{or}\quad |R_2|\le 10^{-6}\,K_0^{nc}. \]

Material Parameters

The parameters must be provided in the props array in the following order:

Index Symbol Unit Description
1 \(E_{50}^{ref}\) F/A Secant stiffness in standard drained triaxial test.
2 \(E_{oed}^{ref}\) F/A Tangent stiffness for primary oedometer loading (target for optimisation).
3 \(E_{ur}^{ref}\) F/A Unloading/reloading stiffness.
4 \(m\) - Power law exponent for stress-dependent stiffness.
5 \(c\) F/A Cohesion.
6 \(\varphi\) deg Friction angle (internally converted to radians).
7 \(\psi\) deg Dilatancy angle (internally converted to radians).
8 \(\nu_{ur}\) - Poisson’s ratio for unloading/reloading.
9 \(p_{\text{ref}}\) F/A Reference mean stress level.
10 \(K_0^{nc}\) - At-rest earth pressure coefficient (normally consolidated).
11 \(R_f\) - Failure ratio.
12 \(E_i^{ref}\) F/A Initial tangent stiffness.
13 \(\alpha\) - Cap aspect ratio parameter.
14 \(H_{pp}\) - Cap hardening parameter.

Automatic determination of cap parameters

If \(\alpha = 0\) and/or \(H_{pp} = 0\), numgeo determines the missing parameter(s) such that the model reproduces the target oedometric stiffness \(E_{oed}^{ref}\) and the target \(K_0^{nc}\) under oedometric loading at the reference stress level.

Internal State Variables (statev):

  1. void_ratio: void ratio (only meaningful if initialised at analysis start)
  2. gammapss: accumulated plastic shear strain measure \(\gamma^p\)
  3. pp: preconsolidation stress \(p_p\)
  4. p: current mean stress \(p\)
  5. q: current deviatoric stress \(q\)

References

  1. Thomas Benz. Small-Strain Stiffness of Soils and Its Numerical Consequences. PhD thesis, Inst. für Geotechnik, 2007. URL: https://www.igs.uni-stuttgart.de/dokumente/Mitteilungen/55_Benz.pdf. ↩