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.

Small-strain stiffness

This module implements the base Hardening Soil model with stress-dependent, but strain-history-independent, stiffness moduli. A small-strain stiffness extension based on the BRICK concept is available as a separate module; see Hardening Soil (MN, Bricks).

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. \]

This shift, \(\text{apex} = c \cot\varphi\), is precomputed once per material point and reused throughout (elasticity, yield surfaces, flow rule, vertex check).

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} = 0.1\) (in the solver's stress units) prevents singular stiffness updates as the stress state approaches the isotropic tensile apex.

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. \]

Both moduli share the same stress-dependent factor \(\left(\frac{\sigma_3^\ast+c\cot\varphi}{p_{\text{ref}}+c\cot\varphi}\right)^m\), which is evaluated once per stress point and reused for \(E_{ur}\), \(E_i\) and the cap hardening modulus \(H_c\) (see Hardening Laws).

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:

\[ \sin\varphi_{\text{mob}} = \frac{3q}{6\chi\,p^\prime + q}, \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\). The hardening increment follows \(\Delta\gamma^p = \Delta\lambda \cdot \tfrac{3}{2}\), i.e. a fixed hardening factor of \(3/2\) relates the plastic multiplier of the cone mechanism to the accumulated plastic shear strain.

  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),
\[ \sin\psi_{\text{mob}}=\frac{\sin\varphi_{\text{mob}}-\sin\varphi_{\text{cs}}}{1-\sin\varphi_{\text{mob}}\sin\varphi_{\text{cs}}}. \]
  • Otherwise (contractive branch), with \(\eta = |q/(p+c\cot\varphi)|\) and \(M_{cs}=\chi\,M_c\) the critical-state stress ratio,
\[ \sin\psi_{\text{mob}} = -\frac{1}{10}\left[ M_{cs}\left(\frac{\eta\,\sin\varphi_{\text{cs}}\,(1-\sin\varphi_{\text{mob}})}{M_{cs}\,\sin\varphi_{\text{mob}}\,(1-\sin\varphi_{\text{cs}})}\right)^{1/15} - \eta \right]. \]

In both branches, \(\sin\psi_{\text{mob}}\) is clamped to \([-1+\epsilon,\,1-\epsilon]\) to prevent the flow potential from inverting.

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\), evaluated at the stress state at the start of the increment (or, within sub-stepping, at the start of the current sub-increment).

  • Sub-stepping trigger: if the trial mobilisation exceeds the friction angle, i.e. \(\varphi_{\text{mob}}^{\text{trial}} \ge \varphi\), the increment is discarded and re-integrated with sub-stepping. 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. Within each sub-increment, the stress-dependent moduli \(E_{ur}\), \(E_i\) and \(H_c\) are recomputed at the current stress level before the local return mapping is invoked. - Return mapping sequence: depending on which yield functions are violated, the algorithm returns to the cone, the cap, or their intersection ("mixed"). If the return-mapped state still exceeds the (unmodified) Matsuoka–Nakai friction angle \(\varphi\), a dedicated failure-surface (or mixed-failure-surface) correction is invoked to pull the stress state back onto the strength envelope; in this correction \(\gamma^p\) is not hardened further for the cone part, and only the cap plastic multiplier continues to be updated in the mixed-failure case. - Tangent operator: the local return-mapping subroutines take the elastic Jacobian \(D^e\) as an input constant (intent(in)) and do not overwrite it with a consistent algorithmic tangent; the tangent returned to the caller is therefore always the current elastic operator, evaluated at the stress state (of the increment, or of the last sub-increment when sub-stepping is active). Users who use dds_dde for a Newton-type outer iteration should be aware that it is not the consistent tangent of the elasto-plastic local solve. - Numerical safeguards: stress components are checked for NaN at the end of the routine, upon which the subroutine terminates with error stop and prints the intermediate mean/deviatoric stresses and plastic state for diagnostics.

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.

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}\).

The algorithm starts with initial guesses of \(\alpha_0=1.20\), \(H_{pp,0}=2E_{oed}^{ref}\). The Jacobian is computed every iteration via central differences with a relative perturbation of \(1\%\) per component of \((\alpha, H_{pp})\). The \(2\times 2\) system is solved with Cramer's rule and the solution is updated by a full Newton step (no line search or damping). The loop exits once both residuals satisfy the relative tolerance (\(10^{-6}\)) simultaneously:

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

Because \(\chi\alpha\) enters the cap yield function only as a squared term \((\chi\alpha)^2\), the sign of the converged \(\alpha\) has no effect on the model response; the routine returns \(|\alpha|\).

Convergence

The optimisation loop runs for up to MAXITER = 10000 Newton iterations and does not fall back to a bracketing or damped scheme if the finite-difference Jacobian becomes singular or the iteration diverges; it exits early with a diagnostic message in that case. For material parameter combinations far from typical soil behaviour (e.g. very low \(R_f\), unrealistic \(K_0^{nc}\)), convergence should be checked by inspecting the returned \((\alpha, H_{pp})\) and, if necessary, by supplying non-zero starting values for \(\alpha\) and \(H_{pp}\) directly in props to bypass the automatic optimisation.

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\)

At least 5 state variables must be provided; the model does not use or require any further entries.

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. ↩