Skip to content

Ta-Ger sand model

The Ta-Ger sand model is an elastoplastic constitutive model for sand of Tasiopoulou and Gerolymos 12 that combines a single-surface bounding-surface concept with Bouc-Wen type smooth hysteresis. The model has no finite purely elastic domain: the transition from elastic to perfectly plastic response is controlled by the mapping variable \(\zeta\) and the hardening exponent \(n\).

The implementation in numgeo follows the Ta-Ger formulation with the progressive calibration concept based on Bolton's relative dilatancy index. The model is intended for monotonic and cyclic loading of sands and includes pressure-dependent stiffness, Lode-angle dependency, Rowe-type stress-dilatancy and an intermediate-principal-stress correction for the peak friction angle.

Availability

The Ta-Ger model will be available with the upcoming release of numgeo

Contents

Overview

The model consists of the following components:

  1. Stress-dependent elasticity: the shear and bulk moduli are functions of the current mean effective stress. Two shear-modulus formulations are available in the implementation.
  2. Explicit elastoplastic stiffness: the constitutive matrix is formulated directly in stress-increment-strain-increment form. No plastic modulus and no return-mapping loading index are required.
  3. Bouc-Wen type mapping: the diagonal plastic matrix is scaled by \(\zeta^n\), where \(\zeta\) measures the distance of the current stress state from the bounding surface and \(n\) controls the smoothness of the transition to failure.
  4. Critical-state-compatible stress-ratio evolution: the bounding stress ratio \(M_s\) and the phase transformation stress ratio \(M_{pt}\) evolve towards the critical state with the accumulated deviatoric strain measure.
  5. Rowe-type dilatancy: the plastic potential gradient uses a scalar dilatancy ratio that is distributed anisotropically to the normal plastic strain increments.
  6. Stress-induced peak-strength correction: the factor \(g_{\alpha b}\) is computed from the intermediate principal stress parameter \(b\) and used in the peak friction angle relation.

Sign Convention and Stress Measures

The stress vector stress follows the solver convention, where compression is negative. Internally, several stress-ratio calculations are carried out in compression-positive form. The Roscoe mean stress is

\[ \begin{equation} p = -\frac{\sigma_{11}+\sigma_{22}+\sigma_{33}}{3} . \end{equation} \]

The deviatoric stress measure \(q\) is computed from the second deviatoric invariant \(J_2\) as

\[ \begin{equation} q = \sqrt{3J_2} . \end{equation} \]

A small compressive stress cut-off is applied before evaluating \(p\) and \(q\) to prevent invalid logarithms and singular stress ratios close to the tensile or hydrostatic limit.

The current void ratio \(e\) evolves with volumetric strain according to

\[ \begin{equation} e_{n+1} = e_n + \Delta \varepsilon_v (1+e_n) , \end{equation} \]

where the sign follows the strain convention used in the implementation.

Elasticity

The implementation supports two shear-modulus formulations, selected by the material parameter Gmodel.

For Gmodel = 1, the shear modulus follows the calibration-oriented form 2

\[ \begin{equation} G = G_0\,(0.13D_{r0}+3.6)\,p^m , \end{equation} \]

where \(D_{r0}\) is the initial relative density. This option follows the stiffness form used in the progressive calibration procedure.

For Gmodel = 2, the pressure- and void-ratio-dependent form 1 is

\[ \begin{equation} G = G_0 p_a \frac{(2.97-e)^2}{1+e}\left(\frac{p}{p_a}\right)^m , \end{equation} \]

with \(p_a = 100\) in the implementation. The bulk modulus \(K\) is obtained from \(G\) and Poisson's ratio \(\nu\):

\[ \begin{equation} K = \frac{2G(1+\nu)}{3(1-2\nu)} . \end{equation} \]

The isotropic elastic matrix is assembled from \(G\) and \(K\).

Elastoplastic Matrix

The Ta-Ger model uses an explicit elastoplastic matrix of the form

