Skip to content

Hypoplasticity + IGS

The IGS-Hypoplasticity model for sands extends the hypoplastic model initially proposed by Wolffersdorff (1996)1 by incorporating the intergranular strain concept (IGS) proposed by Niemunis and Herle, 19972. This enhancement, allows better simulation of the behaviour under unloading-reloading scenarios and cyclic loading, capturing features such as stiffness increase upon reversal loading, and reduced plastic strain rates under repeated loading cycles. The implementation details given in this section refer to the Hypo-IGS version of the two available hypoplastic models with IGS extension. A brief summary of the model and the implementation is given below, covering:


Constitutive Equations

General Formulation

The general hypoplastic equation is:

\[\label{eq:hypo-basic} \dot{\boldsymbol{\sigma}} = \mathsf{L}^{\text{hyp}} : \dot{\boldsymbol{\varepsilon}} + \boldsymbol{N}^{\text{hyp}} \|\dot{\boldsymbol{\varepsilon}}\| = \mathsf{M} : \dot{\boldsymbol{\varepsilon}} \]

Linear and nonlinear stiffness tensors:

\[ \mathsf{L}^{\text{hyp}} = f_b f_e \frac{1}{\hat{\boldsymbol{\sigma}} : \hat{\boldsymbol{\sigma}}}(F^2 \boldsymbol{I} + a^2 \hat{\boldsymbol{\sigma}} \hat{\boldsymbol{\sigma}}) \]
\[ \boldsymbol{N}^{\text{hyp}} = f_d f_b f_e \frac{F a}{\hat{\boldsymbol{\sigma}} : \hat{\boldsymbol{\sigma}}}(\hat{\boldsymbol{\sigma}} + \hat{\boldsymbol{\sigma}}^{\text{dev}}) \]

Therein \(\hat{\boldsymbol{\sigma}}=\dfrac{\boldsymbol{\sigma}}{\text{tr}\boldsymbol{\sigma}}\) and \(\hat{\boldsymbol{\sigma}}^{*}=\hat{\boldsymbol{\sigma}}-\frac{1}{3}\mathbf{I}\) are used. The scalar factors are defined by

\[ \begin{aligned} a=\frac{\sqrt{3}(3-\sin\varphi_c)}{2\sqrt{2}\sin\varphi_c}~,~ f_d=\left(\frac{e-e_d}{e_c-e_d}\right)^\alpha~,~f_e=\left(\frac{e_c}{e}\right)^\beta \end{aligned} \]

and

\[ \begin{aligned} f_{b}=\frac{h_s}{n}\left(\frac{e_{i0}}{e_{c0}}\right)^{\beta}\frac{1+e_i}{e_i}\left(\frac{3p}{h_s}\right)^{1-n}\left[3 +a^2 -a\sqrt{3}\left(\frac{e_{i0}-e_{d0}}{e_{c0}-e_{d0}}\right)^{\alpha}\right]^{-1}. \end{aligned} \]

The pressure-dependent void ratios \(e_i\), \(e_c\) and \(e_d\), describing the loosest, the critical and the densest state, are calculated using the relation proposed by Bauer (1996)3:

\[ \begin{aligned} \frac{e_i}{e_{i0}}=\frac{e_c}{e_{c0}}=\frac{e_d}{e_{d0}}=\exp\left[-\left(\frac{3p}{h_s} \right)^n \right]. \end{aligned} \]

\(p\) is the mean effective stress. The scalar factor \(F\) is given by

\[ \begin{aligned} F=\sqrt{\frac{1}{8}\tan(\psi)^2+\frac{2-\tan(\psi)^2}{2+\sqrt{2}\tan(\psi)\cos(3\theta)}} -\frac{1}{2\sqrt{2}}\tan(\psi), \end{aligned} \]

with \(\tan(\psi)=\sqrt{3}\|\hat{\boldsymbol{\sigma}}^*\|\) and

\[ \begin{aligned} \cos(3\theta)=-\sqrt{6}\frac{\text{tr}(\hat{\boldsymbol{\sigma}}^*\cdot\hat{\boldsymbol{\sigma}}^*\cdot\hat{\boldsymbol{\sigma}}^*)} {\left[\text{tr}(\hat{\boldsymbol{\sigma}}^*\cdot\hat{\boldsymbol{\sigma}}^*)\right]^{3/2}}. \end{aligned} \]

IGS Extension

To improve the performance of the hypoplastic model in the range of small strain cycles Niemunis & Herle(1997)2 introduced a new tensorial state variable, the intergranular strain \(\textbf{h}\), which memorizes the recent deformation history. With this extension the basic equation of the constitutive model reads:

\[ \dot{\boldsymbol{\sigma}} = \mathsf{M}(\boldsymbol{\sigma},\textbf{h},e) : \dot{\boldsymbol{\varepsilon}} \]

where \(\boldsymbol{\dot{\varepsilon}}\) is the strain rate. The stiffness tensor \(\mathsf{M}\) is defined by

