1 Introduction

Hypoplasticity represents a class of incrementally nonlinear constitutive models [9, 11, 16, 17, 28, 42]. Unlike elastic–plastic models, there is no need to decompose the deformation into elastic and plastic deformations. Moreover, hypoplastic models do not make use of yield and potential surfaces [37]. Hypoplastic models are characterised by simple formulation and few material parameters and have been used to simulate soil behaviour in element tests and to solve boundary value problems with FEM [14, 25].

The performance of FEM depends on the efficiency of the numerical integration of constitutive equation. Recently, numerical integration of hypoplastic constitutive equations in finite element analysis is a topic of considerable interest. The fact that hypoplastic model has a single equation makes the implementation more straightforward than elastoplastic models [4, 30]. Since hypoplastic model does not have a yield surface, the stress return mapping algorithms common for elastoplastic models are not needed. Research on the finite element implementation of hypoplastic models can be traced back to the early work by Sikora and Gudehus [26], in which a simple explicit forward Euler scheme with constant step sizes was adopted. Later, this method was used to investigate shear band formation in granular materials by Tejchman and Wu [31]. Roddeman [24] used a generalised midpoint algorithm and Heeres and de Borst [8] considered an implicit integration method together with Newton–Raphson iterative scheme for the stress integration of a hypoplastic model [32]. Tamagnini et al. [30] studied the accuracy of some explicit methods and the generalised midpoint algorithm. Recently, Ding et al. [4] showed that explicit methods with substepping and error control are suitable for the numerical implementation of hypoplastic models.

While the aforementioned integration schemes perform well prior to reaching the failure surface, none of properly handle the stress drift away from the failure surface. The reasons for this are twofold. First, hypoplastic models are very sensitive to step size used in stress integration. A too large step size may impair the convergence and stability of the numerical calculation. Second, hypoplastic model allows some stress state outside the failure surface [40]. For some well-defined hypoplastic models, Wu and Niemunis [39] showed that all accessible stress states are within a bound surface. However, the bound surface is often too far from the failure surface, which gives rise to too high strength [40]. Consequently, the error resulted from stress drift away from the failure surface can accumulate in numerical computations and eventually lead to unphysical behaviour [3, 22] and the loss of stability for a boundary value problem. For numerical calculations, it makes sense not to allow stress to wander outside the failure surface. A solution to this problem is to use the return mapping method originally proposed for elastoplastic models [12]. In this work, we will investigate when the return mapping is needed and how the return mapping can improve accuracy. In doing so, we consider a fairly general and simple hypoplastic model so that our approach can be easily adopted to handle more sophisticated hypoplastic models.

In Sect. 2, the hypoplastic constitutive model and its numerical equations are outlined. In Sect. 3, several commonly used integration schemes are introduced. In Sect. 4 the performance of the integration algorithms is examined by performing a series of triaxial compression tests. The significance of the stress correction at the failure surface is shown by two boundary value problems.

2 A simple hypoplastic constitutive model

2.1 Constitutive model

In the framework of hypoplasticity, the constitutive equation is written in two parts, representing, respectively, the reversible and irreversible behaviour. We start with the formulation by Wu and Kolymbas [38] and write the hypoplastic rate-constitutive equation as the sum of the linear and nonlinear terms of the strain rate \(\dot{{\varvec{\varepsilon }}}\)

$${\mathring{{\varvec{\sigma }}}}={\pmb {\mathscr {L}}}({\varvec{\sigma }}):\dot{{\varvec{\varepsilon }}}-{\varvec{N}}({\varvec{\sigma }})\Vert \dot{{\varvec{\varepsilon }}}\Vert$$
(1)

where the terms \({\pmb {\mathscr {L}}}\) and \({\varvec{N}}\), respectively, denote the linear and nonlinear components, \({\varvec{\sigma }}\) is the Cauchy stress tensor, and \(\dot{{\varvec{\varepsilon }}}\) is the stretching tensor. \(\Vert \dot{{\varvec{\varepsilon }}}\Vert =\sqrt{{\mathrm {tr}}(\dot{{\varvec{\varepsilon }}}^2)}\) stands for the Euclidean norm. The Jaumann rate of the Cauchy stress tensor \({\mathring{{\varvec{\sigma }}}}\) is defined in terms of the time-derivative of the Cauchy stress tensor \(\dot{{\varvec{\sigma }}}\) and the spin tensor \({\varvec{w}}\)

$${\mathring{{\varvec{\sigma }}}}=\dot{{\varvec{\sigma }}}+{\varvec{\sigma }}\cdot {\varvec{w}}-{\varvec{w}}\cdot {\varvec{\sigma }}$$
(2)

The stretching and spin tensors are related to the velocity gradient tensor through

$$\dot{{\varvec{\varepsilon }}}=\frac{1}{2}({\varvec{\nabla v}}+{\varvec{ v \nabla }}) ,\quad {\varvec{w}}=\frac{1}{2}({\varvec{\nabla v}}-{\varvec{ v \nabla }})$$
(3)

where \({\varvec{v}}\) is the velocity and \({\varvec{\nabla }}\) is the gradient operator.

Within the above framework, a simple critical state hypoplastic constitutive equation for granular materials is proposed by Wu [42], which is an improvement of an early version by Wu and Bauer [37]. This model consists of three linear terms and one nonlinear term in the stretching tensor \(\dot{{\varvec{\varepsilon }}}\).

$${\mathring{{\varvec{\sigma }}}}=C_1({\mathrm {tr}}{\varvec{\sigma }})\dot{{\varvec{\varepsilon }}}+ C_2({\mathrm {tr}}\dot{{\varvec{\varepsilon }}}){\varvec{\sigma }}+ C_3\frac{{\mathrm {tr}}({\varvec{\sigma }}\cdot \dot{{\varvec{\varepsilon }}})}{{\mathrm {tr}}{\varvec{\sigma }}}{\varvec{\sigma }}+ C_4({\varvec{\sigma }}+{\varvec{\sigma }}^{*})\Vert \dot{{\varvec{\varepsilon }}}\Vert I_e,$$
(4)

