Hypoplasticity + ISA

The ISA-Hypoplasticity model for sands extends the hypoplastic model initially proposed by Wolffersdorff (1996)¹ by incorporating intergranular strain anisotropy (ISA). This enhancement, proposed by Fuentes and Triantafyllidis (2015)², allows better simulation of the behaviour under cyclic loading, capturing features such as threshold strain amplitude, stiffness increase upon reversal loading, and reduced plastic strain rates under repeated loading cycles. The ISA-hypoplastic model was further extended by Poblete et al. (2016)³ and further improved by Fuentes et al. (2020)⁴ to account explicitly for cyclic mobility effects, enabling simulations of liquefaction phenomena, which occur under undrained cyclic conditions and large strain amplitudes. A brief summary of the model and the implementation is given below, covering:

• Constitutive Equations
  • General Formulation
  • ISA Extension
  • Cyclic Mobility Extension
• Model parameters
• Numerical Integration of the Constitutive Model
  • Forward Euler
  • Euler-Richardson
• Optional features
  • Projection onto the Matsuoka–Nakai surface
  • Enforcement of minimum mean effective stress
• Comparison to the soilmodels.com reference implementation
  • Phantom Elasticity (PE)
  • Evolution Equation for \(\beta_h\)

---

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

ISA modification of \(f_d\) and \(f_e\)

The factors \(f_d\) and \(f_e\) are modified in the ISA formulation to account for cyclic mobility effects, see Section ISA Extension.

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

\[ \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} \]

ISA Extension

The ISA²\(^,\)³\(^,\)⁴ approach, understood as an extension for hypoplastic models to enable simulations on cyclic loading, is based on the concept of intergranular strain (IS), originally introduced by Niemunis and Herle (1997)⁶ in a former formulation. The main difference of the ISA formulation with respect to conventional IS, is the introduction of an elastic locus to describe the small strain behavior. According to the ISA theory, the rate of the intergranular strain \(\boldsymbol{h}\) evolves through an elastoplastic formulation:

\[ \dot{\boldsymbol{h}} = \begin{cases} \dot{\boldsymbol{\varepsilon}} & \text{for } F_H < 0 ~~~\text{(elastic)}\\ \dot{\boldsymbol{\varepsilon}} - \dot{\lambda}_H \boldsymbol{N} & \text{for } F_H = 0 ~~~\text{(plastic)} \end{cases} \]

with

\[ \dot{\lambda}_H = \frac{\langle \boldsymbol{N} : \dot{\boldsymbol{\varepsilon}} \rangle}{1 + \boldsymbol{N} : \bar{\boldsymbol{c}}}, \quad \boldsymbol{N} = \frac{\boldsymbol{h} - \boldsymbol{c}}{\|\boldsymbol{h} - \boldsymbol{c}\|} \]

where \(F_H\) is a yield surface function, \(\dot{\lambda}\) is the plastic multiplier and \(\boldsymbol{N}=\partial F_H/\partial \boldsymbol{h}\) is an associated flow rule. The yield function \(F_H\) was formulated to consider an elastic threshold strain \(\|\Delta \boldsymbol{\varepsilon}\|=R\). Considering that under elastic conditions \(F_H < 0\), increments in strains are equal to increments in intergranular strain, i.e. \(\Delta \boldsymbol{\varepsilon} = \Delta \boldsymbol{h}\), the elastic strain amplitude is simply described through the following yield surface function \(F_H\):

\[ F_H = \| \boldsymbol{h} - \boldsymbol{c} \| - \frac{R}{2} \]

where tensor \(\boldsymbol{c}\) is the intergranular back strain and defines the center of the yield surface (kinematic hardening variable). The intergranular back strain evolves according to

\[ \dot{\boldsymbol{c}} = \dot{\lambda}_H \bar{\boldsymbol{c}}, \quad \bar{\boldsymbol{c}} = \frac{b_h(\boldsymbol{c}_b - \boldsymbol{c})}{R}, \quad \boldsymbol{c}_b = \frac{R}{2}\frac{\dot{\boldsymbol{\varepsilon}}}{\|\dot{\boldsymbol{\varepsilon}}\|} \]

To account for small strain effects, the ISA approach reformulates the continuum tangent stiffness \(\mathsf{M}\) from Equation (\(\ref{eq:hypo-basic}\)) according to the following relation:

\[ \mathsf{M} = \begin{cases} m(\mathsf{L}^{\text{hyp}} + \rho^\chi \boldsymbol{N}^{\text{hyp}} \boldsymbol{N}) & \text{for } F_H = 0 \\ m_R \mathsf{L}^{\text{hyp}} & \text{for } F_H < 0 \end{cases} \]

