Skip to content

Hardening Soil (HS / HSsmall) — Theory & Implementation Notes

This document describes the Hardening Soil model implemented in material_hardening_soil_ (Fortran). It follows the same structure used for other constitutive models in this project and is written to match the driver’s data layout and sign conventions.

Contents

Overview

The Hardening Soil model is a frictional plasticity model with: - stress-dependent elastic stiffnesses, - a hyperbolic shear hardening law that saturates at a fraction \(R_f\) of the Mohr–Coulomb strength, - a triaxial stress-dilatancy rule (mobilized dilatancy), - an optional volumetric cap in compression, - an optional HSsmall shear-modulus reduction at small strains.

The implementation targets general 3D stress states. The rate routine computes only rates (no state mutation). The driver advances stress and internal variables with explicit schemes and, if requested, assembles a numerical tangent by finite differences.

Notation and sign conventions

  • Internal Cauchy stress tensor \(\mathbf{T}\) is compression negative (codebase convention).
  • Mean effective stress \(p\) and deviatoric stress \(q\) are handled and stored as positive magnitudes: [ p=\frac{1}{3}\,\mathrm{tr}(-\mathbf{T})\ge 0,\qquad \mathbf{S}=\mathbf{T}+\tfrac{1}{3}\mathrm{tr}(\mathbf{T})\,\boldsymbol{\delta},\qquad q=\sqrt{\tfrac{3}{2}\,\mathbf{S}:\mathbf{S}}\ge 0 . ]
  • Principal effective stresses used in the dilatancy rule are compression-positive magnitudes: [ \sigma_1\ge\sigma_2\ge\sigma_3>0,\quad \sigma_i = -\lambda_i(\mathbf{T}) . ]
  • The Lode angle \(\theta\) is computed for diagnostics only.
  • Plastic volumetric strain rate \(\dot\varepsilon_v^p\) is positive in compression.

Elastic law (stress-dependent moduli)

Reference moduli at \(p_{\text{ref}}\) are scaled by the current confining magnitude \(\sigma_3\) (triaxial-compression surrogate):

\[ E_{50}=E_{50}^{\text{ref}} \left(\frac{c\cos\varphi+\sigma_3\sin\varphi}{c\cos\varphi+p_{\text{ref}}\sin\varphi}\right)^{m}, \quad E_{\text{oed}}=E_{\text{oed}}^{\text{ref}}\left(\cdots\right)^{m}, \quad E_{\text{ur}}=E_{\text{ur}}^{\text{ref}}\left(\cdots\right)^{m}. \]

Un/reloading moduli:

\[ G_{\text{ur}}=\frac{E_{\text{ur}}}{2(1+\nu_{\text{ur}})},\qquad K=\frac{E_{\text{ur}}}{3(1-2\nu_{\text{ur}})} . \]

Isotropic linear elasticity (with optionally reduced \(G\); see HSsmall):

\[ \mathsf{E}^{\text{el}} = 3K\,\mathsf{I}^{\text{sym}} +2G\left(\mathsf{I}^{\text{sym}}-\tfrac{1}{3}\,\boldsymbol{\delta}\otimes\boldsymbol{\delta}\right). \]

Shear hardening (hyperbolic law)

Mohr–Coulomb TXC failure surrogate:

\[ q_f = \frac{2\,(c\cos\varphi+\sigma_3\sin\varphi)}{1-\sin\varphi},\qquad q_a=\frac{q_f}{R_f},\quad 0<R_f<1 . \]

Hyperbola for mobilized shear strength versus accumulated plastic shear strain \(\varepsilon_q^p\):

\[ q_h(\varepsilon_q^p) = \frac{q_a\,\varepsilon_q^p}{\varepsilon_q^p + q_a/E_i},\qquad E_i=\frac{2E_{50}}{2-R_f}. \]

Perzyna-type plastic shear rate (drives to the hyperbola):

[ \dot{\varepsilon}_q^p = \frac{g_s\,\langle q-q_h\rangle}{\partial q_h/\partial \varepsilon_q^p}, \qquad \frac{\partial q_h}{\partial \varepsilon_q^p} = \frac{q_a2/E_i}{\left(\varepsilon_qp+q_a/E_i\right)^2}, ] with \(g_s>0\) a dimensionless gain and \(\langle\cdot\rangle\) the Macauley bracket.