in which \(C_i\) \((i=1,2,3,4)\) are dimensionless parameters. The deviatoric stress tensor \({\varvec{\sigma }}^*\) in Eq. (4) is defined by \({\varvec{\sigma }}^*={\varvec{\sigma }}-1/3({\mathrm {tr}}{\varvec{\sigma }}){\varvec{\delta _{ij}}},\) with \({\varvec{\delta _{ij}}}\) being the Kronecker delta. \(I_e\) is adopted as the critical state function. It is through \(I_e\) that the model captures the effects of density and confining pressure on the strain–stress behaviour. Several forms for \(I_e\) can be found in the literature [7, 15, 21, 41], but in the present work, a different formulation for the critical state function \(I_e\) is proposed:

$$I_e=\left(\frac{e}{e_{crt}}\right)^{\alpha }$$
(5)

where e and \(e_{crt}\) refer to the current void ratio and critical state void ratio, respectively, and \(\alpha\) is a constitutive constant that controls the degree of strain softening. The critical state function has the value \(I_e=1\) at the critical state, greater than 1 for a loose state, and less than 1 for a dense state. The evolution of the void ratio follows the evolution of the volumetric strain, \(\dot{\varepsilon }_v={\text {tr}}(\dot{{\varvec{\varepsilon }}})\), according to the following relationship:

$$\dot{e}=(1+e)\dot{\varepsilon }_v;$$
(6)

The critical state void ratio is calculated according to Li and Wang [13]. A slightly modified form is used in this work:

$$e_{crt}=e_{co} {\mathrm {exp}} \left[{-\lambda (\frac{p}{p_a}})^{\xi }\right]$$
(7)

where \(e_{co}\), \(\lambda\), and \(\xi\) are parametric constants, and \(p={\mathrm {tr}} {\varvec{\sigma }}/3\) and \(p_a\) denote the hydrostatic pressure and atmospheric pressure (101.325 kPa) for normalisation, respectively. With the additional term \(I_e\), we ensure that both dense and loose sands can be modelled with only one set of material parameters.

Although hypoplastic Eq. (4) was developed mainly for cohesionless soils, in practice, most soils show cohesion to some extent. However, the constitutive model is able to take cohesion into consideration by simply replacing the stress tensor \({\varvec{\sigma }}\) with the following translated stress tensor [33, 42]:

$${\varvec{\sigma _c}}={\varvec{\sigma }}-p_t{\varvec{\delta }}_{ij}$$
(8)

where the translated scale \(p_t=c/{\text {tan}}\phi\) and the parameters c and \(\phi\) are, respectively, the cohesion and friction angle of the cohesive soil.

2.2 Explicit form of the failure and bound surfaces

The failure criterion is defined by vanishing stress rate for a non-vanishing strain rate. By assuming that \({\varvec{N}}\) contains the critical state function \(I_e\), we have:

$${\mathring{{\varvec{\sigma }}}}=({\pmb {\mathscr {L}}}-{\varvec{N}}\otimes \vec {\dot{{\varvec{\varepsilon }}}}):\dot{{\varvec{\varepsilon }}} =0$$
(9)

where \(\vec {\dot{{\varvec{\varepsilon }}}}=\dot{{\varvec{\varepsilon }}}/\Vert \dot{{\varvec{\varepsilon }}}\Vert\) stands for the direction of strain, and the symbol \(\otimes\) denotes an outer product between two tensors. By making use of the fact that \(\vec {\dot{{\varvec{\varepsilon }}}}:\vec {\dot{{\varvec{\varepsilon }}}}=1\), the failure criterion can be readily derived:

$$f({\varvec{\sigma }})=\Vert {\pmb {\mathscr {L}}}^{-1}:{\varvec{N}}\Vert -1=0,$$
(10)

Therefore, the condition of invertibility is \(\Vert {\pmb {\mathscr {L}}}^{-1}:{\varvec{N}}\Vert <1\). This condition is identical with the requirement that a stress state \({\varvec{\sigma }}\) should lie inside the failure surface given in Eq. (10), which means constitutive Eq. (4) is invertible when the stress lies inside the failure surface.

The explicit formula of the failure surface can be obtained using the symbolic computational programme Mathematica, which gives rise to the failure surface:

$$f({\varvec{\sigma }})=\sqrt{J_2}+ \varsigma I_1=0,$$
(11)

where \(J_2\) and \(I_1\) are, respectively, the second deviatoric stress invariant and the first stress invariant, and \(\varsigma\) is a constant determined by the dimensionless parameters \(C_i(i=1,2,3,4)\) and the critical state function \(I_e\), which can be found in “Appendix”.

The hypoplastic model is also characterised by a bound surface that bounds the accessible stress states and which can be derived based on the procedure by Wu [36] and Wu and Niemunis [40]. The explicit formulation of the bound surface is expressed as follows:

$${\mathcal {B}}({\varvec{\sigma }})=\sqrt{J_2}-\frac{C_1}{\sqrt{(I_eC_4)^{2}-3C_1^2}}I_1=0,$$
(12)

With the help of Eqs. (11) and (12), the failure and bound surfaces can be plotted once the model parameters are given.

2.3 Second-order work and stability surface

In view of the complexity of the nonlinear constitutive models, it is desirable to obtain qualitative properties such as stability and uniqueness of the boundary value problems posed with the constitutive models [40]. The problem of stability can be approached based on the analysis of the second-order work. Instability considered in terms of the second-order work virtually means the possibility for spontaneous increase in kinetic energy of the body due to a small disturbance in velocity [19]. According to Hill [10], a sufficient condition for stability is the second-order work

$${\mathcal {W}}_2({\varvec{\sigma }})={\mathrm {tr}}({\mathring{{\varvec{\sigma }}}}\dot{{\varvec{\varepsilon }}}) >0$$
(13)

for all directions of stretching. For hypoplastic models, the second-order work can become negative before the failure surface [35]. Let us consider constitutive equation (4) and search for the boundary between positive and negative second-order work by letting \({\mathcal {W}}_2=0\):

$${\mathcal {W}}_2({\varvec{\sigma }})=\dot{{\varvec{\varepsilon }}}:({\pmb {\mathscr {L}}}-{\varvec{N}}\otimes \vec {\dot{{\varvec{\varepsilon }}}}):\dot{{\varvec{\varepsilon }}}=0$$
(14)

If this boundary builds a surface in the stress space, it will be called stability surface. By using the analytical approach proposed by Niemunis [18], an explicit expression of the stability surface can be readily derived to be:

