Elastoplastic constitutive model for granular soil based on hyperbolic failure surface

Numerical procedure for the development and implementation of a new elastoplastic constitutive model for cohesionless granular materials is presented in the paper. The presented model is based on the hyperbolic failure envelope developed using the theory of incremental plasticity. The governing parameter method (GPM) is used for implicit integration of constitutive relations. The developed algorithm is implemented in the general-purpose finite element program PAK designed for the static, dynamic, linear and non-linear analysis. The model is calibrated and verified through numerical simulation of triaxial test and direct shear test.


Introduction
In numerical analyses of geotechnical problems, it is crucial to choose an appropriate constitutive model that describes mechanical behaviour of analysed material. After adopting the constitutive model, it is essential to determine its parameters. In practical geotechnics, it is convenient to use material parameters that are determined through conventional material testing procedures. Additionally, it is desirable to have a clear physical meaning of constitutive model parameters. In order to include previously defined criteria, a new constitutive model for cohesionless materials, based on hyperbolic failure surface, was developed. The boundary surface of the proposed constitutive model represents surface as the stress state function based on the Mohr-Coulomb failure envelope [1], which results in a more realistic description of granular material mechanical behaviour compared to models with linear boundary surface. The simplest empirical approximation for complex interaction mechanism between particles of different size, form, and mineral composition, expressed by linear dependence, is the stress-strain relation at failure [2]. However, shear resistance is a complex phenomenon in which other mechanisms besides friction are involved [1,3]. There are many proposals for failure envelope description as a function of effective stresses. They can be grouped into parabolic, logarithmic, and hyperbolic expressions [4][5][6][7]. Most of these expressions have disadvantages such as validity in a limited stress range, parameters dependent on stress unit, lack of asymptote of failure envelope, model parameters are devoid of clear physical meaning, etc. Therefore, the failure envelope must satisfy the following conditions: it should be applicable in the entire range of possible stresses and for several types of soil, it should have parameters with clear physical meaning, and it must be consistent with basic concepts of the accepted theory. On the other hand, in the FEM model implementation, stress integration is performed for each integration point of the element, and so it is important for computational algorithm to be efficient. Generally, it is necessary to have a robust algorithm that provides reliable results for all the possible load cases. Experimental investigation of soil shows that the linear failure envelope can be used as an acceptable approximation in a relatively wide range of effective stresses for some types of soil, such as loose sand. Linear failure envelope, such as Mohr-Coulomb and Drucker-Prager [8,9], are often in use because of their simplicity and the relatively easy determination of model parameters. These models have limited possibilities in pre-failure behaviour, especially in cases when the stress level significantly changes, or when the stress path varies considerably [10]. However, the failure envelope of dense materials gives value c on the t axis by extrapolation, which represents a component of a nonfriction material strength. As soil is a set of non-cemented grains, in the absence of normal stress, the soil has no shear strength and cannot receive effective tension stress. In this case, cohesion as the value of shear strength at zero level strain, does not exist, but it is a consequence of the linear approximation of failure envelope. One of the models that eliminate the lack of the linear failure envelope is the Duncan-Chang model [11]. This model is based on the idea that the stress-strain curve in drained triaxial compression tests can be approximated by hyperbola [12], and on the idea that soil stiffness can be formulated as a stress-dependent parameter using a power law formula [13]. Failure of Duncan-Chang model is described by means of the Mohr-Coulomb failure criterion, but this is not properly formulated in the plasticity framework. As a result, this model cannot describe dilatation since this would require a Poisson's ratio greater than 0.5, which is theoretically not admissible [10]. The proposed constitutive model overcomes the abovementioned shortcomings of the most commonly used constitutive models of soil since it properly describes behaviour of granular materials in the case of significant change in normal stress level, and as it has no cohesion as the consequence of failure surface linear approximation. In addition, parameters of the proposed model have clear physical meaning and can easily be obtained from conventional tests of materials. It should be noted that the developed model has no hardening feature nor the previously mentioned models. To calculate the stress, an algorithm for implicit stress integration was developed and implemented using the incremental plasticity theory [14,15] and the governing parameter method [16]. This method represents a generalization of the radial return method used in general plasticity [17], which is based on the fact that the calculation of unknown stresses and internal material variables is reduced to the determination of one (governing) parameter [18][19][20][21]. The return mapping algorithm was used for the implicit stress integration [22][23][24]. The implicit integration method ensures that the yield condition is satisfied at each time step, thereby achieving no deviation from the yield surface. Also, the implicit integration methods allow significantly greater time step than the explicit integration, which leads to faster solution [25]. It should be noted that implicit methods are very often used for solving geotechnical problems [26], as well as other elastoplastic and viscoplastic problems [27].