\[ \begin{equation} \mathbf{E}^{ep}=\mathbf{E}^{e}\left[\mathbf{I}- \frac{\boldsymbol{\Phi}_g\boldsymbol{\Phi}_f^T\mathbf{E}^{e}\mathbf{H}} {\boldsymbol{\Phi}_f^T\mathbf{E}^{e}\boldsymbol{\Phi}_g}\right] , \end{equation} \]

where \(\mathbf{E}^{e}\) is the elastic matrix, \(\boldsymbol{\Phi}_f\) is the gradient of the loading surface and \(\boldsymbol{\Phi}_g\) is the gradient of the plastic potential. The plastic matrix is diagonal and defined as

\[ \begin{equation} \mathbf{H} = \zeta^n \mathbf{I} . \end{equation} \]

The parameter \(\zeta\) controls the position of the loading surface relative to the bounding surface. The hardening exponent \(n\) controls the smoothness of the transition from near-elastic response to the perfectly plastic limit.

In the implementation, the stress increment in a substep is computed from

\[ \begin{equation} \Delta\boldsymbol{\sigma}=\mathbf{E}^{ep}\Delta\boldsymbol{\varepsilon} . \end{equation} \]

Bounding Surface and Mapping Rule

The deviatoric stress-ratio tensor is

\[ \begin{equation} \mathbf{r}=\frac{\mathbf{s}}{p}, \end{equation} \]

where \(\mathbf{s}\) is the deviatoric stress tensor in compression-positive convention. The tensor \(\mathbf{r}_p\) stores the stress-ratio tensor at the last loading reversal. The distance from the reversal point is

\[ \begin{equation} q_p = p\sqrt{(\mathbf{r}-\mathbf{r}_p):(\mathbf{r}-\mathbf{r}_p)} . \end{equation} \]

The normalised direction tensor is

\[ \begin{equation} \mathbf{n}=\frac{\mathbf{s}-p\mathbf{r}_p}{q_p} . \end{equation} \]

The scalar projection of the pivot stress ratio is

\[ \begin{equation} n_p = \mathbf{n}:\mathbf{r}_p . \end{equation} \]

In the implementation, positive values of \(n_p\) are set to zero for numerical robustness. The mapping variable is then calculated as

\[ \begin{equation} \zeta = \frac{q_p/p}{\sqrt{2/3}\,M_s\chi-n_p} . \end{equation} \]

A loading reversal is detected from the sign change of the first-order work measure

\[ \begin{equation} dW=(\mathbf{r}-\mathbf{r}_p):d\boldsymbol{\varepsilon} . \end{equation} \]

If the sign of \(dW\) changes, \(\mathbf{r}_p\) is updated to the current stress ratio.

Lode-Angle Dependency

The Lode-angle dependency is introduced through a polynomial interpolation of the bounding stress ratio. The implementation evaluates the Lode parameter from the deviatoric invariants as

\[ \begin{equation} \cos(3\theta)=\frac{3\sqrt{3}}{2}\frac{J_3}{J_2^{3/2}} . \end{equation} \]

For a given triaxial-compression bounding stress ratio \(M_s\), the equivalent friction angle is obtained from

\[ \begin{equation} \sin\varphi=\frac{3M_s}{6+M_s} . \end{equation} \]

The extension and simple-shear stress ratios are

\[ \begin{equation} M_e=\frac{6\sin\varphi}{3+\sin\varphi}, \qquad M_{ss}=2\sin\varphi . \end{equation} \]

The Lode-dependent stress ratio \(M_{s,\theta}\) is interpolated as

\[ \begin{equation} M_{s,\theta}=\left(\frac{M_s+M_e}{2}-M_{ss}\right)\cos^2(3\theta) +\frac{M_s-M_e}{2}\cos(3\theta)+M_{ss} . \end{equation} \]

The corresponding Lode factor is

\[ \begin{equation} \chi = \frac{M_{s,\theta}}{M_s} . \end{equation} \]

Relative Dilatancy Index

The current relative density is

\[ \begin{equation} D_r = \frac{e_{max}-e}{e_{max}-e_{min}} . \end{equation} \]