$$\begin{aligned} {\mathcal {W}}_2({\varvec{\sigma }})&= (L_1^2-N_1^2)\left[4L_3(L_3+L_1)+N_1^2\right]^2-N_3^2\Big[16L_1^4-32L_1^3L_3 \nonumber \\&\quad +3N_1^4+3N_1^2N_3^2+N_3^4+2L_1L_3(4L_3^2+19N_1^2-4N_3^2) \nonumber \\&\quad -L_3^2(20N_1^2+N_3^2)+4L_1^2(2L_3^2-5N_1^2+2N_3^2)\Big] \end{aligned}$$
(15)

where the variables are expressed as function of the stress invariants \(I_1\) and \(J_2\): \(L_1= C_1I_1\), \(L_3= (C_1+C_2+C_3/3)I_1+2C_3J_2/I_1\), \(N_1=- C_4 I_1 \sqrt{2J_2}/\sqrt{I_1^2+6 J_2}\), and \(N_3=C_4 \left( I_1^2+12 J_2\right) /\sqrt{3I_1^2+18 J_2}\). Note that the stability surface depends on the density through the critical state function \(I_e\) by replacing \(N_1\) and \(N_3\) with \(I_eN_1\) and \(I_eN_3\), respectively. The implication of the stability surface for FE calculations is beyond the scope of this paper. In general, the stability should be considered within a BVP rather than an element. Some stability indicators can be introduced to avoid pitfalls in the FE calculations [23].

By using Eq. (15), we plot the stability surface together with the bound and failure surfaces in the principal stress space. As shown in Fig. 1, the stability surface of the hypoplastic model is a cone with its apex at the origin in the principal stress space. The bound and failure surfaces possess similar geometry, outside the stability surface. Wu and Niemunis [40] show that a stress state may lie outside the failure surface for some paths. We proceed to investigate how far a stress state can drift away from the failure surface for constitutive equation (4).

Fig. 1
figure 1

The bound, failure and stability surfaces in principal stress space

Full size image

Let us consider three stress states \({\varvec{\sigma }}_b\), \({\varvec{\sigma }}_f\) and \({\varvec{\sigma }}_s\) with the same Lode angle lying on the bound, the failure, and the stability surfaces, respectively. The corresponding stress ratio can be obtained as follows:

$$\frac{\sigma _{1}-\sigma _{3}}{\sigma _{1}+\sigma _{3}}={\text {sin}}\phi _\mathrm{{mob}}$$
(16)

where \(\sigma _1\) and \(\sigma _3\) are the maximum and minimum principal stresses, respectively. The stress ratio can be obtained from the mobilised friction angle \(\phi _\mathrm{{mob}}\). In contrast to the Mohr–Coulomb yield criterion, the mobilised friction angle \(\phi _\mathrm{{mob}}\) varies for different Lode angles for the hypoplastic model. Obviously, the mobilised friction angles obtained from \({\varvec{\sigma }}_b\), \({\varvec{\sigma }}_f\) and \({\varvec{\sigma }}_s\) are dependent on the Lode angle. Figure 2 shows the three mobilised friction angles \(\phi _{\mathrm{{mob}}\_b}\) \(, \phi _{\mathrm{{mob}}\_f}\) and \(\phi _{\mathrm{{mob}}\_s}\) for different critical friction angles. For the stress state in triaxial compression (the Lode angle is zero) with a \(20^{\circ }\) frictional angle, one obtains \(\phi _{\mathrm{{mob}}\_b} =21.4^{\circ }\), \(\phi _{\mathrm{{mob}}\_f} =20^{\circ }\) and \(\phi _{\mathrm{{mob}}\_s} =11.7^{\circ }\) while for the stress state in extension (the Lode angle is 60°) the corresponding friction angles are \(\phi _{\mathrm{{mob}}\_b} =28.8^{\circ }\), \(\phi _{\mathrm{{mob}}\_f} =26.3^{\circ }\) and \(\phi _{\mathrm{{mob}}\_s} =16.8^{\circ }\).

Fig. 2
figure 2

The mobilised fiction angle \(\phi _\mathrm{{mob}}\) for different critical frictional angles \(\phi\): a \(\phi =20^{\circ }\), b \(\phi =25^{\circ }\), and c \(\phi =30^{\circ }\)

Full size image
Table 1 The mobilised friction angles \(\phi _\mathrm{{mob}}\) at different stress states
Full size table

The bound surface is an intrinsic property of hypoplastic constitutive equations, which thus some advantage in the numerical integration over most conventional constitutive models, since the stress states lying outside the bound surface will be automatically corrected in the next time step [40]. However, a stress state may also lie between the failure and bound surface for some strain paths. Table 1 shows the mobilised friction angles for the failure surface and bound surface in triaxial tests. We observe that the difference of the mobilised friction angles between the failure and bound surface increases with the critical friction angle and may reach 30° for triaxial extension, which is far too high to be accepted. Obviously, stress corrections are needed to bring the stress back to the failure surface.

3 Time integration of the constitutive model

3.1 Tangential stiffness for FEM

The finite element formulation for the linear continuum model follows the standard Galerkin approximation of the weak form based on the virtual work principle. Thus the global finite element equation can be expressed:

$$\sum \int _{V_e}{\varvec{\mathrm {B}}}^{\mathrm {T}}{\varvec{\sigma }}({\varvec{u}}_n) dV_e= {{\varvec{\mathrm {P}}}}_{n\_\rm{ext}}$$
(17)

where n denotes the nth increment of the nonlinear analysis. \({\varvec{P}}_{n\_\rm{ext}}\) is the nodal force vector, and \({\varvec{u}}\) is the unknown nodal displacement vector. Thus the unbalanced force vector can be expressed as:

$${{\varvec{\mathrm {F}}}}({\varvec{u}}) = \sum \int _{V_e}{{\varvec{\mathrm {B}}}}^{\mathrm {T}}{\varvec{\sigma }}({\varvec{u}}_n) dV_e- {\varvec{\mathrm {P}}}_{n\_\rm{ext}}$$
(18)

For the ith iteration in the framework of the standard Newton–Raphson procedure, the global stiffness matrix reads:

$${{\varvec{\mathrm {K}}}}({\varvec{u}}_n^i)= \frac{\partial {{\varvec{\mathrm {F}}}}({\varvec{u}}_n^i)}{\partial {\varvec{u}}} =\sum \int _{V_e}{{\varvec{\mathrm {B}}}}^{\mathrm {T}} \frac{\partial {\varvec{\sigma }}}{\partial {\varvec{\varepsilon }}}\Big|_{{\varvec{u}}_n^i} {{\varvec{\mathrm {B}}}} dV_e$$
(19)

For the hypoplastic constitutive model, the Jacobian matrix can be readily written out:

$${{\varvec{\mathrm {D}}}}=\frac{\partial \Delta {\varvec{\sigma }}}{\partial \Delta {\varvec{\varepsilon }}} = {\pmb {\mathscr {L}}}+I_e{\varvec{N}} \vec {\dot{{\varvec{\varepsilon }}}}$$
(20)

For constitutive model (4), the tangential stiffness tensor at a given stress \(\sigma _{ij}\) is given below:

$$\frac{\partial \mathring{\sigma }_{ij}}{\partial \dot{\varepsilon }_{mn}}= C_1\Theta \ I_{ijmn}+ C_2 \sigma _{ij}\delta _{mn}+ C_3\frac{\sigma _{ij} \sigma _{mn}}{\mathrm {tr}{\varvec{\sigma }}}+ I_eC_4 (\sigma _{ij}+\sigma _{ij}^{*})\vec {\dot{\varepsilon }}_ {mn}$$
(21)

where \(I_{ijmn}\) is the fourth rank identity tensor with components \(I_{ijmn}=(\delta _{im}\delta _{jn}+\delta _{jm}\delta _{in})/2\). It is worth noting that the Jacobian matrix may become negative definite when using the above tangential operator. A numerical study on the performance of different tangential operators can be found in [4].

3.2 Stress integration algorithms

The constitutive equation can be regarded as an ordinary differential equation, for which the general time integration over an increment step \(t\in [t_n, t_{n+1}]\) can be written as:

$${\varvec{\sigma }}_{n+1}={\varvec{\sigma }}_{n}+\int _{t_n}^{t_{n+1}}h({\varvec{\sigma }},e,\dot{{\varvec{\varepsilon }}})dt = {\varvec{\sigma }}_{n}+{\varvec{\mathrm {D}}}\Delta {\varvec{\varepsilon }}, \quad n=1,2,\ldots$$
(22)

According to Eq. (6), a closed form of integration for the void ratio is available:

$$e_{n+1}= (1+e_n)\cdot {\text {exp}}(\Delta \dot{\varepsilon }_v) -1$$
(23)

The solution for Eqs. (22) and (23) can be obtained stepwise according to one of the integration algorithms in the next subsection.

3.2.1 Generalised midpoint algorithms

Following the method proposed in [5, 6], all stress components and state variables are collected in vector \({\varvec{y}}\) for the simultaneous integration of Eqs. (22) and (23):

$${\varvec{y}}=\{\sigma _{11},\sigma _{22},\sigma _{33},\sigma _{12},\sigma _{13},\sigma _{23},g_1\dots g_m\}^{\mathrm {T}}$$
(24)

where \(g_i(i=1 \ldots m)\) are additional state variables, e.g. void ratio. The integration of Eq. (24) requires the solution of initial value problem:

$${\varvec{y}}'(t)=H({\varvec{y}}(t)),\quad {\varvec{y}}'(0)={\varvec{y}}(0).$$
(25)

Generalised midpoint algorithms are among the most widely used second-order integration methods [2]. The general form of the generalised midpoint method (e.g. [30]) can be written as:

$${\varvec{y}}_{n+1}={\varvec{y}}_{n}+\Delta t_{n+1}[(1-\theta )\dot{{\varvec{y}}}_n+\theta \dot{{\varvec{y}}}_{n+1}]\quad n=1,2,\ldots$$
(26)

where \(\Delta t_{n+1}= t_{n+1}-t_{n}\) is the time step increment and the parameter \(\theta\) is a prescribed constant within [0, 1]. Note that generalised midpoint algorithms with values of \(\theta\) equal to 1 and 0, respectively, correspond to the implicit backward Euler method and explicit forward Euler method with the Crank–Nicolson (midpoint or trapezoidal) method obtained by setting \(\theta = 0.5\). Our implementation tests three widely used schemes: the explicit forward Euler method (FE), the modified Euler method (ME) and the implicit Crank–Nicolson method (CN).

3.2.2 Adaptive explicit algorithms

The accuracy of an integration method can be improved by reducing the size of the time increment. Although this can be carried out in a straightforward manner by dividing the time increment into several equal substeps, the better accuracy is usually gained at the cost of computational time or failure of error control within a tolerant range. A more powerful approach is to employ the adaptive integration schemes described in the literature [27, 34], which enable users to adjust the substep size automatically according to the local truncation error. Studies have revealed that this approach has the merits of being efficient and robust for a wide range of constitutive models. In the present paper, several adaptive explicit integration methods, namely the modified Euler substepping scheme (MEsec), the Richardson extrapolation substepping scheme (REsec) and the Runge–Kutta–Fehlberg substepping scheme (RKFsec), are implemented and compared.

To compute the local error at each substep of the stress integration, two approximate solutions with different orders of accuracy (p, q) are obtained and compared. If the two solutions are in close agreement, the approximation is accepted; otherwise, the substep is rejected and the corresponding step size is further reduced. For the generic substep k in the time interval \([t_n, t_{n+1}]\), with dimensionless size \(\Delta T_k\in (0, 1]\) given by the following equation:

$$\Delta T_k=(t_{k+1}-t_k)/(t_{n+1}-t_n) \leqslant 1 \quad {\text {and}} \quad \sum _{k=1}^{n_s}\Delta T_k=1,$$
(27)

two different approximate solutions of the evolution problem (25) are obtained simultaneously according to

$${\varvec{y}}^{(p)}_{k+1} = {\varvec{y}}_{k}+\Phi _1({\varvec{y}}_{k},\Delta {\varvec{\varepsilon }}_{n+1},\Delta T)\Delta (T_k)^{(p+1)}$$
(28a)
$${\varvec{y}}^{(q)}_{k+1} = {\varvec{y}}_{k}+\Phi _2({\varvec{y}}_{k},\Delta {\varvec{\varepsilon }}_{n+1},\Delta T)\Delta (T_k)^{(q+1)}$$
(28b)

