Acta Polytechnica CTU Proceedings https://doi.org/10.14311/APP.2024.49.0097 Acta Polytechnica CTU Proceedings 49:97–103, 2024 © 2024 The Author(s). Licensed under a CC-BY 4.0 licence Published by the Czech Technical University in Prague FORMULATION AND IMPLEMENTATION OF NONLINEAR ELASTICITY IN ADVANCED CONSTITUTIVE MODELS IN GEOMECHANICS Tereza Žalskáa,∗, Michal Šejnohab, Tomáš Jandab, Alena Zemanováa a Czech Technical University in Prague, Faculty of Civil Engineering, Department of Geotechnics, Thákurova 7, 166 29 Prague, Czech Republic b Czech Technical University in Prague, Faculty of Civil Engineering, Department of Mechanics, Thákurova 7, 166 29 Prague, Czech Republic ∗ corresponding author: tereza.zalska@fsv.cvut.cz Abstract. Several formulations of a nonlinear elastic behavior within the yield surface of the selected advanced constitutive models are described and compared. Since concentrating on geotechnical applications the accompanied numerical simulations are limited to basic laboratory tests such as isotropic compression, oedometer, and drained triaxial compression. The results show a significant influence of the size of the initial load step on final predictions particularly when starting from a very low initial stiffness associated with the assumed zero initial stress. Differences in the predicted response arising from different formulations are also discussed. These might be quite significant and the design engineer should be aware of that when choosing a particular computational software. Keywords: Hardening soil model, generalized cam clay model, soft soil model, nonlinear elasticity, oedometer, isotropic compression, triaxial compression. 1. Introduction Constructing underground structures in densely popu- lated areas, preventing large scale landslides, as well as being able to accurately simulate a complex response of soils observed experimentally, opens the door to the application of advanced constitutive models [1, 2]. Grounding on experimental evidence such models at- tempt to address a nonlinear response of soils already at initial stages of loading. An illustrative example are the critical state models [1, 3] where the nonlinear stress-strain behavior within the yield surface arises naturally. A large group of formulations introduce such a behavior in the spirit of Duncan and Chang hyperbolic model [4–7]. While their appearance in available commercial software is common [8–10] (to cite a few), the details on their numerical implementa- tion within a given constitutive model may differ. This in turn may result in predictions, which are consider- ably influenced by the initial load increment. In some formulations, the load step dependence yields loading- unloading curves which deviate even for loading condi- tions that do not exceed the yield limit, often resulting in an artificial residual strain. It is therefore crucial to understand such potential drawbacks and the user should be aware of that when creating the computa- tional model for a given geotechnical construction. The present paper revisits this issue in light of three popular constitutive models implemented in the above mentioned commercial codes. In particular, the Hard- ening Soil Model [5, 6], the Generalized Cam clay model [10, 11], and the Soft soil model [12] are exam- ined in Section 2 with emphasis on the implementation of the associated nonlinear elastic constitutive model. The results of numerical simulations together with a thorough discussion are presented in Section 3. All calculations are carried out using the GEO5 FEM software [10]. 2. Advanced constitutive models The present section provides details on the formulation of nonlinear elasticity of the selected constitutive mod- els, while the plasticity issues will be presented only briefly with details available in the above mentioned references. 2.1. Hardening soil model The Hardening soil (HS) model developed in [5] and later extended in [6] to account for small strain yield- ing is plotted in Figure 1. As seen in Figure 1a it com- bines two yield surfaces where the shear part, bounded by a limiting yield surface of the Mohr-Coulomb type, is derived on the assumption of the hyperbolic stress- strain response in triaxial compression as depicted in Figure 1c with ε1 being the major principal strain and q representing an equivalent deviatoric stress mea- sure. In the formulation implemented in GEO5 FEM, both the shear yield surface and the cap yield surface, which is assumed elliptic and centered at the origin of the meridian plane, plot as a smooth hexagon in the deviatoric plane in the spirit of the Matsuoka-Nakai (MN) yield surface [13], see Figure 1b. The evolution (hardening) of both surfaces is controlled by the cur- rent value of the equivalent deviatoric plastic strain γ and preconsolidation pressure pc, respectively. 97 https://doi.org/10.14311/APP.2024.49.0097 https://creativecommons.org/licenses/by/4.0/ https://www.cvut.cz/en T. Žalská, M. Šejnoha, T. Janda, A. Zemanová Acta Polytechnica CTU Proceedings (a). Yield surface with zero cohesion in merid- ian plane. (b). Yield surface in deviatoric plane. (c). Response in triaxial compres- sion. Figure 1. Graphical representation of Hardening soil model. (a). Yield surface in meridian plane. (b). Behavior under isotropic compression. (c). Elastic and plastic part of total strain increment. Figure 2. Graphical representation of Generalized Cam clay model. The nonlinear elastic model due to Vermeer [14] writes the elastic modulus Eur in the form of a power law as a function of the current value of mean effective stress σm as: Eur(σm) = Eref ur ( σm − c cot φ σref m − c cot φ )m , (1) where Eref ur , σref m , m are model parameters, c, φ are the shear strength parameters, the cohesion and the angle of internal friction, respectively. In particular, Eref ur is the reference elastic modulus for a reference stress σref m typically equal to −100 kPa. Note that the standard elasticity sign convention is adopted with tension being positive. For simplicity, we shall now consider the volumetric response only and assume a linear variation of the mean stress σm over a given load increment as: σm = ησi m + (1 − η)σi+1 m . (2) Setting η = 1 gives fully explicit formulation (forward Euler (FE) method, σm = σi m in Equation (1)) while η = 0 renders fully implicit (backward Euler (BE) method, σm = σi+1 m in Equation (1)). Henceforth, the former formulation will be referred to as the constant elasticity (CE) return while the latter formulation will be termed the variable elasticity (VE) return. Fur- ther distinction is made by considering either secant formulation in the form: σi+1 m = Kur(σm(η))(εi v + ∆εv), (3) or incremental formulation provided by: σi+1 m = σi m + Kur(σm(η))∆εv, (4) where Kur is the elastic bulk modulus given by: Kur = Eur 3(1 − 2νur) . (5) In particular, cases with η = 1 (FE), η = 0 (BE), and so called mid-point rule (MP), i.e., η = 0.5 → σm(η) = 0.5(σi m+σi+1 m ), will be examined in Section 3. 2.2. Generalized Cam clay model The Generalized Cam clay (GCC) model has been introduced to reconcile the principal drawback of the Modified Cam clay (MCC) model associated with overestimating the shear strength and consequently an excessive softening of highly overconsolidated soils. To this end, the supercritical (dilation) part of the MCC model was suitably adjusted, see Figure 2a, resulting in the dependence on both the critical friction angle φcs and peak friction angle φ via parameter β given by: β = ( sin φ sin φcs 3 − sin φcs 3 − sin φ )2 . (6) 98 vol. 49/2024 Nonlinear elasticity in advanced constitutive models Similar to the HS model, the MN yield surface repre- sents the projection into a deviatoric plane. The nonlinear elastic constitutive model derives from the assumed response in isotropic compression shown in Figures 2b and 2c. When moving along the unloading-reloading κ-line the rate form of the stress-strain law is provided by: σ̇m = −σm κ∗ ε̇el v , κ∗ = κ 1 + e , (7) where κ is called the swelling modulus, e is the void ratio. Integrating over a given time increment gives the evolution of the mean effective stress in the form: σi+1 m = σi m exp [ −∆εel v κ∗ ] . (8) Next, writing the stress increment ∆σm as: ∆σm = σi+1 m − σi m = σi m ( exp [ −∆εel v κ∗ ] − 1 ) = Ks∆εel v , (9) allows us to express the secant bulk modulus Ks as: Ks = σi m exp [ −∆εel v κ∗ ] − 1 ∆εel v . (10) The analysis may simplify if leaving the concept of variable elasticity return represented by Equation (10) and assume the bulk modulus: Ki = −σi m κ∗ , (11) to be constant over the load increment (constant elas- ticity return), recall Equation (7). 2.3. Soft soil model In the original MCC model the shear strength depends solely on the slope of the critical state line g which for the triaxial compression reads: g = 2 √ 3 sin φcs 3 − sin φcs . (12) Brinkgreve showed in [12] that it is not possible to simulate both the oedometer stress path and full tri- axial compression path with good accuracy using the same value of g. This let to the formulation of a Soft soil (SS) model combining the cap yield surface in the form of MCC model to predict a correct K0-path and the shear yield surface, here again assumed in the form of MN model, to ensure a correct shear strength. Such a yield surface is plotted in Figure 3 where Mc for triaxial compression is given in terms of the coeffi- cient of lateral earth pressure of normally consolidated ∙ Figure 3. Yield surface of Soft soil model in meridian plane. soils KNC 0 [12] whereas Ms follows from the solution of nonlinear MN yield function [13]. While the shear yield surface represents a response of an elastic perfectly plastic material, the cap yield surface may undergo isotropic hardening. To avoid potential softening within the cap model, which is not allowed in the present implementation, the Mc-value should be larger than the slope of the Matsuoka-Nakai failure surface Ms. Considering the nonlinear elastic response, we fol- low the implementation in [8] and limit attention to the constant elasticity return of Equation (11). As suggested in [12] it is possible to relate Equation (11) to Equation (1) by writing: K(σm) = Kref ( σm σref m )m . (13) When setting m = 1, the logarithmic compression law is recovered with: 1 κ∗ = −Kref σref m . (14) 3. Simulation of simple laboratory tests Numerical simulation of three basic laboratory tests is performed in this section to examine and compare individual formulations presented in the previous sec- tion. In all cases, the computational model consists of two constant strain triangular elements. While the oedometric test is run in plane strain regime, the axisymmetric state of stress is assumed for isotropic and triaxial compression. The loading and boundary conditions are displayed in Figure 4 for individual cases. All simulations assumed zero initial stress. 3.1. Hardening soil model The response in isotropic compression and oedometer is simulated to compare the secant and incremental formulation represented by Equations (3) and (4). As for Equation (2), only the cases with η = 0 (BE) and η = 0.5 (MP) are tested. Although purely elastic behavior is addressed, the nonzero values of c = 10 kPa and φ = 30° are adopted together with Eref ur = 30 MPa 99 T. Žalská, M. Šejnoha, T. Janda, A. Zemanová Acta Polytechnica CTU Proceedings (a). Oedometer. (b). Isotropic compression. (c). Drained triaxial com- pression. Figure 4. Computational models of simple laboratory tests. 0 0,01 0,02 0,03 0,04 0,05 0 2 4 6 8 10 Step 0.0025 Step 0.001 Step 0.0001 Secant approach Incremental approach -ε 3 [-] -σ m [M P a] M ea n st re ss Vertical strain (a). Comparing predictions provided by secant and incre- mental (BE) formulations. 0 0,01 0,02 0,03 0,04 0,05 0 2 4 6 8 10 Step 0.005 MP Step 0.001 MP Step 0.0001 MP Step 0.003 BE Step 0.001 BE Step 0.0001 BE -ε 3 [-] -σ m [M P a] M ea n st re ss Vertical strain (b). Comparing predictions provided by MP and BE implementation of incremental formulation. Figure 5. Response of HS model in isotropic compression. and σref m = −100 kPa to set the initial stiffness of Eur(σm = 0) according to Equation (1). The results for isotropic compression appear in Fig- ure 5. Point out that achieving convergence of the underlying nonlinear equation with σm(η) provided by Equations (2)–(4): R = σi+1 m − σm(η), (15) required a relatively small initial load increment. So only a minor difference in the nonlinear elastic re- sponse is observed for incremental formulation for the tested load steps. On the other, the secant formula- tion shows response, see Figure 5a, which is step size independent. This conclusion is further supported via results presented in Figure 6a. Point out that both formulations adopt in the present study the same model parameters in the calculation of Eur. However, it is clear from Figure 5a that arriving at identical predictions would require adjusting at minimum the power law coefficient m for individual formulations based on the experimental measurements. Figure 5b compares the results associated with MP and BE implementation of Equation (1) in the frame- work of incremental formulation. It appears that MP approach is less sensitive to the size of the loading step. This approach is further promoted by predicting an identical stress-strain curve in both loading and un- loading if keeping the same load increment. Note that such an approach is implemented in [8] and currently also in [10]. The above discussion is fully supported by simu- lating an oedometric test as seen in Figure 6. The influence of the load step size, particularly for BE approach, is clearly evident. 3.2. Comparing generalized Cam clay and soft soil models This section demonstrates potential differences be- tween the variable elasticity approach, adopted with GCC model, and the constant elasticity approach, adopted with SS model, in light of incremental formu- lation. The material data were set such as to arrive at identical predictions in isotropic compression when enforcing the constant elasticity return also with GCC model. The material parameters used in all subse- quent simulations are: κ∗ = 0.0094, λ∗ = 0.0678, ν = 0.26 and g = Mc = 1.074. The initial value of K0 = 1 κ∗ was considered. Before proceeding with our discussion on individual results it is worth mentioning that the GEO5 FEM 100 vol. 49/2024 Nonlinear elasticity in advanced constitutive models 0 0,0025 0,005 0,0075 0,01 0,0125 0,015 0 100 200 300 400 500 Step 0.1 Step 0.01 Step 0.001 Secant approach Incremental approach -ε 3 [-] -σ 3 [k P a] V er tic al s tr es s Vertical strain (a). Comparing predictions provided by secant and incre- mental (BE) formulations. 0 0,004 0,008 0,012 0,016 0 100 200 300 400 500 Step 0.1 MP Step 0.01 MP Step 0.001 MP Step 0.1 BE Step 0.01 BE Step 0.001 BE -ε 3 [-] -σ 3 [k P a] V er tic al s tr es s Vertical strain (b). Comparing predictions provided by MP and BE im- plementation of incremental formulation. Figure 6. Response of HS model in oedometer. software allows for gradually increasing the initial size of the load step with a number of iterations needed to achieve global equilibrium. In the present study, the elastic simulations essentially adopted two times larger step size with every new load step when this option, denoted here as VS, was exploited. In this case, the parameter step0, see e.g. Figure 7, represents just the initial step being gradually increased in the course of analysis. All other calculations considered constant step size throughout the analysis. Exposition to the derived results begins with Fig- ure 7 to show the importance of a relatively small ini- tial load step size to avoid excessive strains attributed to initially very low stiffness. The influence of stiffness variation over the load increment (VE), Equation (10), in comparison to CE return, Equation (11), appears for isotropic compres- sion in Figure 8. Note that the ∗-lines require 10 steps only in comparison to 100 and 1 000 steps associated with red and blue simulations, respectively. Promoting VE return over CE return might seem justifiable when leaving purely elastic response and allowing for plasticity as shown in Figure 8b where simulations based on variable step size match those with the constant step size in the case of VE rather well. Similar conclusions can be drawn from the results in Figure 9 derived for oedometer. Figure 9a further shows that, unlike the HS model, the assumed stiff- ness evolution leads to artificial residual strains upon complete unloading even in the case of pure elasticity, which may become significant for larger load steps. On the other hand, this issue will become less impor- tant when loading the soil beyond its relatively low elastic limit, because the generated permanent plastic strains will then well exceed their elastic counterparts. The last simulations are concerned with the plastic response predicted by the GCC and SS models for the drained triaxial compression test. As seen in Fig- ure 10 the analysis begins with isotropic compression, Figure 4b, when the specimen is first loaded to acquire 0 0,025 0,05 0,075 0,1 0,125 0,15 0 2 4 6 8 10 step = 0.1 step = 0.01 step = 0.001 step 0 = 0.001 (VS) K 0 =1/κ∗ pc 0 = 1kPa VS - variable step size -ε 3 [-] -σ 3 [M P a] V er tic al s tr es s Vertical strain Figure 7. Nonlinear elastic response of GCC model in isotropic compression: influence of initial step size. an initial pressure of 50 kPa. Clearly, no deviatoric stresses are generated during this loading stage. The evolution of deviatoric stresses begins with the second stage of loading corresponding to triaxial compression, Figure 4c. Note that the deviatoric stress measure J is defined as the square root of the 2nd invariant of deviatoric stresses sij as: J = √ 1 2sijsij . (16) Similarly to previous studies the CE predictions linked to variable step size (VS) slightly deviate from those obtained with the constant step size. This is not the case for the approach based on VE return. However, it is fair to mention, that smooth conver- gence with VE required a slightly larger number of iterations, thus a smaller load step size, in contrast to CE return (star and diamond symbols identify individ- ual load steps). Also note that for sufficiently small load steps the difference between the two approaches becomes insignificant. 101 T. Žalská, M. Šejnoha, T. Janda, A. Zemanová Acta Polytechnica CTU Proceedings 0 0,005 0,01 0,015 0,02 0,025 -ε 3 [-] 0 2 4 6 8 10 -σ 3 [M P a] step = 0.01 step = 0.001 step 0 = 0.001 (VS) K 0 =1/κ∗ pc 0 = 1kPa gcc model ss model VS - variable step sizeV er tic al s tr es s Vertical strain (a). Nonlinear elastic response. 0 0,02 0,04 0,06 0,08 0,1 0 2 4 6 8 10 step = 0.001 step = 0.001 (VS) K 0 =1/κ∗ pc 0 = 1kPa gcc model ss model VS - variable step size -ε 3 [-] -σ 3 [M P a] V er tic al s tr es s Vertical strain (b). Plastic response. Figure 8. Comparing GCC and SS models in isotropic compression. 0 0,02 0,04 0,06 0,08 0,1 0 100 200 300 400 500 Loading (500 steps) Unloading (500 steps) Loading (VS, 10 steps) Unloading (VS, 9 steps) gcc model step in = 0.002 VS - variable step size K 0 = 1/κ∗ ss model -ε 3 [-] -σ 3 [k P a] V er tic al s tr es s Vertical strain (a). Nonlinear elastic response. 0 0,1 0,2 0,3 0,4 0,5 0 100 200 300 400 500 gcc model (500 steps) gcc model (VS, 11 steps) ss model (500 stes) ss model (VS, 9 steps) K 0 = 1/κ∗ VS - variable step size) -ε 3 [-] -σ 3 [k P a] V er tic al s tr es s Vertical strain (b). Plastic response. Figure 9. Comparing GCC and SS models in oedometer. 0 0,05 0,1 0,15 0,2 0,25 0 15 30 45 60 75 90 S tr es s J, -σ m [k P a] gcc ss gcc (VS) ss (VS) VS - variable step size σ m J step in = 0.01 -ε 3 [-]Vertical strain Figure 10. Comparing GCC and SS models in drained triaxial compression. 4. Conclusion The paper presented several illustrative examples via simulations of simple laboratory tests to show dif- ferences in the implementation of nonlinear elastic constitutive laws arising in the formulation of three particular advanced plasticity models. It has been confirmed that the chosen method of implementation may considerably affect the predicted response. The user of a given commercial software should be aware of that, since using the same model in different softwares does not assure the same results. Common to all examined models is their dependence on the mean effective stress which may generate a very low initial stiffness when starting from zero initial stress. This in turn leads to excessive strains providing the initial load steps are not sufficiently small. To avoid a large number of loading steps to reach the prescribed load the program GEO5 FEM allows for a gradual increase of the initially small step size in dependence on the number of iterations needed to converge for the previous load increment. It has also been observed that with the GCC and SS models irreversible elastic strains occur as the stiff- ness at the end of the loading and at the beginning of unloading is different. On the other hand, this issue does not seem to play a major role as the influ- ence of irreversible elastic strains become negligible in plastic analyses. In this regard, the approach based on variable elasticity return provides response which 102 vol. 49/2024 Nonlinear elasticity in advanced constitutive models seems less dependent on increasing the load step in simulations adopting the variable step size (VS). To conclude we remind that the GEO5 FEM soft- ware, adopts, similar to [8], the incremental formu- lation with MP approach for the HS model and the concept of constant elasticity return for the SS model. On the other hand, the MCC and GCC models employ the variable elasticity return approach. Acknowledgements The support provided by the SGS project No. SGS23/032/OHK1/1T/11 and by the GAČR grant No. 22-12178S is gratefully acknowledged. References [1] K. H. Roscoe, J. B. Burland. On the generalised stress-strain behaviour of ‘wet’ clay. In J. Heyman, F. A. Lechie (eds.), Engineering plasticity, pp. 535–609. Cambridge University Press, 1968. [2] D. Mašín. A hypoplastic constitutitve models for clays. International Journal for Numerical and Anlytical Methods in Geomechanics 29(4):311–336, 2005. https://doi.org/10.1002/nag.416 [3] D. M. Potts, L. Zdravkovič. Finite element analysis in geotechnical engineering: Volume one – theory. Thomas Telford Publishing, UK, 1999. https://doi.org/10.1680/feaiget.27534 [4] J. M. Duncan, C.-Y. Chang. Nonlinear analysis of stress and strain in soils. Journal of the Soil Mechanics and Foundations Division 96(5):1629–1653, 1970. https://doi.org/10.1061/JSFEAQ.0001458 [5] T. Schanz, P. A. Vermeer, P. G. Bonnier. The hardening soil model: Formulation and verification. In R. B. J. Brinkgreve (ed.), Beyond 2000 in Computational Geotechnics, pp. 281–290. 1999. https://doi.org/10.1201/9781315138206-27 [6] T. Benz. Small-strain stiffness of soils and its numerical consequences. Ph.D. thesis, University of Stuttgart, 2007. [7] J. Veselý. The use of the advanced material models for the numerical modelling of the underground structures in clays. Ph.D. thesis, CTU in Prague, Faculty of Civil Engineering, 2016. [8] Seequent.com. PLAXIS. [2023-11-03]. http://www.plaxis.nl/ [9] ZSOIL. [2023-11-03]. http://www.zsoil.com/ [10] Fine.cz. GEO5 FEM. [2023-11-03]. http://www.fine.cz/ [11] T. Janda, M. Šejnoha. Formulation of generalized Cam clay model. Engineering Mechanics 13(5):367–384, 2006. [12] R. B. J. Brinkgreve. Geomaterial model and numerical analysis of softening. Ph.D. thesis, Delft University of Technology, 1994. [13] J. P. Bardet. Lode dependences for isotropic pressure-sensitive elastoplastic materials. Journal of Applied Mechanics 57(3):498–506, 1990. https://doi.org/10.1115/1.2897051 [14] P. A. Vermerr. Formulation and analysis of sand deformation problems. Ph.D. thesis, Delft University of Technology, 1980. 103 https://doi.org/10.1002/nag.416 https://doi.org/10.1680/feaiget.27534 https://doi.org/10.1061/JSFEAQ.0001458 https://doi.org/10.1201/9781315138206-27 http://www.plaxis.nl/ http://www.zsoil.com/ http://www.fine.cz/ https://doi.org/10.1115/1.2897051 Acta Polytechnica CTU Proceedings 49:97–103, 2024 1 Introduction 2 Advanced constitutive models 2.1 Hardening soil model 2.2 Generalized Cam clay model 2.3 Soft soil model 3 Simulation of simple laboratory tests 3.1 Hardening soil model 3.2 Comparing generalized Cam clay and soft soil models 4 Conclusion Acknowledgements References