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.

Overview

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

Shear mechanism: a non-associative hyperbolic hardening cone, modified by the Matsuoka–Nakai shape factor $\chi$.
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)¹), 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

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.
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 cap constants $\alpha$ and $H_{pp}$ are internal model parameters. When automatic determination is requested, the implementation identifies both constants from the prescribed reference oedometric stiffness $E_{oed}^{ref}$ and normally consolidated earth-pressure coefficient $K_0^{nc}$. The procedure is a local nonlinear root-finding problem, rather than a global material-parameter optimisation.

For a trial pair

\[ \boldsymbol{x}= \begin{bmatrix} \alpha\\ H_{pp} \end{bmatrix}, \]

the routine performs a virtual oedometer test. Starting from zero stress, an axial strain increment

\[ \Delta\varepsilon_v=-5\times10^{-6} \]

is applied repeatedly while the lateral strains remain zero. The constitutive response is integrated until the vertical stress reaches the reference pressure $p_{ref}$. The response of the final increment provides

\[ E_{oed}^{sim}(\boldsymbol{x}) = \frac{\Delta\sigma_v}{\Delta\varepsilon_v}, \qquad K_0^{sim}(\boldsymbol{x}) = \frac{\Delta\sigma_h}{\Delta\sigma_v}. \]

The calibration residual is therefore

\[ \boldsymbol{R}(\boldsymbol{x})= \begin{bmatrix} E_{oed}^{ref}-E_{oed}^{sim}(\boldsymbol{x})\\[2mm] K_0^{nc}-K_0^{sim}(\boldsymbol{x}) \end{bmatrix}. \]

Scaled logarithmic formulation

The two unknowns differ substantially in magnitude and both must remain positive. The iteration is therefore carried out in logarithmic, dimensionless variables

\[ \boldsymbol{y}= \begin{bmatrix} \ln\alpha\\[1mm] \ln\!\left(H_{pp}/E_{oed}^{ref}\right) \end{bmatrix}, \qquad \boldsymbol{x}= \begin{bmatrix} 1\\ E_{oed}^{ref} \end{bmatrix} \odot\exp(\boldsymbol{y}), \]

where $\odot$ denotes component-wise multiplication. This parameterisation guarantees $\alpha>0$ and $H_{pp}>0$ throughout the iteration. The initial values are

\[ \alpha^{(0)}=1.20, \qquad H_{pp}^{(0)}=2E_{oed}^{ref}. \]

The residual equations are normalised by their target values,

\[ \boldsymbol{r}(\boldsymbol{y})= \begin{bmatrix} \dfrac{E_{oed}^{ref}-E_{oed}^{sim}}{E_{oed}^{ref}}\\[3mm] \dfrac{K_0^{nc}-K_0^{sim}}{K_0^{nc}} \end{bmatrix}, \]

and the scalar merit function used to globalise the Newton iteration is

\[ \Phi(\boldsymbol{y})= \frac{1}{2}\,\boldsymbol{r}^{\mathsf T}\boldsymbol{r}. \]

The convergence measure is the maximum component-wise relative residual,

\[ e=\|\boldsymbol{r}\|_{\infty} = \max\!\left( \frac{|R_E|}{E_{oed}^{ref}}, \frac{|R_K|}{K_0^{nc}} \right). \]

Convergence is reached for

\[ e\le 10^{-6}, \]

which is equivalent to requiring both relative residuals to satisfy the tolerance simultaneously.

Globalised Newton iteration

At every accepted iterate, the Jacobian of the scaled simulated response with respect to $\boldsymbol{y}$ is recomputed by central differences. With

\[ \boldsymbol{g}(\boldsymbol{y})= \begin{bmatrix} E_{oed}^{sim}/E_{oed}^{ref}\\ K_0^{sim}/K_0^{nc} \end{bmatrix}, \]

the implemented approximation is

\[ J_{ij}\simeq \frac{ g_i\!\left(\boldsymbol{y}+\tfrac{h}{2}\boldsymbol{e}_j\right) -g_i\!\left(\boldsymbol{y}-\tfrac{h}{2}\boldsymbol{e}_j\right) }{h}, \qquad h=10^{-2}. \]

The Newton search direction $\boldsymbol{p}$ follows from

\[ \boldsymbol{J}\boldsymbol{p}=\boldsymbol{r}. \]

The $2\times2$ system is solved by Gaussian elimination with partial pivoting. Before the line search, the direction is scaled, if necessary, such that

\[ \|\boldsymbol{p}\|_{\infty}\le\ln(10). \]

Thus, an undamped accepted update cannot change either physical parameter by more than one order of magnitude.

A backtracking line search subsequently seeks

\[ \boldsymbol{y}^{(k+1)}= \boldsymbol{y}^{(k)}+\lambda\boldsymbol{p}, \qquad \lambda=1,\frac12,\frac14,\ldots, \]

and accepts the first trial state satisfying the Armijo condition

\[ \Phi\!\left(\boldsymbol{y}+\lambda\boldsymbol{p}\right) \le \Phi(\boldsymbol{y}) +c_1\lambda\,\nabla\Phi^{\mathsf T}\boldsymbol{p}, \qquad c_1=10^{-4}. \]

The backtracking factor is $1/2$ and at most 40 reductions are attempted. If the Newton system is singular, the computed Newton direction is not a descent direction, or the Newton direction cannot be globalised by the line search, the algorithm retries with the negative gradient of the merit function,

\[ -\nabla\Phi=\boldsymbol{J}^{\mathsf T}\boldsymbol{r}, \]

subject to the same logarithmic step bound and Armijo line search. The calibration is terminated with a diagnostic if no admissible descent step can be found or if convergence has not been achieved after 200 accepted iterations.

Nature of the calibration

The line search and fallback direction improve numerical robustness, but the procedure remains a local root search. The implementation does not imply global uniqueness of $(\alpha,H_{pp})$ for arbitrary combinations of $E_{oed}^{ref}$, $K_0^{nc}$ and the remaining material parameters.

Convergence history

The calibration routine accepts an optional history_file argument. If supplied, the initial state and every accepted iterate are written to a CSV file; finite-difference evaluations and rejected line-search trials are not included. The recorded quantities comprise the two residuals, their relative values, the convergence measure $e$, the tolerance, $\alpha$, $H_{pp}$, the merit function $\Phi$, the accepted line-search step length and the number of backtracking reductions.

The following example illustrates the convergence obtained for the parameter set used in the accompanying verification study. Both scaled residuals decrease below the prescribed tolerance after four accepted Newton updates.

Figure 1: Convergence history of the automatic determination of $\alpha$ and $H_{pp}$. The two component-wise relative residuals and their maximum $e$ are shown together with the implemented tolerance of $10^{-6}$. Only accepted iterates are plotted.

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

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

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

Reference manual

How to use it
Benchmark manual

Validation
Small-strain extension

Hardening Soil (MN, Bricks)

References

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

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.