\[ \mathsf{M}=\Big(\rho^{\chi}m_{T} + (1-\rho^{\chi})m_{R}\Big)\mathsf{L} \\ + \Bigg\{ \begin{array}{ll} \rho^{\chi}(1-m_{T})\mathbf{\mathsf{L}}:\vec{\mathbf {h}}\otimes\vec{\mathbf{h}} + \rho^{\chi}\mathbf{N}\vec{\mathbf{h}} \\ \rho^{\chi}(m_{R}-m_{T})\mathsf{L}:\vec{\mathbf{h}}\otimes\vec{\mathbf{h}} \end{array} ~~~\textrm{for}~~~ \begin{array}{ll}\vec{\mathbf{h}}:\dot{\boldsymbol{\varepsilon}}>0 \\ \vec{\mathbf{h}}:\dot{\boldsymbol{\varepsilon}}\le 0\end{array} \]

with the intergranular strain tensor \(\textbf{h}\), its degree of mobilization \(\rho\) and its direction \(\vec{\mathbf{h}}\) defined as

\[ \rho=\frac{\|\textbf{h}\|}{R}~~\textrm{and}~~\vec{\textbf{h}}=\frac{\mathbf{h}}{\|\textbf{h}\|}. \]

\(\chi\), \(m_T\), \(m_R\) and \(R\) are material parameters controlling the influence of the intergranular strain. The evolution law of the intergranular strain \(\textbf{h}\) is

\[ \dot{\textbf{h}} = \Bigg\{ \begin{array}{ll} (\mathbf{\sf{I}}-\vec{\textbf{h}}\otimes\vec{\textbf{h}}\rho^{\beta_r}) : \boldsymbol{\dot{\varepsilon}}\\ \boldsymbol{\dot{\varepsilon}} \end{array} ~~~\textrm{for}~~~ \begin{array}{ll}\vec{\mathbf{h}}:\boldsymbol{\dot{\varepsilon}}>0 \\ \vec{\mathbf{h}}:\boldsymbol{\dot{\varepsilon}}\le 0\end{array} \]

where \(\beta_R\) is another material parameter and \(\mathbf{\sf{I}}\) is the fourth-order identity tensor.

Model Parameters

Symbol Description Unit
\(\varphi_c\) Critical state friction angle degrees (°)
\(h_s\) Granular hardness kPa
\(n_B\) Barotropy exponent
\(e_{d0}\) Minimum void ratio at p = 0
\(e_{c0}\) Critical void ratio at p = 0
\(e_{i0}\) Maximum void ratio at p = 0
\(a\) Dilatancy exponent
\(b\) Pyknotropy exponent
\(m_R\) Stiffness factor for elastic response for changes of the direction of the strain path by \(180^\circ\)
\(m_T\) Stiffness factor for elastic response for changes of the direction of the strain path by \(90^\circ\)
\(R\) Elastic strain amplitude threshold
\(\beta_R\) IS hardening parameter
\(\chi\) IS exponent

These parameters are calibrated based on experimental data, commonly including drained/undrained triaxial tests and shear modulus degradation curves.


Numerical Integration of the Constitutive Model

The constitutive update of the IGS-hypoplastic model is performed via explicit time integration. Two integration schemes are available:

  • Forward Euler integration with fixed subincrements
  • Euler–Richardson integration with adaptive stepsize control

For a discussion of the integration schemes it is referred to the documentation of the Hypo-ISA section of the Theory Manual.


Optional features

numgeo implements a couple of optional features to improve numerical stability of the implementation. The use of these features is not mandatory and can be adjusted by the user.

Projection onto the Matsuoka–Nakai surface

If the trial stress state exceeds an auxiliary Matsuoka–Nakai (MN) failure surface, the stress is projected radially within the deviatoric plane back onto this surface, while maintaining a constant mean effective stress \(p=-\mbox{tr}(\boldsymbol{T})/3\), thus maintaining consistency with the volumetric behaviour prescribed by the hypoplastic formulation. The MN surface is defined by a cut-off friction angle \(\varphi^{\text{MN}}\), which is selected in a state-dependent manner to reflect the evolving strength conditions. More details are provided in the Hypo-ISA section of the Theory Manual or in 4.

In the current implementation of numgeo, users can select between different strategies for stress projection onto the Matsuoka–Nakai (MN) failure surface via the phi_cut parameter:

  • projection_dev <= 0:

    Deactivates the stress projection (default). The stress state is not constrained by the MN surface.

  • projection_dev = 1:

    Activates the projection onto the MN failure surface using a variable friction angle defined as \(\varphi^{\text{MN}} = \max\left(\varphi_{\text{mob}}^{\text{hyp}}, ~\tan^{-1}(f_{e0} \tan \varphi_c)\right)\) , where \(\varphi_c\) is the critical state friction angle, \(f_d\) and \(f_e\) are hypoplastic factors depending on void ratio and relative density.

  • projection_dev = 2:

    Activates the projection onto the MN failure surface using a variable friction angle defined as \(\varphi^{\text{MN}} = \max\left(\varphi_{\text{mob}}^{\text{hyp}}, ~\tan^{-1}(f_{e0} f_{d0} \tan \varphi_c)\right)\) , where \(\varphi_c\) is the critical state friction angle, \(f_d\) and \(f_e\) are hypoplastic factors depending on void ratio and relative density.

  • projection_dev > 2:

    Activates the projection onto the MN surface using a fixed user-defined cut-off friction angle \(\varphi^{\text{MN}} = \varphi^{\text{cut}}\) , with \(\varphi^{\text{cut}}\) given in degrees.


Enforcement of minimum mean effective stress

To prevent unphysical stress states with excessively low or tensile mean effective stress, the implementation includes an optional stress projection that enforces a user-defined lower bound on the mean effective stress, denoted as \(p_{\text{min}}\) (compression is positive).

This is particularly relevant for numerical robustness under undrained conditions, where liquefaction or large pore water pressure accumulation may cause the effective stress to drop to zero or become tensile.

After each constitutive update, the stress tensor \(\boldsymbol{\sigma}\) is checked for compliance with the pressure condition:

\[ p = -\frac{1}{3} \text{tr}(\boldsymbol{\sigma}) = -\frac{1}{3} (\sigma_{11} + \sigma_{22} + \sigma_{33}) \]

If the current mean effective stress satisfies:

\[ p \geq p_{\text{min}}, \]

no modification is applied. Otherwise, the projection step adjusts the normal stress components:

\[ \sigma_{ii} \leftarrow \sigma_{ii} + (p - p_{\text{min}}) \quad \text{for } i = 1,2,3 \]

This operation shifts the stress state isotropically to ensure:

\[ p = p_{\text{min}} \]

The deviatoric stress state remains unchanged.

  • The projection is only applied if the computed mean effective stress violates the lower bound.
  • The user defines the value of \(p_{\text{min}}\), per default \(p_{\text{min}}=10^{-2}\) F/A.

Comparison to the soilmodels.com reference implementation

This section outlines key differences between the ISA-hypoplastic implementation presented here (e.g. in numgeo) and the freely available implementation provided on soilmodels.com.

Phantom Elasticity (PE)

The implementation on soilmodels.com includes a continuously active elastic component known as Phantom Elasticity (PE), which contributes additively to the stress rate throughout the entire calculation. This additional stress increment is given by:

\[ \dot{\boldsymbol{\sigma}} = \dot{\boldsymbol{\sigma}}^{\text{ISA}} + \dot{\boldsymbol{\sigma}}^{\text{PE}} \quad \text{with} \quad \dot{\boldsymbol{\sigma}}^{\text{PE}} = \mathsf{C} : \dot{\boldsymbol{\varepsilon}}, \]

where:

  • \(\dot{\boldsymbol{\sigma}}^{\text{ISA}}\) is the rate from the ISA-hypoplasticity model,
  • \(\dot{\boldsymbol{\sigma}}^{\text{PE}}\) is the elastic stress rate,
  • \(\mathsf{C}\) is a linear elastic stiffness tensor defined by a phantom Young's modulus \(E^{ph}\) and Poisson's ratio \(\nu^{ph}\).

In the soilmodels.com version (as of 22.05.2025), the default parameters are:

  • Young’s modulus: \(E^{ph} = 20\) kPa
  • Poisson’s ratio: \(\nu^{ph} = 0.45\)

numgeo

In contrast, the ISA model implemented in numgeo provides Phantom Elasticity as an optional feature. By default, PE is disabled (i.e. \(E = 0\)). Users can enable PE by specifying the elastic parameters using the keyword *Optional mechanical parameter. Further details on this keyword and its usage can be found in the Reference Manual. Furthermore, deviating from the original implementation from W. Fuentes, the phantom elasticity in the numgeo implementation is only active when the mean effective stress is equal to or drops below \(p_{min,ph}\).


References


  1. P.-A. Wolffersdorff. A hypoplastic relation for granular materials with a predefined limit state surface. Mechanics of Cohesive-frictional Materials, 1(3):251–271, 1996. doi:10.1002/(SICI)1099-1484(199607)1:3<251::AID-CFM13>3.0.CO;2-3

  2. A. Niemunis and I. Herle. Hypoplastic model for cohesionless soils with elastic strain range. Mechanics of Cohesive-frictional Materials, 2(4):279–299, 1997. doi:10.1002/(SICI)1099-1484(199710)2:4<279::AID-CFM29>3.0.CO;2-8

  3. Erich Bauer. Calibration of a Comprehensive Hypoplastic Model for Granular Materials. Soils and Foundations, 36(1):13–26, 1996. URL: https://linkinghub.elsevier.com/retrieve/pii/S0038080620313433, doi:10.3208/sandf.36.13

  4. Jan Machaček, Merita Tafili, and Carlos Eduardo Grandas Tavera. Overshooting of hypoplastic models with in- tergranular strain: manifestation, implications and remedy. International Journal for Numerical and Analytical Methods in Geomechanics, 2025-08.