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 accoring 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}} + q^v \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

In hypoplasticity, the failure and bounding surfaces are geometrically similar. However, studies have shown that stress paths can lead to states lying outside the failure surface ⁸\(^,\)⁹. This may result in unrealistically high mobilised friction angles, reaching values up to 70° ⁸ - which are not just questionable from a soil mechanics standpoint, but can also introduce severe numerical issues. Specifically, such unbounded states are typically associated with an ill-conditioned material Jacobian, which can degrade global convergence behaviour in finite element simulations.

To mitigate these issues, numgeo implements an optional stress projection onto the Matsuoka–Nakai (MN) failure surface based on the work of ⁹. If activated, stress states that violate an imaginary failure condition defined by a cut-off friction angle \(\varphi^{MN}\) are projected back onto the MN failure surface. This projection is performed radially in the deviatoric plane while preserving the mean effective stress \(p\).

The MN surface in terms of the stress tensor \(\boldsymbol{\sigma}\) is defined by a cubic function in the deviatoric plane:

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

where \(\boldsymbol{s} = \boldsymbol{\sigma} + p \boldsymbol{I}\) is the deviatoric stress and \(\gamma\) is a function of the friction angle \(\varphi^{MN}\):

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

Details on the algorithm for the stress projection are provided in the expendable box below.

Projection algorithm

1. Compute mean stress and deviatoric stress:

   \[ p = -\frac{1}{3} \text{tr}(\boldsymbol{T}) \]
   \[ \boldsymbol{T}_{\text{dev}} = \boldsymbol{T} + p \boldsymbol{I} \]

2. Normalize the Deviatoric Tensor:

   \[ \boldsymbol{Z} = \frac{\boldsymbol{T}_{\text{dev}}}{\text{tr}(\boldsymbol{T})}, \quad Z_{ii} \leftarrow Z_{ii} - \frac{1}{3} \]
   \[ \|\boldsymbol{Z}\| = \sqrt{Z_{ij} Z_{ij}}, \quad \boldsymbol{Z} \leftarrow \frac{\boldsymbol{Z}}{\|\boldsymbol{Z}\|} \]

3. Compute Determinant and Solve Cubic Equation:

   \[ \text{Let } \delta = \det(\boldsymbol{Z}) \]

   Define:

   \[ r = \frac{\gamma}{2} - \frac{1}{6}, \quad tmp = \frac{1}{27} - \frac{\gamma}{3} \]

   • If \(\delta \approx 0\), use the analytic approximation:

     \[ r = \sqrt{\max\left(0, -\frac{tmp}{r}\right)} \]

   • Otherwise, solve the cubic equation using Cardano's formula:

     Define:

     \[ a = -\frac{r^2}{3}, \quad b = \frac{2r^3}{27} + \frac{tmp}{\delta}, \quad R = \sqrt{\left|\frac{a}{3}\right|}, \quad \omega = \arccos\left (\frac{b}{2 R^3}\right) \]

     Then the real roots are:

     \[ r_i = -2 R \cos\left(\frac{\omega + 2\pi i}{3}\right) - \frac{r}{3}, \quad i = 0,1,2 \]

     Choose the largest admissible positive root.

4. Apply Radial Return Mapping:

   \[ \boldsymbol{T}_{\text{proj}} = -\text{tr}(\boldsymbol{T}) \cdot (r \boldsymbol{Z} + \frac{1}{3} \boldsymbol{I}) \]

If the computed radius is smaller than the current normalized radius or zero, the original stress is left unchanged.

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:

• phi_cut = -1:

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

• phi_cut = 0:

  Activates the projection onto the MN failure surface using a variable friction angle defined as \(\varphi^{\text{MN}} = \tan^{-1}\left( f_e \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.

• phi_cut > 0:

  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.

The effectiveness of the implemented stress projection methods is demonstrated by reproducing an over-shooting effect reported by Stastny et al. (2025)¹⁰.

The test setup involves a drained monotonic triaxial compression test on an initially dense sand specimen (\(D_{r0} = 0.95\)). The monotonic shearing is interrupted after 0.5 % axial strain by 10 consecutive cycles of un- and reloading. After this cyclic phase, monotonic shearing is resumed until an axial strain of 15 % is reached. Stastny et al. (2025)¹⁰ observed a pronounced overshooting in deviatoric stress during post-cyclic monotonic loading when no stress projection was applied. This phenomenon is replicated in the numgeo simulations without projection.

Figure 3 shows simulation results for all three projection options. The overshooting is effectively suppressed when either of the projection schemes is applied. A reference simulation without any cyclic interruption is also included for comparison.

Figure 3: Simulation of the overshooting effect in drained triaxial compression following cyclic loading on initially dense sand. Comparison of different stress projection strategies implemented in numgeo. Results aligned with the findings of Stastny 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.

Here is a new section for your documentation, clearly highlighting the deviations from the reference implementation available on soilmodels.com. The section integrates the equations from your uploaded images and maintains consistency with the rest of your mkdocs-material structure.

---

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. URL: https://onlinelibrary.wiley.com/doi/abs/10.1002/%28SICI%291099-1484%28199607%291%3A3%3C251%3A%3AAID-CFM13%3E3.0.CO%3B2-3, 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. URL: https://onlinelibrary.wiley.com/doi/abs/10.1002/%28SICI%291099-1484%28199710%292%3A4%3C279%3A%3AAID-CFM29%3E3.0.CO%3B2-8, 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. Shun Wang, Wei Wu, Chong Peng, Xuzhen He, and Deshan Cui. Numerical integration and FE implementation of a hypoplastic constitutive model. Acta Geotechnica, 13(6):1265–1281, 12 2018. URL: http://link.springer.com/10.1007/s11440-018-0684-z, doi:10.1007/s11440-018-0684-z. ↩↩

9. S. Chrisopoulos, V. A. Osinov, and T. Triantafyllidis. Dynamic Problem for the Deformation of Saturated Soil in the Vicinity of a Vibrating Pile Toe. In Theodoros Triantafyllidis, editor, Holistic Simulation of Geotechnical Installation Processes: Benchmarks and Simulations, Lecture Notes in Applied and Computational Mechanics, pages 53–67. Springer International Publishing, 2016. URL: https://doi.org/10.1007/978-3-319-23159-4_3, doi:10.1007/978-3-319-23159-4_3. ↩↩

10. A. Stastny, G. Medicus, V. Galavi, M. Tafili, and F. Tschuchnigg. Evaluation of advanced soil models for the cyclic soil-structure interaction of integral bridges. Acta Geotechnica, 2025. URL: https://link.springer.com/10.1007/s11440-025-02632-9, doi:10.1007/s11440-025-02632-9. ↩↩