The two function \(\Phi _1\) and \(\Phi _2\) are constructed as follows according to Table 2, which summarises different integration methods:

$$\Phi _1 {:=} \sum _{i=0}^{p} k^{(p)}_i h_i({\varvec{y}}_{k},\Delta {\varvec{\varepsilon }}_{n+1},\Delta T)$$
(29a)
$$\Phi _2 {:=} \sum _{i=0}^{q} k^{(q)}_i h_i({\varvec{y}}_{k},\Delta {\varvec{\varepsilon }}_{n+1},\Delta T)$$
(29b)

where

$$h_i({\varvec{y}}_{k},\Delta {\varvec{\varepsilon }}_{n+1},\Delta T):=h({\varvec{y}}_{k}+\Delta T \sum _{j=0}^{l-1}\vartheta _{lj}h_j,\Delta {\varvec{\varepsilon }}_{n+1} )$$
(30)

For the sake of simplicity, we usually set \(q = p + 1\). The constants \(k^{(p)}_i\), \(k^{(q)}_i\) and \(\vartheta _{lj}\) are used to obtain the approximated solutions of order p and q, respectively. Then the local truncation error of the lower-order method at time \(T_{k+1}\) can be obtained by using the difference in the above two approximate solutions:

$${\varvec{R}}_{k+1}= {\varvec{y}}^{(p)}_{k+1}-{\varvec{y}}^{(q)}_{k+1}, \quad {\text {and}} \quad R_{k+1}=\frac{\Vert {\varvec{R}}_{k+1} \Vert }{\Vert {\varvec{y}}^{(q)}_{k+1} \Vert }$$
(31)

The integration over the k-th substep is assumed to be successful when, for a given stress error tolerance STOL,

$$R_{k+1}\leqslant {\text {STOL}},$$
(32)

with the new substep size then estimated using the following extrapolation formula:

$$\Delta T_{k+1}=\Delta T_{k} \left[\frac{{\mathrm {STOL}}}{R_{k+1}}\right]^{1/(p+1)}$$
(33)

If the estimated error is less than the prescribed accuracy tolerance STOL, the step is accepted and we enlarge our step size according to

$$\Delta T_{k+1}=\Delta T_{k}\cdot {\text {min}} \left\{1.1, 0.9 \left[\frac{{\mathrm {STOL}}}{R_{k+1}}\right]^{1/(p+1)} \right\}$$
(34)

If condition (32) is not satisfied, the k-th substep will be re-calculated with a smaller step size \(\Delta T_{k}^{*}\):

$$\Delta T_{k}^{*}=\Delta T_{k} \cdot {\text {max}} \left\{ 0.25 , 0.9 \left[\frac{{\mathrm {STOL}}}{R_{k+1}}\right]^{1/(p+1)}\right\}$$
(35)

Note that the right side of Eq. (33) is multiplied by a factor that is typically set to 0.9. An upper bound of 1.1 and a lower bound of 0.25 are also taken for each new substep in order that the extrapolation is not carried too far. More details regarding the substepping algorithm can be found in the literature [1, 29, 30]. After the integration process is complete, the stress tensor can be extracted from the vector \({\varvec{y}}\).

3.2.3 Correction of stresses to failure surface

At the end of each increment in the integration process, the stresses may diverge from the failure function so that \(f({\varvec{\sigma }})> FTOL\). The extent of this violation, which is commonly known as failure surface ‘drift’, depends on the accuracy of the integration scheme and the nonlinearity of the constitutive relations. Sloan [27] suggested that, provided the integration is performed accurately, the extent of drift from the failure surface tends to be small and no remedial action is required. Wu and Niemunis [40] and Niemunis [18], on the other hand, have reported that some stress states may surpass the failure surface irrespective of the accuracy of the used integration method. In such cases, the stress state does not satisfy the consistent condition.

Fig. 3
figure 3

Sketch of return mapping scheme, a return to the cone, and b return to the apex, the direct of \({\varvec{\sigma }}\) is negative

Full size image

Let us consider a stress state inside the failure surface, e.g. stress \({\varvec{\sigma }}_{n}\) at the nth step. As demonstrated in [39, 40], the hypoplastic model allows some stress state lying outside the failure surface. For the strain path shown in Fig. 3, no matter how accurate the integration method, the stress defined by \({\varvec{\sigma }}^{\rm{trial}}_{n+1}\) at the \((n+1)\)th step of analysis will violate the consistent condition, so that

$$f({\varvec{\sigma }}^{\rm{trial}}_{n+1})= \sqrt{J_2({\varvec{\sigma }}^{*\rm{trial}}_{n+1})} + \varsigma (C_i, I_e) I_1 >FTOL,$$
(36)

with the explicit formulation defined in “Appendix”.

As any deviations from the failure surface are cumulative and may result in unacceptable errors in subsequent computations, the stresses should be corrected to satisfy the current consistent condition. As shown in Fig. 3a, for \(p^{\rm{trial}}_{n+1} < 0\), the stress is corrected along the radial direction to the failure surface [35]. With the radial return scheme, the corrected stress state takes the following form:

$${\varvec{\sigma }}^*_{n+1}=\eta {\varvec{\sigma }}^{*\rm{trial}}_{n+1} , \quad p_{n+1}=p^{\rm{trial}}_{n+1},$$
(37)

where \(\eta\) is an unknown multiplier. Using the previous definition of \(J_2\), it follows that:

$$\frac{{\varvec{\sigma }}^{*}_{n+1}}{J_2({\varvec{\sigma }}^{*}_{n+1})}=\frac{{\varvec{\sigma }}^{*\rm{trial}}_{n+1}}{J_2({\varvec{\sigma }}^{*\rm{trial}}_{n+1})}.$$
(38)

In order to return the stress state to the failure surface, it is desirable that the total strain increment, \(\Delta {\varvec{\varepsilon }}\), remains unchanged, since this is consistent with the displacement-based finite element procedure. The corrected stress state in Eq. (37) satisfies the consistency condition. Using Eqs. (36) and (38), together with the assumption that departures from the failure surface are sufficiently small that only one return step is required, the consistent condition is expressed as follows:

$$f= \eta \sqrt{J_2({\varvec{\sigma }}^{*\rm{trial}}_{n+1})} + \varsigma (C_i,I_e) I_1 =0,$$
(39)

which yields the unknown multiplier \(\eta\).