Therein, \(m\) and \(\rho\) are factors controlling stiffness reduction and plastic rate, respectively:

\[ m = m_R + (1 - m_R)y_h, \quad y_h = \rho^\chi \langle \boldsymbol{N} : \dot{\boldsymbol{\varepsilon}} \rangle \]

\[ \rho^\chi = 1 - \frac{\|\boldsymbol{d}^b\|}{2R}, \quad \boldsymbol{d}^b = \boldsymbol{h}_b - \boldsymbol{h}, \quad \boldsymbol{h}_b = R\boldsymbol{N} \]

Accumulated strain evolution:

\[ \dot{\varepsilon}_{\text{acc}} = \frac{\epsilon^{acc}}{R}(1 - y_h - \varepsilon_{\text{acc}})\|\dot{\boldsymbol{\varepsilon}}\| \]

Exponent \(\chi\) controlling plastic strain rate:

\[ \chi = \chi_0 + \varepsilon_{\text{acc}}(\chi_{\text{max}} - \chi_0) \]

Cyclic Mobility Extension

Additional state variable \(\boldsymbol{Z}\):

\[ \dot{\boldsymbol{Z}} = c_z \langle F_d \rangle (\boldsymbol{N} - \boldsymbol{Z}) \|\dot{\boldsymbol{\varepsilon}}\| \]

Function \(F_d\) detecting dilatant behavior:

\[ F_d = \frac{q/p}{M_c F f_{d0}} - 1 \]

Factor \(f_z\) controlling cyclic mobility:

\[ f_z = \langle -\boldsymbol{Z}:\boldsymbol{N} \rangle \]

Modified factors \(f_d\) and \(f_e\) for cyclic mobility paths:

\[ f_e = f_{e0} - f_z \langle f_{e0}-1 \rangle, \quad f_{e0} = \left(\frac{e_c}{e}\right)^\beta \]

\[ f_d = f_{d0} + f_z \langle 1 - f_{d0} \rangle, \quad f_{d0} = \left(\frac{e - e_d}{e_c - e_d}\right)^\alpha \]

Factor \(b_h\) adaptation:

\[ \beta_h = \beta_{hmax} + (\beta_{h0}-\beta_{hmax})(1-f_z)f_h, \quad f_h = |\dot{\boldsymbol{\varepsilon}} : \boldsymbol{d}^b| \]

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	–
\(R\)	Elastic strain amplitude threshold	–
\(\beta_{h0}\)	Minimum IS hardening parameter	–
\(\beta_{h\text{max}}\)	Maximum IS hardening parameter	–
\(\chi_0\)	Minimum IS exponent	–
\(\chi_{\text{max}}\)	Maximum IS exponent	–
\(\epsilon^{acc}\)	Accumulation rate factor	–
\(c_z\)	Cyclic mobility factor	–

These parameters are calibrated based on experimental data, commonly including drained/undrained triaxial tests and shear modulus degradation curves. Refer to Fuentes et al. (2020)⁴ for details on calibration strategies.

---

Numerical Integration of the Constitutive Model

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

Forward Euler integration with fixed subincrements

A simple and robust approach is the use of the forward Euler method with constant subincrement size. Given a total strain increment \(\Delta \boldsymbol{\varepsilon}\), the increment is divided into a predefined number of substeps, and the stress is updated incrementally.

Subincrement Size Determination The number of subincrements \(n_{\text{sub}}\) is computed based on the norm of the strain increment:

\[ n_{\text{sub}} = \min\left( \max\left( \left\lfloor \frac{\| \Delta \boldsymbol{\varepsilon} \|}{\varepsilon_{\text{tol}}} \right\rfloor, 1 \right), N_{\text{max}} \right) \]

Where:

• \(| \Delta \boldsymbol{\varepsilon} |\) is the Euclidean norm of the strain increment,
• \(\varepsilon_{\text{tol}} = 10^{-6}~\) is the strain subincrement tolerance,
• \(N_{\text{max}} = 50000\) is the maximum allowed number of subincrements.

Each subincrement is then integrated using the basic forward Euler update:

\[ \boldsymbol{\sigma}_{n+1} = \boldsymbol{\sigma}_n + \mathsf{M}(\boldsymbol{\sigma}_n, \boldsymbol{\varepsilon}_n) : \Delta \boldsymbol{\varepsilon}_{\text{sub}} \]

where \(\mathsf{M}\) is the tangent operator evaluated at the current state. The Forward Euler integration scheme is straightforward to implement, robust, and computationally inexpensive. However, it offers limited control over numerical accuracy. In particular, the integration result can be sensitive to the size of the applied strain increment.

Figure 1 illustrates this effect by showing the simulation results of an undrained cyclic triaxial test using three different strain increment sizes \(\Delta \varepsilon_2\).