Dilatancy rule (triaxial states)

Plastic volumetric rate is tied to the plastic shear rate through the mobilized dilatancy:

\[ \dot{\varepsilon}_v^p = \frac{3}{2}\,\sin\psi_m\,\dot{\varepsilon}_q^p . \]

Mobilized friction (Rowe form) from principal magnitudes:

\[ \sin\varphi_m = \frac{\sigma_1-\sigma_3}{\sigma_1+\sigma_3-2c\cot\varphi},\qquad \cot\varphi=\frac{\cos\varphi}{\sin\varphi} . \]

Critical-state friction associated with input \((\varphi,\psi)\):

\[ \sin\varphi_{\text{cv}} = \frac{\sin\varphi - \sin\psi}{1-\sin\varphi\,\sin\psi}. \]

Mobilized dilatancy rule:

  • If \(\sin\varphi_m < \tfrac{3}{4}\sin\varphi\): \(\psi_m=0\).
  • Else if \(\psi>0\): [ \sin\psi_m = \max!\left( \frac{\sin\varphi_m-\sin\varphi_{\text{cv}}}{1-\sin\varphi_m\,\sin\varphi_{\text{cv}}}, \,0\right). ]
  • Else (\(\psi\le 0\)): \(\psi_m=\psi\).

All sines are clamped to \((-1,1)\) for robustness.

Volumetric cap (compression)

A simple cap regularizes the compressive volumetric response:

\[ F_c = p - p_c,\qquad p_c=\text{cap size (state)}. \]

If \(F_c>0\), activate cap plasticity with hardening:

\[ \dot{\varepsilon}_{v,\text{cap}}^p = -\,\frac{g_c\,\mathrm{tr}\,\dot{\boldsymbol{\varepsilon}}}{1+H_c/K}, \qquad \dot{p}_c = H_c\,\dot{\varepsilon}_{v,\text{cap}}^p . \]

Here \(H_c = E_{\text{oed}}\) (clipped positive), \(K\) is the current bulk modulus, and \(g_c>0\) a dimensionless gain. The sign yields \(\dot{\varepsilon}_{v,\text{cap}}^p>0\) in compression when \(\mathrm{tr}\,\dot{\boldsymbol{\varepsilon}}<0\).

Rate formulation used by the integrators

Decompose the plastic part into deviatoric and volumetric contributions:

  1. Deviatoric plastic flow aligned with the current deviatoric stress: [ \dot{\boldsymbol{\varepsilon}}_s^p = \dot{\varepsilon}_q^p\,\widehat{\mathbf{S}}, \qquad \widehat{\mathbf{S}} = \begin{cases} \mathbf{S}/|\mathbf{S}|, & |\mathbf{S}|>0,\[2pt] \mathbf{0}, & \text{otherwise.} \end{cases} ]

  2. Volumetric plastic flow: [ \dot{\boldsymbol{\varepsilon}}vp=\tfrac{1}{3}\,\dot{\varepsilon}_vp\,\boldsymbol{\delta}, \qquad \dot{\varepsilon}_v^p = \dot{\varepsilon}^p + \dot{\varepsilon}_{v,\text{cap}}^p . ]}

  3. Elastic stress rate: [ \dot{\mathbf{T}} = \mathsf{E}^{\text{el}}:\Big( \dot{\boldsymbol{\varepsilon}}

  4. \dot{\boldsymbol{\varepsilon}}_s^p
  5. \dot{\boldsymbol{\varepsilon}}_v^p \Big). ]

The rate routine returns \(\dot{\mathbf{T}}\), \(\dot{\varepsilon}_q^p\), \(\dot{\varepsilon}_v^p\), and \(\dot{p}_c\).

Numerical integration

Modified-Euler (Heun) with error control