On the other hand, if the mean stress \(p^{\rm{trial}}_{n+1} > 0\), as shown in Fig. 3b, the stress is corrected to the apex of the failure surface

$${\varvec{\sigma }}_{n+1}=0 , \quad p_{n+1}=0,$$
(40)

After solving the above equation, we correct the violated stress either to the cone or to the apex of the failure surface via Eqs. (37) and (40).

Fig. 4
figure 4

The sketch of corrected Stress Response Envelope (SRE)

Full size image

The stress correction scheme can also be easily adapted to incorporate the effect of critical state and cohesion. To this end, the constant \(\varsigma\) including the critical state function \(I_e\) can be found in “Appendix”. Figure 4 shows the corrected Stress Response Envelope (SRE). A detailed study regarding SRE can be found in [30]. For any stress lying between the failure and bound surfaces, the updated stress is forced back to the failure surface and thus the function of the bound surface is abandoned.

4 Numerical tests for different integration strategies

To obtain an overall assessment of the integration methods presented in Sect. 3, a set of numerical tests is conducted for constitutive equation (4). Firstly, drained and undrained triaxial compression tests are modelled. We study the influence of stress correction on the stress–strain relation in drained and undrained triaxial tests, together with an analysis of the stress drift from the failure surface. Secondly, two boundary problems, namely a rigid footing, and the safety factor of a homogeneous slope, are considered using the finite element code Abaqus Standard with user-defined materials (Umat). The performance of different integration methods and stress correction scheme in solving these problems is then compared in terms of accuracy, efficiency, and robustness.

4.1 Performance of integration methods

Two groups of integration methods are analysed in the numerical triaxial compression tests:

  1. (1)

    Constant substep method: the forward Euler method (FE), the modified Euler method (ME), and the Crank–Nicolson method (CN);

  2. (2)

    Adaptive explicit method: the adaptive modified Euler method (MEsec), the adaptive Richardson extrapolation method (REsec), and the adaptive Runge–Kutta–Fehlberg method (RKF23sec).

The formulation of different integration methods and error estimations are shown in Table 2.

Table 2 Different integration methods for generating estimated errors
Full size table
Table 3 Parameters for numerical simulation of the drained triaxial test
Full size table
Fig. 5
figure 5

Stress–strain relations of triaixal compression tests with a, b 10 increments, 2 substeps, STOL = \(10^{-4}\) and c, d 20 increments, 1 substeps, STOL = \(10^{-4}\)

Full size image

To assess their numerical performance, a benchmark solution is obtained using the adaptive RKF45 method, in which the integration error tolerance is set to \(10^{-9}\). The relative error is calculated for every step as follows:

$$R_n=\frac{\Vert {\varvec{y}}^n_{\rm{exact}}-{\varvec{y}}^n\Vert }{\Vert {\varvec{y}}^n_{\rm{exact}}\Vert }. \qquad n=n \mathrm{{th}} \;{\text {step}}$$
(41)

For the triaxial compression tests, an initial isotropic stress state with \(\sigma _{11}=\sigma _{22}=\sigma _{33}=100\) kPa is assumed. The initial void ratio is set to \(e_i\) = 0.78 for the drained triaxial test and \(e_i\) = 0.93 for the undrained triaxial test; both tests are strain-controlled with a maximum axial (vertical) strain of \(10\%\) applied, and the horizontal strain increment is calculated by the constitutive model for a given axial strain increment. The parameters used in the simulations are shown in Table 3. In the numerical procedures, two types of increments are adopted. In the first calculation, the loading process is divided into 10 equal increments, representing a large increment size scheme. In the second calculation, the loading process is divided into 20 equal increments, representing a relatively small increment size scheme. For each scheme, different substeps are performed for the explicit Euler method and the implicit CN method, and different STOLs applied for the CN method. Similarly, the integration error tolerance STOL is changed for the adaptive explicit methods. For each method, the integration results are accepted once convergence is obtained or the iteration number limit reached. The numerical results of the drained and undrained triaxial tests carried out using the various integration methods with 10 increments (2 substeps) and 20 increments (1 substeps) are shown in Fig. 5. The maximum errors obtained from computations with 10 increments under different substeps (substep = 2, 20, 100) for constant substep methods and under different STOLs for adaptive explicit methods are shown in Fig. 6.

Fig. 6
figure 6

The maximum stress error obtained from computaion with 10 increments using a Constant substep methods, and b Adaptive explicit methods

Full size image

As shown in Figs. 5 and 6, the integration strategies exhibit very different behaviour. The forward Euler (FE) method with 2 substeps provides the roughest estimation of the stress–strain response in both the drained and undrained triaxial tests. Indeed, the relative error produced by this scheme reaches 0.029 and 0.6491 in the drained and undrained tests, respectively, which can easily lead to unacceptable results in finite element calculations. The poor performance in the undrained tests is due to the used hypoplastic constitutive model which is very sensitive to both time step and stress path. Comparison of Figs. 5b and 4d reveals that the main error results from the first increment, and this error is accumulated to the rest increments. However, by increasing the number of substeps, the maximum error resulted from FE method can be largely decreased, as shown in Fig. 6a.

The performance of the simple Modified Euler method (ME) is rather poor in undrained triaxial test. The implicit Crank–Nicolson method (CN) with 2 substeps is more accurate than either of the above explicit methods. Among the three adaptive explicit methods, the modified Euler method and the RKF23 method achieve the highest accuracy, with relative errors around of \(10^{-5}\), whereas the accuracy of the REsec method with STOL = \(10^{-1}\) is relatively low. When STOL is decreased, all adaptive explicit methods perform well; the stress errors decrease as STOL is reduced, with the best accuracy achieved with the error tolerance set to \(10^{-6}\), as shown in Fig. 6b. In all, the performance of the integration methods is better in the drained tests than in the undrained tests, which implies that the integration methods are dependent on the stress path.

4.2 The effect of stress correction

To evaluate the effect of stress correction, three different initial isotropic stress states with \(\sigma _{11}=\sigma _{22}=\sigma _{33}=\) 100/200/300 kPa are assumed. The initial void ratio is set to \(e_i\) = 0.78 for the drained triaxial test and \(e_i = 0.95\) for the undrained triaxial test, with both tests strain-controlled with a maximum axial (vertical) strain of \(20\%\). The parameters used in these simulations are shown in Table 3. We consider only the adaptive RKF23sec method. The integration error tolerance is set at STOL = \(10^{-4}\). Stress correction is employed in each step to bring the stress back to the failure surface.