Bolton's relative dilatancy index is used as a state parameter:

\[ \begin{equation} I_r = D_r\left[Q-\ln(p)\right]-R . \end{equation} \]

The logarithm is evaluated with the same pressure unit convention used in the model calibration. The typical calibration in the cited papers uses stresses in kPa. Values \(I_r>0\) indicate denser, dilative states, whereas \(I_r<0\) indicates looser, contractive states.

The initial values \(p_0\), \(e_0\), \(D_{r0}\) and \(I_{r0}\) are stored at istep == 1 and iinc == 1 and remain constant afterwards.

Initial-State Calibration

The hardening exponent is computed from the initial relative density:

\[ \begin{equation} n = 0.4D_{r0}+0.14 . \end{equation} \]

The critical-state stress ratio in triaxial compression is

\[ \begin{equation} M_{cs}=\frac{6\sin\varphi_{cs}}{3-\sin\varphi_{cs}} . \end{equation} \]

The initial friction angle used for the initial bounding stress ratio is

\[ \begin{equation} \varphi_{s0}=\kappa\varphi_{cs}+5I_{r0} . \end{equation} \]

To avoid numerical problems, the implementation limits the corresponding initial bounding stress ratio to a value slightly below the critical-state value if \(\varphi_{s0}\geq\varphi_{cs}\). Otherwise,

\[ \begin{equation} M_{s0}=\frac{6\sin\varphi_{s0}}{3-\sin\varphi_{s0}} . \end{equation} \]

The evolution parameter is

\[ \begin{equation} c = 6+\delta I_{r0} . \end{equation} \]

The coefficients \(\kappa\) and \(\delta\) account for intrinsic fabric effects and must be calibrated for the sand under consideration.

Evolution of Bounding and Phase Transformation Stress Ratios

The peak friction angle is computed as

\[ \begin{equation} \varphi_{max}=\max\left(\varphi_{cs}+g_{\alpha b}I_r,\varphi_{cs}\right) . \end{equation} \]

The peak bounding stress ratio is

\[ \begin{equation} M_{s,peak}=\frac{6\sin\varphi_{max}}{3-\sin\varphi_{max}} . \end{equation} \]

The auxiliary stress ratio \(M_{sp}\) is obtained from the calibration relation used to place the peak of the bounding-ratio evolution. In the implementation,

\[ \begin{equation} M_{sp}=2M_{s,peak}-M_{cs}+\frac{1}{2}\sqrt{D_M}, \end{equation} \]

with

\[ \begin{equation} D_M=\max\left((2M_{cs}-4M_{s,peak})^2-16M_{s0}(M_{s,peak}-M_{cs})-4M_{cs}^2,0\right) . \end{equation} \]

The current bounding stress ratio evolves with accumulated deviatoric strain \(\sum d\varepsilon_q\):

\[ \begin{equation} M_s=M_{cs}+\left[M_{sp}+(M_{s0}-M_{sp})\exp(-c\sum d\varepsilon_q)-M_{cs}\right] \exp(-c\sum d\varepsilon_q) . \end{equation} \]

For dense states, \(I_r>0\), the initial phase transformation stress ratio is computed as

\[ \begin{equation} M_{pt0}=M_{s,peak}\zeta^n-\frac{3(0.3I_r)}{(3+0.3I_r)\chi} . \end{equation} \]

For loose states, \(I_r\leq0\), the implementation uses

\[ \begin{equation} M_{pt0}=M_{s,peak}\zeta^n . \end{equation} \]

The phase transformation stress ratio evolves as

\[ \begin{equation} M_{pt}=M_{cs}+(M_{pt0}-M_{cs})\exp\left(-\frac{1}{2}c\sum d\varepsilon_q\right) . \end{equation} \]

Stress-Induced Peak-Strength Correction

The implementation computes the intermediate principal stress parameter from the principal stresses in compression-positive convention, ordered as \(\sigma_1\geq\sigma_2\geq\sigma_3\):

\[ \begin{equation} b=\frac{\sigma_2-\sigma_3}{\sigma_1-\sigma_3} . \end{equation} \]