Figure 1: Simulation of a cyclic undrained triaxial test - influence of strain increment size on the result of Forward Euler integration.

Reducing the strain increment size generally increases the accuracy of the model integration. The deviations in the stress response observed in Figure 1 clearly demonstrate the importance of choosing an appropriate increment size when using the Forward Euler method.

However, due to the use of a constant step size, smaller increments inevitably lead to longer simulation times. In the present case, reducing the strain increment from \(\Delta \varepsilon_2 = 10^{-5}\) to \(10^{-6}\) increases the runtime from approximately 3 seconds to 13 seconds. A further reduction to \(\Delta \varepsilon_2 = 10^{-7}\) results in a total runtime of about 120 seconds. While the absolute values depend on the hardware and implementation, the trend is consistent and highlights the computational trade-off between accuracy and efficiency.

Euler–Richardson integration with adaptive stepsize control

An enhanced integration scheme is available using the Euler–Richardson method with adaptive step size control to maintain numerical accuracy and efficiency as outlined in Fellin et al (2009)⁷.

The method follows these steps per subincrement:

• Compute a provisional stress update \(\boldsymbol{\sigma}_1\) using the forward Euler step.
• Compute the midpoint strain and use it to estimate a midpoint stress rate, resulting in \(\boldsymbol{\sigma}_2\).
• Estimate the local integration error as:

\[ \text{TOL}^\sigma = \frac{\| \boldsymbol{\sigma}_2 - \boldsymbol{\sigma}_1 \|}{\min\left(1.0, \; 2.0 \cdot \| \boldsymbol{\sigma}_0 \| \right)} \]

Where \(\boldsymbol{\sigma}_0\) is the initial stress in the current increment.

If \(\text{TOL}^\sigma > 10^{-4}\), the step is rejected and a smaller step is attempted. Otherwise, the step is accepted and the algorithm proceeds with the next strain portion.

The Euler–Richardson integration scheme offers a more flexible and efficient approach by adaptively adjusting the step size based on the local nonlinearity of the response. It permits larger steps in regions of smooth behaviour and automatically reduces the step size when high accuracy is required, such as near stress reversals or phase transitions.

Figure 2 compares the performance of the Euler–Richardson method with the conventional Forward Euler scheme. In this example, the Euler–Richardson scheme uses a maximum strain increment of \(\Delta \varepsilon_2 = 10^{-5}\) and targets a maximum allowable integration error of \(\text{TOL}^\sigma < 10^{-4}\), defined in terms of stress. As a reference solution, the Forward Euler scheme is used with a much smaller subincrement size of \(\Delta \varepsilon_2 = 10^{-7}\) to ensure high accuracy.

The results demonstrate that the Euler–Richardson scheme achieves comparable accuracy while requiring significantly fewer substeps. The total simulation time for the adaptive scheme was approximately 2.5 seconds—nearly 50 times faster than the reference Forward Euler solution, which required around 120 seconds. This highlights the potential of adaptive integration strategies to significantly improve computational efficiency without compromising accuracy.

Figure 2: Simulation of a cyclic undrained triaxial test - comparison of the Euler–Richardson scheme (\(\Delta \varepsilon_2 = 10^{-5}\), \(\text{TOL}^\sigma < 10^{-4}\)) with the Forward Euler scheme (\(\Delta \varepsilon_2 = 10^{-7}\)).

---

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

The MN surface is expressed in terms of the stress tensor \(\boldsymbol{\sigma}\) through a cubic condition on the deviatoric stress tensor \(\boldsymbol{\sigma}^* = \boldsymbol{\sigma} + p \boldsymbol{\delta}\):

\[ \det(\boldsymbol{\sigma}^*) = \gamma p^3, \]

where the scalar parameter \(\gamma\) depends on the cut-off friction angle \(\varphi^{\text{MN}}\) and is given by:

\[ \gamma = \frac{1 - \sin^2 \varphi^{\text{MN}}}{9 - \sin^2 \varphi^{\text{MN}}}. \]

To capture the nonlinear evolution of the strength envelope of the base hypoplastic model, a dynamic definition of \(\varphi^{\text{MN}}\) based on the friction angle mobilised by a pseudo-hypoplastic stress state \(\boldsymbol{\sigma}^{\text{hyp}}\), computed in each sub-increment of the integration algorithm is proposed herein. This stress state is obtained by evaluating a modified hypoplastic tangent stiffness \(\mathsf{M}^{\text{hyp}}\) with the intergranular strain terms deactivated (\(m_R = 1\), \(\rho = 1\)), thereby yielding a purely hypoplastic response:

\[ \boldsymbol{\sigma}^{\text{hyp}} = \boldsymbol{\sigma}_0 + \mathsf{M}^{\text{hyp}} : \Delta \boldsymbol{\varepsilon}, \]

where \(\boldsymbol{\sigma}_0\) is the stress at the beginning of the current sub-increment. From \(\boldsymbol{\sigma}^{\text{hyp}}\), the mobilised friction angle \(\varphi_{\text{mob}}^{\text{hyp}}\) is computed and used to define the MN cut-off friction angle:

\[ \varphi^{\text{MN}} = \varphi_{\text{mob}}^{\text{hyp}}. \]

The mobilised friction angle is defined as

\[ \varphi_{\text{mob}} = \arcsin\left( \frac{3 \eta}{6 + \eta \cos(3\theta)} \right), \]

where \(\eta = q/p\) is the stress ratio, and \(\theta\) is the Lode angle.

Formulating \(\varphi^{\text{MN}}\) depending on \(\varphi_{\text{mob}}\) ensures that the deviatoric projection surface evolves consistently with the material state and loading history. It captures both the initial stiffness and the peak strength development, while avoiding artificial saturation or abrupt yield behaviour. Importantly, the projection is only activated if the stress predicted by the Hypo-IGS/ISA model exceeds the strength envelope implied by the base hypoplastic model. In such cases, the stress increment is regularised by anchoring it to the strength level that the base model would have predicted under identical strain input. This provides physically meaningful and numerically efficient regularisation mechanism for the Hypo-ISA formulation, requiring only minor modifications to existing explicit integration schemes.

To avoid underestimating shear strength in dense states or at the onset of loading, where \(\varphi_{\text{mob}}^{\text{hyp}}\) may be low, a lower bound is enforced:

\[ \varphi^{\text{MN}} = \max\left(\varphi_{\text{mob}}^{\text{hyp}}, ~\tan^{-1}(f_{e0} f_{d0} \tan \varphi_c)\right). \]

This guarantees that the projection surface does not fall below the minimum strength level predicted by the base hypoplastic model, thereby preserving a realistic stress response in early or small-strain loading stages. Further details on the numerical implementation of this stress projection procedure are provided in Machacek et al. (2025)⁸.

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.

Examples showcasing the effectiveness of the proposed stress projection methods are provided in Machacek et al. (2025)⁸.

---

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

Evolution Equation for \(\beta_h\)

Another key deviation lies in the evolution of the intergranular strain hardening parameter \(\beta_h\). The numgeo implementation adheres to the equation presented in the original ISA-hypoplasticity publication:

\[ \beta_h = \beta_{h\text{max}} + (\beta_{h0} - \beta_{h\text{max}})(1 - f_z)f_h \]

However, the soilmodels.com version modifies this expression. Instead, the following evolution rule is used:

\[ \beta_h = \beta_{h\text{max}} + (\beta_h - \beta_{h\text{max}})(1 - f_z)f_h \]

The key difference is that the soilmodels implementation evolves \(\beta_h\) recursively from its current value, whereas numgeo uses a fixed reference value \(\beta_{h0}\) as the lower bound. This may lead to differences in cyclic behaviour, especially in how stiffness degradation and recovery are modelled across load reversals.

---

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. W. Fuentes and Th. Triantafyllidis. ISA model: A constitutive model for soils with yield surface in the intergranular strain space: ISA MODEL: A CONSTITUTIVE MODEL FOR SOILS. International Journal for Numerical and Analytical Methods in Geomechanics, 39(11):1235–1254, 2015. URL: https://onlinelibrary.wiley.com/doi/10.1002/nag.2370, doi:10.1002/nag.2370. ↩↩

3. M. Poblete, W. Fuentes, and Th. Triantafyllidis. On the simulation of multidimensional cyclic loading with intergranular strain. Acta Geotechnica, 11(6):1263–1285, 2016. URL: http://link.springer.com/10.1007/s11440-016-0492-2, doi:10.1007/s11440-016-0492-2. ↩↩

4. William Fuentes, Torsten Wichtmann, Melany Gil, and Carlos Lascarro. ISA-Hypoplasticity accounting for cyclic mobility effects for liquefaction analysis. Acta Geotechnica, 15(6):1513–1531, 2016. URL: http://link.springer.com/10.1007/s11440-019-00846-2, doi:10.1007/s11440-019-00846-2. ↩↩↩

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

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

7. Wolfgang Fellin, Markus Mittendorfer, and Alexander Ostermann. Adaptive integration of constitutive rate equations. Computers and Geotechnics, 36(5):698–708, 2009. URL: https://linkinghub.elsevier.com/retrieve/pii/S0266352X08001547, doi:10.1016/j.compgeo.2008.11.006. ↩

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