Constitutive model formulation
The use of constitutive models based on a linear failure envelope is accurate enough for numerical simulation of mechanical behaviour in granular material for higher values of normal stress. However, for lower values of normal stress, the linear failure envelope does not reflect real mechanical behaviour of granular (cohesionless) material. Therefore, the hyperbolic soil model was developed in order to realistically GRAĐEVINAR 72 (2020) 2, 115-125 Elastoplastic constitutive model for granular soil based on hyperbolic failure surface describe mechanical behaviour of granular material for all possible stress states. The failure envelope of this model resulted from the model with the non-linear deformable sawteeth [28], as presented in Figure 1. The failure surface of this constitutive model describes quite realistically mechanical behaviour of granular materials especially for lower values of normal stress. The hyperbolic constitutive model is a modification of the Mohr-Coulomb model in which the internal friction angle is defined as a function of stress state. The hyperbolic failure surface of the model was defined using three material parameters with clear physical meaning, as explained in the paper. Equations of mechanical behaviour of soil as a porous media describe the laws referring to solid skeleton as it is the case for other non-porous materials. The effective stress is used in the analysis of partially or fully saturated media, i.e. in the determination of constitutive relations. Effective values of stress and strain are usually marked with the sign ('). However, in the case of coarse type of soil such as sand and gravel soil, the pore pressure dissipation takes place very rapidly, and so the effective stress and total stress are equal, which is why this sign will be omitted in this paper. Starting from the fact that there is no cohesion in granular unbound materials [29,30], shear strength of the soil can be defined using the following equation (1).
where maximal value of internal friction angle is defined as (2) and f cv represents the angle of shear resistance with constant volume, whereas ψ represents material dilatation. The angle of internal friction of a material is a function of normal effective stress and can be defined as (3) According to [1], the second part of the eqn. (3) is formulated as (4) and so, the internal friction angle of a material (3) assumes the following form (5) Using eqn. (5), the shear strength of material (1) can be finally formulated as (6) which represents the equation of the failure shear stress envelope [31] as a function of normal stress.
In the eqns. (3) to (6), the following material parameters are used: The failure shear stress envelope is defined with eqn. (6) and presented in Figure 2. As can be seen from this figure, when normal stress reaches zero (σ n → 0),the slope of the failure envelope tangent in the origin (initial slope) is  For the stress state at failure, described using Mohr's circles, the tangent to the Mohr's circle from the origin can be defined (Figure 3), and its angle is: (8) where as the corresponding normal stress is (9)

Figure 3. Failure envelope of the model and conversion of material parameters
It is obvious that, due to curvature of failure envelope, point K where the tangent touches Mohr's circle and forms the angle f S with axis σ, does not match point F, where the tangent is cutting the failure envelope. However, according to [32], this deviation is small enough and can be neglected, and so it is adopted that p N ≈ p F . The change of the Mohr's circle tangent angle can be formulated as a function of normal stress σ ff , and so, using (5), the equation of internal friction angle can be formulated as (10) The hyperbolic failure envelope defined in stress space σ-t can be transformed into the elastoplastic constitutive model in principal stress space if a generalization is made so that, instead of σ ff stress, the mean stress is used (11) Therefore, the internal friction angle (5) can be calculated as (12) where, parameter p AV can be calculated using the following equation It should be emphasized that the failure envelope defined using internal friction angles (5), (10) and (12) does not match, except when the normal stress reaches zero or infinity or when the normal stress matches the central angle f M . According to [26], the difference in friction angle has the range of ±0.2° and reduces with the decrease in angle difference Df. Therefore, the introduced approximation can be used in numerical analyses with sufficient accuracy for finding solutions to engineering problems.

Yield surface
The analysis of the Mohr-Coulomb and hyperbolic soil model failure envelope leads to the conclusion that eqn. (1) (14) The yield surface of the hyperbolic soil model, defined using eqn. (14), is shown in Figure 4 in principal stress space. This constitutive model does not contain a hardening feature. In other words, this model is elastic-perfectly plastic. Therefore, the failure surface of the model represents, at the same time, the yield surface. The term failure surface is used in this paper, but it should be noted that the failure surface represents the yield surface as well. This paper discusses the associative yield condition or the case when the failure surface and the plastic potential surface match (g=f). So, these two surfaces will be described using the same equation, which does not represent a general case. GRAĐEVINAR 72 (2020) 2, 115-125 Elastoplastic constitutive model for granular soil based on hyperbolic failure surface In eqn. (14), I 1 stands for the first stress invariant, J 2D is the second deviatoric stress invariant, q is the Lode's angle, whereas f (σ n ) represents the internal friction angle defined by eqn. (12). Required parameters of the presented constitutive model for granular soil based on hyperbolic failure surface are summarized in Table 1.

Derivatives of yield function
The equation of yield surface in the hyperbolic soil model (14) represents a composite stress function, and so its derivative can be calculated using chain rules [33] according to (15) where σ is the stress tensor for Cartesian components which, in case of isotropic materials, contain six components and can be written in vector form (16) Individual parts of eqn. (15) represent derivatives of yield function (14) with respect to stress invariants, Lode's angle and internal friction angle. Individual derivatives can be calculated as Derivatives of the first stress invariant and the second deviatoric stress invariant are: The derivative of the internal friction angle with respect to stress is (23) where as the derivative of Lode's angle can be calculated using chain rule [33] as follows For a shorter form, the following change is introduced in eqn.
The above equations are used in implicit stress integration of the hyperbolic soil model. Steps for stress integration using this model are summarized in form of the algorithm presented in the following section.
A flaw of this method is its potential inability to achieve the numerical solution convergence in local iterations for higher increment values as well as in the cases when the yield function or the plastic potential function is highly nonlinear [45][46][47][48].
Using the developed theory, the trial elastic solution or so-called elastic predictor is calculated at the beginning of each time step. The yield condition is then checked. If plastic strain occurs in the current time step, the yield condition will not be fulfilled and so it will be necessary to correct total strains by calculating the plastic corrector. The plastic corrector represents the part of plastic strains in total strain (increment). In order to calculate the accurate value of plastic strain increment in the time step, local iterations are performed with the correction of the scalar, or the intensity of plastic strain increment is corrected whereas the direction of plastic strain vector is changed. In some constitutive models, derivatives of yield function and plastic potential can be very complex, which complicates calculation of these derivatives. However, this can be overcome by using numerical derivatives of yield and plastic potential function instead of analytical derivatives. The developed algorithm for implicit stress integration of the hyperbolic soil model is presented in Table 2. The presented algorithm is incorporated in the PAK program [49] and verified via that program through numerical examples.

Model validation and algorithm verification
The developed algorithm for implicit stress integration using the hyperbolic soil model was verified on two test examples. The first verification example is a numerical simulation of triaxial test aimed at verifying whether the developed model accurately describes the strength of the material sample for given material parameters. The second example represents a numerical simulation of the direct shear test and its purpose is to verify whether the developed model accurately describes mechanical behaviour of real samples during shear load. Material parameters were identified using back analysis [50]. Numerical simulation results were compared with analytical and experimental results.

Numerical simulation of triaxial test
Numerical modelling of triaxial test is a simple way to check whether the developed constitutive model describes the strength of the material in accordance with the theoretical failure criterion for given model parameters. Therefore, numerical solution of triaxial test was not compared with experimental results. The stress path can generally be classified according to the type of loading and direction of loading. In this case, the performance  Material parameters used in numerical simulation of triaxial test were taken from report [51]. Material parameters of the model are shown in Table 3.

Table 3. Material parameters used in triaxial test simulation
Triaxial test simulation results, for different values of confining pressure in case of compression and extension are shown in Figure 7.
Numerical results are shown in the form of stress paths for both analysed cases and for all confining stresses in the σ m -q stress space, where σ m is the mean stress (confining stress) while q is the second invariant of deviatoric stress shown in (29).

Figure 7. Compression and extension stress path in triaxial test simulation
According to the analysis of results, it can be concluded that the developed stresses in the model using this method coincide with the failure surface obtained from the hyperbolic soil model. In other words, the developed model describes the strength of the material that corresponds to theoretical values of failure stress.

Numerical simulation of direct shear test
The next verification example of the developed algorithm for implicit stress integration of hyperbolic soil model is a numerical simulation of the direct shear test. This relatively simple test of the developed elastoplastic model for granular soil based on hyperbolic failure surface was checked for axial compression and axial extension. Four different confining stresses were used. In the compression test, the sample was subjected to load using the hydrostatic stress state, after which the stress was increased in one direction, while it remained constant in the other two directions. In the tension test, after set the hydrostatic stress state, the stress was reduced in one direction, while it remained constant in the other two directions. The FE model used consists of one 3D hexagonal finite element with unit dimensions. The model geometry, boundary conditions, and loads are shown in Figure 5. The analysed model has three planes of symmetry and so appropriate boundary conditions of symmetry were used. Model loads were applied using the prescribed pressure in three coordinate directions.
The load was increased in multiple steps until sample failure (inability to achieve convergence of numerical solutions).

Figure 5. FE model for triaxial test simulation
In order to confirm that the model provides analytical stress values in failure for different stress states, the procedure was repeated for four levels of confining stress: σ m = 0.213 MPa, 0.421 MPa, 0.839 MPa, and 1.665 MPa. In the compression test, after reaching the initial stress state, the vertical pressure was increased until failure. In the extension test, after reaching the initial stress state, the vertical pressure was reduced until failure. Load function useds are shown in Figure 6. is often used for identification of material parameters, and so its numerical simulation is convenient for validation of the constitutive model as well. The analysed material is the rockfill from downstream slope of a dam (granular material), and so the application of the hyperbolic soli model is suitable for numerical simulation of mechanical behaviour of this material. Experimental results of the direct shear test of the supporting body material from an embankment dam [51] are presented in Table 3 and used for identification of constitutive model parameters. The same normal stress values were used in numerical simulation of the direct shear test. A large-scale graphic representation of direct shear test results, and the failure surface obtained by identification of parameters, are shown on Figure 8. Estimated material parameters from the hyperbolic soil model are presented in Table 5. These estimated parameters were used in numerical simulation of shear test. The FE model used consists of one finite element of unite dimensions with boundary conditions and loads. It is presented in Figure 9. Boundary conditions used in numerical simulation correspond to boundary conditions that exist in the shear layer of the specimen.  (Table 4). . Load functions used in the test device were the same as the load function used in numerical simulation (Figure 9b).

Figure 9. FE model for direct shear test simulation and load functions
Numerical simulation results and test results are represented in the form of t xy -e x ( Figure 10). By comparing numerical results obtained using the developed algorithm with experimental results shown in Figure 10, it can be noticed that the developed constitutive model significantly follows the trend of experimental results. Significant deviations can be observed for lower strain values. This may be due to the fact that teh developed model does not contain hardening feature, which could be dealt with in the scope of further development of the model.

Strip footing
This example represents numerical simulation of a strip footing resting on cohesionless soil taken from the literature [53].
In order to calibrate the model, th eload-settlement curve of laboratory scale experiment on strip footing was obtained. Numerical simulation was performed using the Mohr-Coulomb constitutive model, as well as the soil model based on the hyperbolic failure surface. The obtained numerical solution was compared to the laboratory testing solution from the same literature.
A schematic representation of sample testing in laboratory, with boundary conditions and loads, is shown in Figure 11, while Figure 12 shows the model of finite elements used in the footing test numerical simulation.  Figure 12 also shows boundary conditions and loads used in numerical simulation, according to boundary conditions and loads used in the laboratory test. Due to the symmetry of the problem, a half of the specimen is modelled using an appropriate boundary condition. The model loading was conducted in two stages: in the initial stage a dead weight was assigned in order to establish the initial stress state, while in the second stage the footing load was applied, using prescribed displacement. As previously mentioned, material parameters of the Mohr-Coulomb model and the hyperbolic soil model were estimated using experimental results by nonlinear curve fitting, as shown in Table 6.  Numerical simulation results are shown in Figure 13. It can be seen from the figure that the results obtained using both constitutive models are in accordance with laboratory test results.  However, it can be seen that, for the estimated parameters, the results obtained by using the soil model based on hyperbolic failure surface better correspond to the experimental results of the footing test. This fact shows that the soil model based on hyperbolic failure surface describes more accurately mechanical behaviour of non-cohesive granular materials such as dry sand.

Conclusion
The development of constitutive model for cohesionless granular material based on the hyperbolic failure surface, using incremental plasticity theory, is discussed in the paper. The yield (boundary) surface of the model is defined by modifying the Mohr-Coulomb yield surface and by introducing a variable internal friction angle as a function of stress state. The model describes a more realistic mechanical behaviour of cohesionless soil, especially for lower values of normal stress. The model formulation is given and the constitutive relation development for implicit stress integration is presented in detail. The yield surface of the model is defined using three material parameters whose physical meaning is presented in the paper. These material parameters can be obtained using either the direct shear test or triaxial test. A return mapping algorithm is applied to implement the model in the generalpurpose element method program called PAK. The algorithm is verified through several test examples. The developed model ensures good correspondence between numerical results and analytical results and significantly follows the trend of the experimental results. Some deviations can however be observed for lower values of strain, which is due to the fact that the developed model does not have a hardening feature. These deviations could be studied in the scope of further development of the model. The above properties of the model confirm its applicability in real-life geotechnical problem solving. The suitability of the hyperbolic constitutive model based on nonlinear failure envelope for problem solving in a variety of engineering applications is reflected in the fact that parameters can be obtained directly using standard laboratory tests. The developed model can be improved by introducing a non-associated yield condition. In addition, the model can be modified by introducing kinematic hardening, and so the model can be suitable for dynamic analysis of granular materials. Due to the simplicity of reducing the shear stress envelope, this constitutive model is suitable for application of the shear strength reduction (SSR) method.