The peak-strength coefficient is then

\[ \begin{equation} g_{\alpha b}=3+2\sin(\pi b), \qquad 3\leq g_{\alpha b}\leq 5 . \end{equation} \]

This means that triaxial compression and extension tend towards \(g_{\alpha b}=3\), whereas plane-strain-like states with \(b\approx0.5\) tend towards \(g_{\alpha b}=5\).

Scope of stress-induced anisotropy in this implementation

The current numgeo implementation uses \(g_{\alpha b}\) for the peak friction angle. The additional reduction factor \(a_b\) involving the principal stress rotation angle is not included.

Flow Rule

The scalar dilatancy ratio used in the plastic potential is

\[ \begin{equation} d = \sqrt{\frac{2}{3}}M_{pt}\chi - \mathbf{n}:\mathbf{r} . \end{equation} \]

The plastic potential gradient follows the anisotropic Rowe-type form

\[ \begin{equation} \boldsymbol{\Phi}_g = \mathbf{n}+d\mathbf{n}^2 . \end{equation} \]

The loading-surface gradient is

\[ \begin{equation} \boldsymbol{\Phi}_f = \mathbf{n}-\frac{1}{3}(\mathbf{n}:\mathbf{r})\mathbf{I} . \end{equation} \]

The actual implementation expands these tensor equations explicitly in Voigt notation, including the required factors for shear components.

State Variables

The model stores both evolving and reference state variables. The reference state variables \(p_0\), \(e_0\), \(D_{r0}\) and \(I_{r0}\) are initialised only in the first increment of the first step. They are then reused to compute the reference calibration quantities \(n\), \(\varphi_{s0}\), \(c\), \(M_{cs}\) and \(M_{s0}\).

The evolving state variables are the current void ratio \(e\), the first-order work value from the previous substep, the pivot stress-ratio tensor \(\mathbf{r}_p\), the accumulated deviatoric strain measure \(\sum d\varepsilon_q\) and the auxiliary water bulk modulus \(K_w\).

Stress Integration

The model is integrated explicitly with adaptive substepping:

  1. The global strain increment is split into subincrements of size \(dt\,\Delta\boldsymbol{\varepsilon}\).
  2. A first stress estimate is computed from \(\mathbf{E}^{ep}\) at the beginning of the substep.
  3. A second stress estimate is computed from the updated trial stress state.
  4. The two estimates are averaged.

  5. The relative error measure

    \[ \begin{equation} R_n=\frac{1}{2}\frac{\|\Delta\boldsymbol{\sigma}_2-\Delta\boldsymbol{\sigma}_1\|} {\|\boldsymbol{\sigma}_{n+1}\|} \end{equation} \]

    controls the substep size.

  6. If the stress state exceeds the bounding surface, the stress ratio is scaled back towards the admissible surface.

The returned tangent matrix is an elastic diagonal approximation assembled from the current elastic stiffness. It is not a consistent elastoplastic tangent.

Implementation Notes

  • The model supports ntens = 4 and ntens = 6.
  • The stress vector is expected in the standard numgeo Voigt ordering [11,22,33,12,(13,23)].
  • The initial void ratio must be prescribed through the first state variable before the first increment.
  • The pressure used in \(\ln(p)\) must be positive; the implementation applies a small cut-off to maintain numerical robustness.
  • The variable \(K_w\) is stored as an auxiliary quantity for locally undrained calculations. It is updated from the current porosity and shear modulus.

References

  1. Panagiota Tasiopoulou and Nikos Gerolymos. Constitutive modeling of sand: Formulation of a new plasticity approach. Soil Dynamics and Earthquake Engineering, 82:205–221, March 2016. doi:10.1016/j.soildyn.2015.12.014. ↩↩

  2. P. Tasiopoulou and N. Gerolymos. Constitutive modelling of sand: a progressive calibration procedure accounting for intrinsic and stress-induced anisotropy. Géotechnique, 66(9):754–770, September 2016. doi:10.1680/jgeot.15.P.284. ↩↩