Numerical integration and FE implementation of a hypoplastic constitutive model

Hypoplastic constitutive equation based on nonlinear tensor functions possesses a failure surface but no yield surface. In this paper, we consider the numerical integration and FE implementation of a simple hypoplastic constitutive equation. The accuracy of several integration methods, including implicit and explicit methods, is examined by performing a set of triaxial compression tests. Adaptive explicit schemes show the best performance. In addition, the stress drift away from the failure surface is corrected with a predictor-corrector scheme, which is verified by two boundary value problems, i.e. rigid footing tests and slope stability.


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.

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 _ e r ¼ LðrÞ : _ e À NðrÞk _ ek ð 1Þ where the terms L and N, respectively, denote the linear and nonlinear components, r is the Cauchy stress tensor, and _ e is the stretching tensor. k _ ek ¼ ffiffiffiffiffiffiffiffiffiffi ffi trð _ e 2 Þ p stands for the Euclidean norm. The Jaumann rate of the Cauchy stress tensor r is defined in terms of the time-derivative of the Cauchy stress tensor _ r and the spin tensor w The stretching and spin tensors are related to the velocity gradient tensor through where v is the velocity and r 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 _ e.
in which C i ði ¼ 1; 2; 3; 4Þ are dimensionless parameters. The deviatoric stress tensor r Ã in Eq. (4) is defined by r Ã ¼ r À 1=3ðtrrÞd ij ; with d 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: where e and e crt refer to the current void ratio and critical state void ratio, respectively, and a 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, _ e v ¼ trð _ eÞ, according to the following relationship: The critical state void ratio is calculated according to Li and Wang [13]. A slightly modified form is used in this work: where e co , k, and n are parametric constants, and p ¼ trr=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 r with the following translated stress tensor [33,42]: where the translated scale p t ¼ c=tan/ and the parameters c and / are, respectively, the cohesion and friction angle of the cohesive soil.

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 N contains the critical state function I e , we have: where _ ẽ ¼ _ e=k _ ek stands for the direction of strain, and the symbol denotes an outer product between two tensors. By making use of the fact that _ ẽ : _ ẽ ¼ 1, the failure criterion can be readily derived: Therefore, the condition of invertibility is kL À1 : Nk\1. This condition is identical with the requirement that a stress state r 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: where J 2 and I 1 are, respectively, the second deviatoric stress invariant and the first stress invariant, and 1 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: With the help of Eqs. (11) and (12), the failure and bound surfaces can be plotted once the model parameters are given.

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 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 W 2 ¼ 0: 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: where the variables are expressed as function of the stress invariants I 1 and 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 e N 1 and I e N 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).
Let us consider three stress states r b , r f and r 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: where r 1 and r 3 are the maximum and minimum principal stresses, respectively. The stress ratio can be obtained from the mobilised friction angle / mob . In contrast to the Mohr-Coulomb yield criterion, the mobilised friction angle / mob varies for different Lode angles for the hypoplastic model. Obviously, the mobilised friction angles obtained from r b , r f and r s are dependent on the Lode angle. Figure 2 shows the three mobilised friction angles / mob b ; / mob f and / mob s for different critical friction angles. For the stress state in triaxial compression (the Lode angle is zero) with a 20 frictional angle, one obtains / mob b ¼ 21:4 , / mob f ¼ 20 and / mob s ¼ 11:7 while for the stress state in extension (the Lode angle is 60°) the corresponding friction angles are / mob b ¼ 28:8 , / mob f ¼ 26:3 and / mob s ¼ 16:8 . 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.

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: where n denotes the nth increment of the nonlinear analysis. P n ext is the nodal force vector, and u is the unknown nodal displacement vector. Thus the unbalanced force vector can be expressed as:   For the ith iteration in the framework of the standard Newton-Raphson procedure, the global stiffness matrix reads: For the hypoplastic constitutive model, the Jacobian matrix can be readily written out: For constitutive model (4), the tangential stiffness tensor at a given stress r ij is given below: where I ijmn is the fourth rank identity tensor with compo- 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].

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 2 ½t n ; t nþ1 can be written as: According to Eq. (6), a closed form of integration for the void ratio is available: The solution for Eqs. (22) and (23) can be obtained stepwise according to one of the integration algorithms in the next subsection.

Generalised midpoint algorithms
Following the method proposed in [5,6], all stress components and state variables are collected in vector y for the simultaneous integration of Eqs. (22) and (23): where g i ði ¼ 1. . .mÞ are additional state variables, e.g. void ratio. The integration of Eq. (24) requires the solution of initial value problem: y 0 ðtÞ ¼ HðyðtÞÞ; y 0 ð0Þ ¼ yð0Þ: 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: where Dt nþ1 ¼ t nþ1 À t n is the time step increment and the parameter h is a prescribed constant within [0, 1]. Note that generalised midpoint algorithms with values of h 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 h ¼ 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).

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 DT k 2 ð0; 1 given by the following equation: two different approximate solutions of the evolution problem (25) are obtained simultaneously according to The two function U 1 and U 2 are constructed as follows according to Table 2, which summarises different integration methods: For the sake of simplicity, we usually set and # 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: The integration over the k-th substep is assumed to be successful when, for a given stress error tolerance STOL, with the new substep size then estimated using the following extrapolation formula: If the estimated error is less than the prescribed accuracy tolerance STOL, the step is accepted and we enlarge our step size according to If condition (32) is not satisfied, the k-th substep will be recalculated with a smaller step size DT Ã k : 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 y.

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 ðrÞ [ 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. Let us consider a stress state inside the failure surface, e.g. stress r 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 r trial nþ1 at the ðn þ 1Þth step of analysis will violate the consistent condition, so that 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 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: where g is an unknown multiplier. Using the previous definition of J 2 , it follows that: In order to return the stress state to the failure surface, it is desirable that the total strain increment, De, remains unchanged, since this is consistent with the displacementbased 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: which yields the unknown multiplier g.
On the other hand, if the mean stress p trial nþ1 [ 0, as shown in Fig. 3b, the stress is corrected to the apex of the failure surface 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). The stress correction scheme can also be easily adapted to incorporate the effect of critical state and cohesion. To this end, the constant 1 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.

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 userdefined 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.

Performance of integration methods
Two groups of integration methods are analysed in the numerical triaxial compression tests: (1) Constant substep method: the forward Euler method (FE), the modified Euler method (ME), and the Crank-Nicolson method (CN); (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.
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 ¼ ky n exact À y n k ky n exact k For the triaxial compression tests, an initial isotropic stress state with r 11 ¼ r 22 ¼ r 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. 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.

The effect of stress correction
To evaluate the effect of stress correction, three different initial isotropic stress states with r 11 ¼ r 22 ¼ r 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. 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.
As discussed in Sect. 2, stress drift is particularly relevant for large critical friction angle, e.g. / ¼ 30 . 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 / ¼ 15 to 4.5% for / ¼ 45 . Obviously, stress correction is important in the numerical implementation of the hypoplastic constitutive models.

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,  Fig. 8 a Stress path of triaxial tests, b the value of failure function f in drained triaxial test 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.

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

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 socalled 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  FEs (100 substeps) 23.7 6 9 10 6 10,000 -ME (100 substeps) 41. RKF45sec method STOL = 10 -9 29.5 644,161 1429 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]. 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 / 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, 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. 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 / ¼ 16:7 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 nonconvergence 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.
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   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.

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) 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   displacements (a, b) and strains (c, d) upon slope failure by strength reduction using the RKF23sec method, (I) with and (II) without stress correction friction angle, is so large that a stress correction is necessary.
(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) 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.