Fig. 7
figure 7

Results of stress correction in triaxial tests: Stress–strain relations of a drained triaxial test, b undrained triaxial test

Full size image

Figure 7 presents the results of stress correction in both the drained and undrained triaxial tests. As can be observed from Fig. 7a, stress correction takes effect only after the peak in the stress–strain curve, which implies that stress correction is particularly relevant in the softening regime. However, stress correction does not make noticeable difference in the undrained test, as shown in Fig. 7b. The stress path recorded in the drained and undrained tests are presented in Fig. 8a. It can be observed from Fig. 8a that whereas the undrained stress paths do not surpass the critical state line, the drained stress paths not only exceed the critical state line but also reach the peak stress state line. This reveals that stress correction mainly takes place in the domain between the critical state line and the peak stress state line. Figure 8b shows that the failure function is violated, i.e. \(f>0\) beyond the peak. However, with the adoption of the stress correction scheme, the failure function f becomes null, thus guaranteeing that the stress lies on the failure surface.

Fig. 8
figure 8

a Stress path of triaxial tests, b the value of failure function f in drained triaxial test

Full size image
Fig. 9
figure 9

The effect of stress correction on drained triaxial tests (p = 50 kPa): a Stress–strain relation with various internal friction angles, and b relative error after stress correction

Full size image

As discussed in Sect. 2, stress drift is particularly relevant for large critical friction angle, e.g. \(\phi =30^{\circ }\). In order to explore this phenomenon, the effect of stress correction for materials with various friction angles is thus investigated. Figure 9a, which illustrates the stress–strain relationship for different friction angles, reveals that the magnitude of the corrected stress increases with an increase in the friction angle. Figure 9b shows the relative error of stress correction in the drained triaxial tests. The relative error increases from 0.5% for \(\phi =15^{\circ }\) to 4.5% for \(\phi =45^{\circ }\). Obviously, stress correction is important in the numerical implementation of the hypoplastic constitutive models.

4.3 Boundary value problems

We now consider two typical boundary value problems, i.e. rigid footing and slope stability. In the rigid footing test, attention is focused not only on accuracy and robustness, but also on the computational efficiency of the numerical schemes. The second problem of slope stability is particularly relevant to bring out the effect of stress correction for the numerical calculations. As shown in the last section, the simple explicit and implicit methods with large step sizes can produce large errors, a constant substepping method is considered together with the forward Euler method, modified Euler method and Crank–Nicolson method. For the latter, different stress tolerances are used and the maximum number of local iterations is set to 10,000. For the adaptive explicit methods, the maximum number of substeps is less than 10,000 and the minimum substeps size is less than 10\(^{-7}\) of the current increment size. For both problems, benchmark solutions are obtained via the RKF45sec method with STOL = 10\(^{-9}\), which is compared to the numerical solutions obtained by the different methods with and without stress correction. To accommodate tensile stresses, a cohesion c is assigned to the soil, thereby allowing the development of tensile stresses during the computations.

4.3.1 Rigid footing test

Further investigation of the above numerical methods is carried out for the boundary value problem of a rigid footing. The computation domain, as shown in Fig. 10, is 4.0 m deep by 12 m wide and the width of the footing is \(w = 1.2\) m.

For the sake of simplicity, an asymmetric model is chosen using a total number of 150 four-node plane strain elements, and 600 Gauss integration points. The maximum vertical displacement is \(d = 0.5\) m, at which point the vertical force reaches its peak value, with the displacement divided into 100 equal increments. Prior to loading of the footing, an initial geostatic stress (120 kPa) is applied. The parameters used in this simulation are shown in Table 4 for an initial void ratio of \(e_i=0.78\). Stress integration errors are evaluated at the end of calculation. Note that the explicit Euler method with stress correction is not involved in the evaluation of stress error. The numerical results displayed in Table 5, in which the “Total number of substeps” is calculated according to the accumulated number of substeps for all Gauss points across all increments, while the “Maximum number of substeps” represents the substeps of one Gauss point with the maximum number of substeps across all increments.

Fig. 10
figure 10

Finite element meshes of the rigid footing

Full size image
Table 4 Parameters for the rigid footing test
Full size table
Table 5 Results of different methods for rigid footing test (100 increments)
Full size table
Fig. 11
figure 11

a The contour of the number of substep at the 3rd increament, and b the strain contour of the foundation (RKF45sec method with STOL = \(10^{-9}\))

Full size image

Table 5 shows that while the FE method and ME method have similar average errors, the latter method requires twice the CPU time of the former. In addition, the stress correction accounts for a very small proportion of the total CPU time. In contrast, the CN method with 20 substeps produce less accurate result than the CN method with 100 substeps, although the former took less time than the latter. As expected, the adaptive explicit methods are able to control the integration error and CPU time cost effectively for a given STOL. Among the two adaptive methods, the MEsec method is more efficient and the RKF23sec method more accurate, both with excellent performance.

A colour contour plot of the number of substeps for each element is shown in Fig. 11a, which displays a well-defined shear band developed near the footing. As expected, the average number of substeps is much higher in the shear band, which indicates that the substepping scheme reduces the calculation effort. Figure 11b shows a triangular zone of large strain directly under the foundation, as well as a radial zone and a Rankine passive zone, forming together the three zones assumed by Terzaghi.

Fig. 12
figure 12

Relations of vertical force and vertical displacement

Full size image

The relationship between vertical force and vertical displacement is shown in Fig. 12. First, let us look at the calculations with one substep. The calculation without stress correction becomes instable near the ultimate load, while the calculation with stress correction remains stable even beyond the peak. The calculations with 10 constant substeps show that the performance can be improved by substepping, and the stress correction gives rise to slightly lower limit load.

4.3.2 Stability of homogeneous slope

The stress correction scheme is further validated by evaluating the safety factor of a homogeneous slope and simulating the subsequent failure process. In slope stability analysis, the safety factor is typically evaluated using the so-called shear strength reduction technique, in which the shear strength (friction angle and cohesion) is reduced by a reduction factor until slope failure occurs; the safety factor is thus defined by this reduction factor. Here the safety factor and failure of a homogeneous slope are studied using different integration methods, i.e. the implicit CN method, the FE method, the adaptive RKF23 method and each of these methods with stress correction treatment. Recently, the shear strength reduction technique is used in SPH to study slope stability with our hypoplastic constitutive model [20].