Per sub-step \(h\):

  1. Evaluate rates at the start; predict stress.
  2. Re-evaluate at the predicted state; apply Heun corrector for stress.
  3. Update hardening variables with trapezoid rule.
  4. Estimate stress error from start vs. end rates; accept or shrink \(h\).
  5. Optionally assemble numerical tangent over the accepted sub-step.

Euler–Richardson with adaptive step size

Per sub-step \(h\):

  1. Forward Euler trial from start rates.
  2. Mid-point evaluation; Richardson extrapolation for stress.
  3. Local stress error from trial vs. extrapolated stress; accept or shrink.
  4. Trapezoid update of hardening variables.
  5. Optionally assemble numerical tangent at the accepted end state.

Numerical Jacobian

  • num_diff = 0: elastic stiffness only (HSsmall aware).
  • num_diff = 1: forward differences \(O(h)\).
  • num_diff = 2: central differences \(O(h^2)\).

A small diagonal “floor” is added to avoid singular tangents near stress-free states.

Optional features

HSsmall small-strain stiffness

If \(G_0^{\text{ref}}>0\) and \(\gamma_{0.7}>0\), the shear modulus is reduced from \(G_0\) to \(G_{\text{ur}}\) using a hyperbolic decay versus an equivalent shear strain \(\gamma_{\text{eq}}\):

\[ G_0 = G_0^{\text{ref}}\left(\frac{p}{p_{\text{ref}}}\right)^m,\qquad \gamma_r = \frac{\gamma_{0.7}}{0.7},\qquad G = \max\!\left(\frac{G_0}{1+|\gamma_{\text{eq}}|/\gamma_r},\,G_{\text{ur}}\right). \]

Minimum mean effective stress

After each accepted sub-step a projection enforces \(p\ge p_{\min}\):

\[ p<p_{\min}\ \Rightarrow\ \sigma_{ii}\leftarrow \sigma_{ii}+(p-p_{\min}),\ i=1,2,3 . \]

This isotropic shift preserves the deviatoric stress.

Material parameters (order in input)

Exactly the order expected by the implementation:

Idx Symbol Name Unit Notes
1 \(\varphi\) friction angle deg converted to rad
2 \(\psi\) dilatancy angle deg converted to rad
3 \(c\) effective cohesion kPa \(c_{\text{eff}}\)
4 \(E_{50}^{\text{ref}}\) secant modulus at 50% strength kPa at \(p_{\text{ref}}\)
5 \(E_{\text{oed}}^{\text{ref}}\) oedometer modulus kPa at \(p_{\text{ref}}\)
6 \(E_{\text{ur}}^{\text{ref}}\) un/reloading modulus kPa at \(p_{\text{ref}}\)
7 \(m\) stress dependency exponent – HS exponent
8 \(\nu_{\text{ur}}\) Poisson ratio (un/reload) – \(0<\nu_{\text{ur}}<0.5\)
9 \(G_0^{\text{ref}}\) small-strain shear modulus kPa HSsmall (optional)
10 \(\gamma_{0.7}\) shear strain at \(G/G_0=0.7\) – HSsmall (optional)
11 \(p_{\text{ref}}\) reference pressure kPa default \(100\)
12 \(K_w\) bulk modulus of water kPa \(0 \Rightarrow\) drained

Derived helpers kept by the code:

\[ M_\varphi=\frac{6\sin\varphi}{3-\sin\varphi},\qquad M_\psi=\frac{6\sin\psi}{3-\sin\psi}. \]

Driver options (opt_props)

These map to lib%material(i)%opt_props:

Idx Meaning Typical value
1 Integration scheme (1=Modified-Euler, 2=Euler–Richardson) 1.0
2 Stress error tolerance (dimensionless) 1.0e-4
3 Cap error tolerance (kept for parity) 1.0
4 Tangent assembly (1=end of sub-inc, 2=within sub-inc) 1.0
5 Numerical differentiation order (1 \(O(h)\), 2 \(O(h^2)\)) 2.0
6 Perturbation factor for finite differences 1.0e-6
7 Minimum mean effective stress \(p_{\min}\) (positive) 1.0e-2
8 Allow global cutback (1=yes, 0=no) 0
9 Minimum sub-increment size fraction of \(\Delta t\) 1.0e-5