Table 6 Parameters for the homogeneous slope
Full size table
Fig. 13
figure 13

Geometry and boundary conditions of the slope

Full size image

The geometry and boundary conditions of the considered slope are shown in Fig. 13. The considered slope is discretised by 349 four-node plane strain elements with the bottom of the slope fixed in the horizontal and vertical directions, and the lateral boundaries fixed only in the horizontal direction. The slope is assumed to consist of cohesive soil with the material parameters listed in Table 6, including an initial void ratio of \(e_i=0.88\). The friction angle \(\phi\) and cohesion c are the two shear strength parameters subjected to strength reduction. In the searching process, the actual shear strength is reduced by a factor, that is,

$$\phi _f = \frac{\phi }{F_s}, \quad c_f = \frac{c}{F_s},$$
(42)

with the reduced shear strength parameters then used to compute the corresponding hypoplastic parameters \(C_1 ,C_2,C_3, C_4\) according to the procedure given in [33]. To obtain the initial stress state, a geostatic step is performed by applying gravity loading to the soil. During the geostatic loading, the factor \(F_s\) is maintained at a constant value of 0.5 to avoid failure. In the second step, the shear strength parameters are reduced by increasing the factor \(F_s\) from 0.5 to 2.0, with the initiation of slope failure defined as occurring at the time when the computation is non-converging. This procedure enabled the adoption of a feasible global increment in the simulations, the numerical results are presented in Table 7.

Table 7 Results of different methods for safety factor of the slope
Full size table

It can be observed from Table 7 that all the tested explicit methods, with or without stress correction treatment, produced a safety factor of approximately 1.2, corresponding to the shear parameters of about \(\phi = 16.7^{\circ }\) and \(c = 10\) kPa. However, the various integration methods exhibited very different levels of performance. Although the forward Euler (FE) and forward Euler with stress correction (FEs) methods both complete around 480 increments of the computation, with every loading increment being divided into 100 equal substeps, the FEs method took approximately 100 s longer than the FE method and also produced a smaller safety factor. In contrast, the implicit CN method with 100 equal substeps took more than 500 increments to obtain the safety factor, while an increase in error tolerance did not result in a corresponding increase in computation time for this method. The CN method with stress correction failed at the first increment, with non-convergence recorded during the local iterations, implying that local iteration is very sensitive to the stress state. Obviously, stress correction may result in a large difference between the current and last substep of the local iterations. Similarly, the adaptive explicit method with or without stress correction required more than 500 increments to obtain the safety factor. Using this method the total number of substeps for all Gauss points is much higher than that for the FE method and the implicit CN method, although the total CPU time is lower than that required for the other methods without substepping schemes. This further indicates that adaptive explicit methods can effectively save CPU time due to their adaptive nature. Table 7 also reveals that the adaptive explicit method with stress correction, e.g. the RKF23secs method, took less time to complete the calculations than the RKF23sec method without stress correction.

Fig. 14
figure 14

Change of the horizontal displacement at the top of the slope by strength reduction using different integration strategies

Full size image
Fig. 15
figure 15

The displacements (a, b) and strains (c, d) upon slope failure by strength reduction using the RKF23sec method, (I) with and (II) without stress correction

Full size image

The above results imply that stress correction has a significant influence on numerical computations using hypoplastic models. Figure 14 presents the change in horizontal displacement at the top of the slope (point A in Fig. 13) under different reduction factors and using different integration strategies. It can be observed from this figure that the computations produced using a stress correction scheme recorded horizontal displacements of 0.1 and 0.4 m for the adaptive explicit method and forward Euler method, respectively. In contrast, the same methods without stress correction produced non-convergence with negligible displacement. In addition, there are noticeable difference in the safety factors with and without stress correction (1.26 with correction and 1.31 without correction), which shows the impact of stress correction on the numerical results.

The effect of the stress correction scheme can be further interpreted based on the shear surface of the slope. Figure 15 shows contour plots of the shear surface at the final increment obtained using the RKF23sec method, (I) with and (II) without stress correction. Analysis of this figure reveals that whereas a failure surface, as depicted by the displacement in Fig. 15a, can be observed in the data produced with a stress correction scheme, no failure surface is generated in the computation undertaken without a stress correction scheme, as shown in Fig. 15b. Similarly, Fig. 15c also shows a shear band (as depicted by the equivalent strain) produced in the computation with stress correction, but no such a failure band is generated in Fig. 15d. Note that the above different results are obtained from computation under different safety factors.

5 Conclusions

This paper presents the numerical implementation of a simple hypoplastic constitutive model using finite element method. Various commonly used integration methods, including both implicit and explicit methods, are enhanced by stress correction, with such influence factors as the load increment size and the specified error tolerance on the performance of the different integration strategies studied. The main conclusions of this work can be summarised as follows:

  1. (1)

    The hypoplastic model is characterised by a failure surface and a bound surface, which restricts the accessible stresses. For some loading directions, the stress may surpass the failure surface but is still within the bound surface. bound surface largely increases with an increase in the critical friction angle. However, the difference between the failure surface and the bound surface, in term of mobilized friction angle, is so large that a stress correction is necessary.

  2. (2)

    Compared with the implicit Crank–Nicolson method, the adaptive explicit methods shows better performance in the numerical computations, being less sensitive to the loading direction and increment size since the incremental step size can be changed automatically according to the prescribed accuracy requirement. Moreover, the implicit integration techniques are also sensitive to the specified error tolerance due to the strong nonlinearity of the hypoplastic model. The accuracy can thus be effectively improved by tightening the error tolerance, without increasing the number of substeps nor the CPU time.

  3. (3)

    Although the tested adaptive explicit methods can produce accurate numerical results, the intrinsic “shortcoming” of the hypoplastic model, that some stresses may lie outside the failure surface, cannot be overcome. In addition, the bound surface is often too far from the failure surface. In order to make sure that the stress does not surpass the failure surface, a stress correction must be employed. Such a stress correction can guarantee that the stress does not go beyond the failure surface. Moreover, the stress correction stabilizes the numerical computation for large increment size.