Acta Polytechnica doi:10.14311/AP.2017.57.0058 Acta Polytechnica 57(1):58–70, 2017 © Czech Technical University in Prague, 2017 available online at http://ojs.cvut.cz/ojs/index.php/ap NUMERICAL MODELLING OF THE SOIL BEHAVIOUR BY USING NEWLY DEVELOPED ADVANCED MATERIAL MODEL Jan Vesely Department of Geotechnics, Faculty of Civil Engineering, Czech Technical University in Prague, Thakurova 7, Prague, Czech Republic correspondence: veselyjan@centrum.cz Abstract. This paper describes a theoretical background, implementation and validation of the newly developed Jardine plastic hardening-softening model (JPHS model), which can be used for numerical modelling of the soils behaviour. Although the JPHS model is based on the elasto-plastic theory, like the Mohr-Coulomb model that is widely used in geotechnics, it contains some improvements, which removes the main disadvantages of the MC model. The presented model is coupled with an isotopically hardening and softening law, non-linear elastic stress-strain law, non-associated elasto-plastic material description and a cap yield surface. The validation of the model is done by comparing the numerical results with real measured data from the laboratory tests and by testing of the model on the real project of the tunnel excavation. The 3D numerical analysis is performed and the comparison between the JPHS, Mohr-Coulomb, Modified Cam-Clay, Hardening small strain model and monitoring in-situ data is done. Keywords: numerical modelling; advanced material model; soils; small strain stiffness; tunnel excavation. 1. Introduction In the last few years, the numerical analyses are in- creasingly used in the design of the underground struc- tures [1–8]. The designers can work with different user- friendly commercial software and relatively easy simu- late the behaviour of the soil mass. Unfortunately, one important fact is overlooked. Numerical analyses are only approximation of real behaviour and are highly dependent on the input parameters [9, 10]. One area, which can affect the calculated results, is the correct choice of the material model [11–13]. Many laboratory tests show that the behaviour of soils is highly non- linear [14–16] and therefore the use of a common linear elasto — perfectly plastic material models based on the Mohr-Coulomb yield criterion is not appropriate. Instead, the advanced material model like a newly developed Jardine plastic hardening-softening model (the JPHS model) should be used. 2. Theoretical background of the model The Jardine plastic hardening-softening model (the JPHS model) is based on the elasto — plastic theory as in geotechnics widely used Mohr-Coulomb model, but it contains some other features, which improve its capabilities and allow better simulations of the ground behaviour. Firstly, the Mohr-Coulomb fail- ure yield criterion is replaced by the Willam-Warnke failure yield criterion, which eliminates the singular tips from the Mohr-Coulomb surface and has a better agreement with the data from the experimental tests. Secondly, the model contains non-associated material description and the plastic flow is controlled by the Drucker-Prager potential function. Thirdly, the model has the ability to simulate the non-linear isotropic hardening and softening of the soils. Fourthly, it is capable to increase the stiffness modulus with an in- creasing depth or stress level. Fifthly, the volumetric hardening is controlled by a cap yield surface and finally, the model has the possibility to calculate the non-linear stress-strain dependence in a small strain range. 2.1. Yield failure criterion The Mohr-Coulomb failure criterion for soils is one of the oldest failure criteria. It is experimentally verified in triaxial compression and extension, but it is also very conservative for intermediate principal stress states between the triaxial compression and extension as can be seen in Figure 1. To eliminate this conservative behaviour and approximate the JPHS model more to the reality, the Willam-Warnke failure criterion [17] is implemented in the model. If rc is the distance from the hydrostatic axis to the failure surface at the compressive meridian and rt at the tension meridian, then at any intermediate position, the distance r is given by (see Figure 2): r = 2rc(r2c −r2t ) cos θ + rc(2rt −rc) √ D1 4(r2c −r2t ) cos2 θ + (rc − 2rt)2 , (1) where D1 = 4(r2c −r 2 t ) cos 2 θ + 5r2t − 4rtrc. (2) 58 http://dx.doi.org/10.14311/AP.2017.57.0058 http://ojs.cvut.cz/ojs/index.php/ap vol. 57 no. 1/2017 Numerical Modelling of Soil Behaviour Figure 1. The three-dimensional failure surface of Kaolin Clay in octahedral plane, lab. data from [18]. After several mathematical operations, the locus of yield surface in deviatoric plane can be written as: g(θ) = 2(1 −e2) cos θ + (2e− 1) √ D2 4(1 −e2) cos2 θ + (2e− 1)2 , (3) where D2 = 4(1 −e2) cos2 θ + 5e2 − 4e, (4) e = rc rt = 3 − sin φ′ 3 + sin φ′ . (5) And the whole Willam-Warnke failure criterion can be expressed as follows: F = A1p + q g(θ) + A2 = 0, (6) where A1 = − 6 − sin φ′ 3 − sin φ′ , A2 = − 6c− sin φ′ 3 − sin φ′ . (7) On the one hand, it is evident that the formula is relatively complex which can cause some difficulties during evaluation of differential and secondary order differential form, but on the other hand, this complex- ity allows that the locus of yield surface in deviatoric plane has the following features: • fits with the Mohr-Coulomb key points; i.e., the experimental extension and compression meridian points are on the curve; • it is differentiable at the compression and exten- sion meridian points; i.e., the singularities in the corner, which cause the numerical difficulties, are eliminated; • it is convex in the whole range of 0 < φ ≤ π/2 (0.5 < e ≤ 1). Figure 2. The geometrical interpretation of Willam- Warnke failure criteria. 2.2. Plastic flow The JPHS model contains the non-associated mate- rial description and the plastic flow is controlled by the Drucker-Prager yield criterion [19] that has the following form: F = √ J2 −αDPI1 −k = 0 (8) The coefficients αDP and kDP are not commonly used parameters in geotechnics and a proper calibration can lead to some difficulties in practice. Since the Drucker-Prager yield surface is a smooth version of the Mohr-Coulomb yield function, the Drucker-Prager constants can be expressed in terms of cohesion and friction angle. By relating hydrostatic stress p and deviatoric stress q to the invariants I1 and J2 the form of Drucker-Prager yield criterion can be rewritten as follows : F = q √ 3 − 3pαDP −k = 0, (9) → F = −3 √ 3pαDP + q − √ 3k = 0, (10) where the coefficient αDP and kDP can be expressed for the triaxial compression and extension: αCTC = 2 sin ψ √ 3(3 − sin ψ ,αCTE = 2 sin ψ √ 3(3 + sin ψ , (11) kCTC = 6c sin ψ √ 3(3 − sin ψ ,kCTE = 6c sin ψ √ 3(3 + sin ψ . (12) And the general form of the Drucker-Prager poten- tial function is then: F = A1p + q g(θ) + A2 = 0, (13) 59 Jan Vesely Acta Polytechnica Figure 3. The hardening-softening law — depen- dence of friction angle on equivalent deviatoric plastic deformation. where A1 = − 6 − sin ψ 3 − sin ψ ,A2 = − 6c− sin ψ 3 − sin ψ ,g(θ) = 1. (14) 2.3. Non-linear hardening/softening In the JPHS model, the isotropic hardening rule is used to describe the hardening and softening law. The change in yield criterion can be described as follows: dF = ∂F ∂σij dσij + h(d�p). (15) According to the strain hardening hypothesis [20], hardening process can be work hardening or strain hardening. To cover both cases, it is useful to define the hardening variable dκ and the parameter H, which express the dκ in terms of equivalent plastic strain d�p: H = dκ d�p . (16) The change in the yield function can then be rewritten as follows: dF = ∂F ∂σij dσij + ∂F ∂κ dκ = 0, (17) → dF = ∂F ∂σij dσij + ∂F ∂κ dκ d�p d�p dλ dλ = 0, (18) where the equivalent plastic strain can also be charac- terised as follows: d�p = √ 2 3 d�p : d�p = √ 2 3 dλ √ dQ. (19) Finally, substituting equations (16), (17) and (19) gives the expression of the hardening function hκ, which is directly related to the input parameter H: hκ = √ 2 3 H √ dQ ∂F ∂κ . (20) Figure 4. Cap yield surface in p-q plane. The JPHS model allows to assume the variation of friction angle φ′ with accumulated plastic strain as shown in Figure 3. There are three zones. In zone 1, φ′ is assumed to increase from the initial value (φ′crit) to the peak value (φ′peak) while in zone 2, φ ′ is reduced from the peak value to the residual value (φ′res) and in zone 3, φ′ remains constant and equal to the residual value (φ′res). In each of these zones, mathematical expressions can be assigned to the variation of φ′ with equivalent plastic strain d�p, and therefore the hardening/softening rules can be expressed in a piece- wise manner and the parameter H is defined in the JPHS model as follows: Hhardening = (φ′peak −φ ′)δ �hardening , (21) Hsoftening = (φ′ −φ′res)δ �softening . (22) Note that φ′peak, φ ′ res, δ, �hardening, �softening are in- put parameters, which are based on the result of the triaxial tests 2.4. Cap yield surface The shear yield surface defined in § 2.1 does not ex- plain the plastic volume strain that is measured in the isotropic compression. Therefore, the second type of yield surface is implemented in the JPHS model to close the elastic region for compressive stress path. The shape of cap yield surface in p-q plane (see Fig- ure 4) is chosen in the way to get the best agreement with the laboratory test data and is defined as follows: F = p + kcap − qc −q qc −pc = 0. (23) The initial pre-consolidation is calculated according to the following equation: pc,ini = OCR(pini + kcapqini), (24) qc,ini = 6 sin φ′peak 3 − sin φ′peak pc,ini + 6c sin φ′peak 3 − sin φ′peak . (25) During the volumetric strain hardening, the value of pre-consolidation pressure is then updated as follows: pc = pc,ini + K�v. (26) 60 vol. 57 no. 1/2017 Numerical Modelling of Soil Behaviour 2.5. Non-linear elasticity The strain range, in which the soils can be considered truly elastic, is very small and with increasing strain amplitude, soil stiffness decays non-linearly. To simu- late this soil behaviour, non-linear elastic stress-strain law, based on the Jardine function, is implemented in the model. Jardine et al. [14] proposed a periodic logarithmic function to express the non-linear relationship between the normalized secant Young’s modulus Eu and axial strain �a as follows: Eu cu = A + B cos ( α log (�a C )γ) , (27) where A, B, C, α, γ are the Jardine material model constants based on triaxial tests. The normalized tangent Young’s modulus Eut cor- responding to the secant Young’s modulus Eu can be derived by a differentiating and rearranging (27): Eut cu = A + B cos ( α log (�a C )γ) − Bαγ log ( �a C )γ−1 2.303 . (28) For the purpose of the numerical modelling, it is ap- propriate to replace the axial strain �a by deviatoric strain invariant �dev that is defined as: �dev = √ 2 3 ( (�1 − �2)2 + (�2 − �3)2 + (�3 − �1)2 ) . (29) By substituting the stress state of the undrained tri- axial test (�1 = �a and �2 = �3 = −11�a) into the equation (29), the dependence between axial strain �a and deviatoric strain invariant is �dev = √ 3�a and the normalized tangent Young’s modulus Eut can be rewritten as follows: Eut cu = A + B cos ( α log ( �dev√ 3C )γ) − Bαγ log ( �dev√ 3C )γ−1 2.303 . (30) Finally, to simulate the dependency of the modulus with the depth or stress level, the equation is modified as follows [21]: Eut = σm ( A + B cos(αIγ) − BαγIγ−1 2.303 sin(αIγ) ) , (31) where σm is mean stress pressure and I = log �dev√3C . Due to the Jardine model’s trigonometric nature, it is also mandatory to specify the exact strain range, in which the non-linear stress-strain law is to be applied. When exceeding the upper (�max) or lower (�min) limit of the strain range specified, stiffness is set as constant. A practical value for �min is the smallest strain for which the test data is available. For �max it is required to ensure compatibility with onset of plastic yield and check that high value will not gain to the negative tangent Young’s modulus (see Figure 5). Figure 5. Jardine function — strain range. 3. Implementation of the model The software PLAXIS allows users to implement a wide range of material models into the program. These models must be programmed in FORTRAN language, compiled as a Dynamic Link Library (DLL), and then added to the PLAXIS program directory. These mod- els simulate the soil behaviour in a single material point and the global behaviour is then governed by the Finite Element Method implementation in PLAXIS. The whole flow chart containing all steps needed for implementation of the model is presented in Figure 7. The first action that must be done is the declara- tion and initialization of the material properties and state variables and the specification of the undrained behaviour. The second action is to calculate the current value of the E-modulus. This is done by calling the following three subroutines: STRAIN_CALCULATION: This subroutine is calculating the deviatoric strain invariant �dev ac- cording to the formula (29). The input parameters are all components of the total strains and the out- put is the deviatoric strain invariant �dev. JARDINE: This subroutine is calculating the cur- rent E-modulus according to the formula (31) and rules mentioned in Section 2.5. The input parame- ters are the deviatoric strain invariant �dev and the JPHS model parameters. The output is the actual E-modulus UNLOADING: This subroutine is controlling the value of the E-modulus during unloading and reload- ing cycles according to the formula: Eunload = Ekunload, (32) where kunload is an unloading coefficient The fourth action is the calculation of the actual hard- ening / softening parameter by using the user defined subroutine HARDENING. The input of this subrou- tine is the JPHS model parameters. The output is a new value of the Hardening / Softening parameter H, determined according to the formulas (21) and (22). 61 Jan Vesely Acta Polytechnica Figure 6. Jardine model — secant stiffness defined as a logarithmic function of strain [14]. The last action is the definition of the constitutive stresses. After the calculation of the actual E-modulus, the new stiffness matrix is determined and the trial stresses are calculated as follows: predσn+1 = σn + D�n+1. (33) These stresses are one of the main inputs into the user-defined subroutine STRESS_INTEGRATION. This complex subroutine is the core element of the implementation and the backward return-mapping algorithm is implemented there. The first step in the subroutine is the check of the trial stresses. The JPHS model contains three different yield surfaces (Willam-Warnke yield surface Fs, Cap yield surface Fc and Tension cut-off yield surface Ft), and therefore the trial stresses can occur in 6 different regions: • elastic region (Fs < 0, Fc < 0 and Ft < 0); • Willam-Warnke region (Fs ≥ 0, Fc < 0 and Ft < 0); • cap region (Fs < 0, Fc ≥ 0 and Ft < 0); • tension cut-off region (Fs < 0, Fc < 0 and Ft ≥ 0); • Willam-Warnke and Tension cut-off region (Fs ≥ 0, Fc < 0 and Ft ≥ 0); • Willam-Warnke and Cap region (Fs ≥ 0, Fc ≥ 0 and Ft < 0). If the trial stress occurs in the elastic region (Fs < 0, Fc < 0 and Ft < 0), then the strain increment is elastic and the trial stress is the real stress for the increment: σn+1 = predσn+1. (34) If the trial stress occurs in other regions (Fs ≥ 0 or Fc ≥ 0 or Ft ≥ 0), then the trial stress is an inad- missible stress state and return-mapping algorithm is called. The general implementation of the return mapping algorithm is as follows [22]: (1.) Calculate the starting point: (a) calculate the derivatives (according to the re- gion) n+1nmn = ( ∂F ∂σmn ) n+1 , (35) n+1mpq = ( ∂Q ∂σpq ) n+1 ; (36) (b) update the plastic multiplier λ (yield function predF according to the region) dλ(0) = predF prednmnDijklpredmpq −hκ ; (37) (c) update stress and state variable (n+1)σ(0)mn = predσmn − dλ(0)Dmnpqpredmpq, (38) (n+1)κ(0) = predκ + dλ(0)hκ. (39) (2.) Calculate Backward Euler algorithm: DO WHILE n+1F (i) < tol (a) calculate the yield function n+1F (i) (according to the region) n+1F (i) = f ( n+1σ (i) ij , n+1κ(i) ) ; (40) (b) calculate the derivatives (according to the re- gion) n+1n(i)mn = ( ∂F ∂σmn ) n+1 , (41) n+1m (i) kl = ( ∂Q ∂σkl ) n+1 , (42) ∂mkl ∂σmn ∣∣∣∣(i) n+1 , ∂mkl ∂κ ∣∣∣∣(i) n+1 ; (43) 62 vol. 57 no. 1/2017 Numerical Modelling of Soil Behaviour Figure 7. JPHS model — flow chart. 63 Jan Vesely Acta Polytechnica (c) update the plastic multiplier λ (yield function predF according to the region) dλ(i+1) = n+1F (i) −n+1n(i)mnoldrij(n+1Tijmn)−1 D3 −n+1hκ , (44) where D3 = n+1n(i)mnDijkl n+1mkl(n+1Tijmn)−1 (45) n+1Tijmn = δimδnj + λ(i)Dijkl ∂mkl ∂σmn ∣∣∣∣(i) n+1 , (46) rij = σij − (predσij −λ(i)Dijkln+1m (i) kl ); (47) (d) update stress and state variable (n+1)σ(i+1)mn = predσmn − dλ(i)Dmnpqn+1mpq, (48) (n+1)κ(i) = predκ + dλ(i)hκ. (49) (e) i = i + 1 END DO 4. Validation of the model The following chapter is focused on the validation and verification of theJPHS model. The basic soil tests (the triaxial test and oedometer test) are simulated by using the numerical analyses and the results are compared with the real measured data from laboratory tests done on the Brno clay [15]. 4.1. The triaxial test The 2D analysis (PLAXIS v.2016) of the triaxial test is simulated by using the axisymmetric model. The input parameters for the JPHS model are summarized in Table 1 and Table 4. The coarseness of the mesh, geometry and dimensions of the model are presented in Figure 8. The left and the bottom boundaries are fixed in horizontal and vertical directions respectively. The isotropic / axial loading is represented by a distributed load and is applied on the top and right boundary. The modelling sequence is as follows: (1.) Initial conditions. (2.) Isotropic loading (The loads applied to the top and right boundary are activated, the undrained behaviour is ignored. After the isotropic loading, the prescribed displacements are set to zero). (3.) Axial compression to the failure (The load applied to the top is increased to a value which causes the failure, the undrained behaviour is active). Figure 9 and 10 presents the results of the triaxial tests for different values of initial stress conditions (275, 500 and 750 kPa). The results confirm that the JPHS model predicts a good match with the labora- tory data. Slightly over predicted is only the peak strength for 750 kPa, but this initial value corresponds Figure 8. The Triaxial Test — 2D model. Figure 9. The Triaxial Test Results — strain versus deviatoric stress, lab. data from [15]. to the depth of approximately 100 m that is not a com- mon area for the construction of geotechnical struc- tures. The comparison of water pressure development is shown in Figure 11 and again confirms a good match. The curve’s shape is almost identical with the labora- tory data and the slight underprediction of maximum excess pore pressure is only 15%. 4.2. Oedometer test The numerical simulation of the oedometer test is done by a 3D analysis (PLAXIS v.2013). The input paramters for the JPHS model are summarized in Table 1 and Table 4. The coarseness of the mesh, geometry and dimensions of the model are presented in Figure 12. The side boundaries are fixed in horizontal directions; the vertical fixities are used on the bottom 64 vol. 57 no. 1/2017 Numerical Modelling of Soil Behaviour Figure 10. The Triaxial Test Results — hydrostatic pressure versus deviatoric stress, lab. data from [15]. Figure 11. The Triaxial Test Results — strain versus water pressure, lab. data from [15]. side. The axial load σ1 is represented by a distributed load on the top plane. The modelling sequence is as follows: (1.) Initial conditions. (2.) Axial compression (the load applied on the top is gradually increased). (3.) Unloading (the load applied on the top is gradu- ally decreased). The results are shown in Figure 13 and confirm a good match between the JPHS model simulation and the laboratory data for both the loading and unloading stages. If the absolute values of axial strain are com- pared, then the JPHS model slightly underestimates the maximum value (approx. 10%). Figure 12. Oedometer Test — 3D model. Figure 13. Oedometer Test Results — strain versus stress, lab. data from [15]. 5. Practical use of the model — numerical modelling of tunnel excavation The 3D numerical analysis of the shallow tunnel in stiff clay is performed to verify the JPHS model and to compare the results calculated by this model with the real data from the geotechnical monitoring. For the calculation, the Kralovopolske tunnels (exploration adit) is chosen. The Kralovopolske tunnels are a part of the ring road of Brno town in the Czech Republic and are excavated in difficult geological conditions. The overburden varies between 6 m to 22 m, the tun- nels are excavated in Brno clay and there is an urban area on the surface. 5.1. Model discretization and boundary conditions The numerical analysis is carried out in the software PLAXIS 3D v.2013. Model represents 100 m wide, 65 Jan Vesely Acta Polytechnica Loess and clay loams Sand Deposits Stiff Clay Unit weight γ [kNm−3] 19.0 19.0 18.0 E-modulus E [MPa] 10 65 15 Poisson’s ratio ν [–] 0.35 0.35 0.40 Cohesion c’ [kPa] 10 5.0 5.0 Friction angle φ′ [°] 20 30 26 Dilatation angle ψ [°] 4.0 8.0 1.0 Table 1. MC model — geotechnical parameters. Cam-Clay swelling index Cam-Clay compression index Tangent of critical state line Initial void ratio Poisson’s ratio κ [–] λ [–] M [–] e [–] νur [–] Stiff Clay 0.2 0.09 1.55 0.5 0.2 Table 2. Modified Cam Clay model — geotechnical parameters. 48 m high and 90 m deep section of the soil mass, where the top boundary of the model represents the ground surface. The bottom and the side model boundaries are set at a distance required to reliably predict stress redistribution and ground deformation around the tunnel. The tunnel has a triangular shape with 4.8 m width and 4.2 m height and the crown of the tunnel is situated in a stiff clay approx. 22 m under the sur- face. The model geometry including the mesh and coordination system is presented in Figure 14. The model boundary conditions are set that the vertical displacement is fixed at the bottom boundary and the horizontal displacements are fixed at the side bound- aries. 5.2. Ground properties and primary stress conditions The area of Kralovopolske tunnels is formed by Miocene marine deposits of the Carpathian fore deep. The main geological strata, taken from the ground surface to bedrock, are loess and clay loams with a thickness between 3–10 m, followed by the layer of sand deposits (thickness approx. 10 m) and stiff clay (locally called Brno clay). The thickness of this clay is expected to be several hundreds of meters and most of the route of the tunnels is located there. The ground water level is connected with a sand layer and observed in the depth of 13–20 m. According to the study done by Svoboda et al. [15], material properties of the loams and sand deposits have an only small influence on the prediction of surface settlement, and therefore they are modelled only by using the Mohr-Coulomb model in this study. However, the stiff clay is modelled using the JPHS model and its prediction is then compared with the prediction calculated by the Mohr-Coulomb model, Modified Cam Clay model and Hardening Small Strain Figure 14. Model geometry. model. The model parameters are calibrated based on laboratory tests [15], which have been carried out during the site investigation. The parameters used for the numerical analysis are summarized in Table 1 to Table 4 and the calibration of the stiff clay for the JPHS model is shown in Figure 9, Figure 11, Figure 13 and Figure 15. The calibration of other models is mentioned in [23]. To set the initial stress conditions, it is necessary to determine the coefficient of earth pressure K0. The Brno clay can be characterized as over-consolidated clay with OCR = 6.5 [15]. By using the formula by Mayne and Kulhawy [24] the K0 coefficient is: K0 = (1 − sin φ)OCRsin φ = (1 − sin 26.5)6.5sin 26.5 = 1.25. (50) 5.3. Primary lining and modelling stages The Kralovopolske exploration adit is excavated by the New Austrian Tunnelling Method (NATM) with 66 vol. 57 no. 1/2017 Numerical Modelling of Soil Behaviour E-modulus Eref50 [MPa] Erefoed[MPa] E ref ur [MPa] Gref0 [MPa] γ0.7 [–] Stiff Clay 3.0 2.9 9.0 11.0 2.0 · 10−4 Power stress dependency Poisson’s ratio K0 value for NC Reference pressure m [–] νur [–] KNC0 [–] pref [kPa] Stiff Clay 1.0 0.2 0.54 100 Table 3. HSS model — geotechnical parameters. Stiff clay Jardine Stiffness Parameters A [–] 180 B [–] 165 C [–] 1.7 · 10−5 α [–] 0.85 γ [–] 1.05 �min [–] 1.701 · 10−5 �max [–] 0.02 Critical Friction angle φ′crit [°] 15.0 Peak Friction angle φ′peak [°] 28.0 Residual Friction angle φ′res [°] 18.0 Hardening constant �hardening [–] 0.07 Softening constant �softening [–] 0.04 Hard./soft. constant δ[–] 1.075 Unloading coefficient kunload [–] 2.7 Cap coefficient kcap [–] 0.45 Table 4. JPHS model — geotechnical parameters. Figure 15. JPHS model calibration- small strain stiffness, laboratory data from [15]. a full-face adit excavation. The numerical analysis corresponds to a construction process with the round length of 1.0 m that leads to 90 calculation phases. Each phase represents excavation progress of 1.0 m and installation of the tunnel lining in period of 8 hours. The tunnel lining consists of 100 mm of sprayed concrete, lattice girders (1 m span) and wire meshes. The lining is modelled using shell elements, which are capable of taking normal forces and bending moments. The lining is directly connected to the soil mass and is modelled as linearly elastic. To simulate the in- fluence of the lattice girders on the stiffness of the lining, the homogenization of a steel-concrete lining is done according to [25]. This procedure converses the cross-section of a lining consisting of two components with different E-modulus to a substitute-homogenized cross-section with only one modulus of elasticity. To simulate the increase in stiffness of the sprayed con- crete with the time, the Young’s modulus of shotcrete is updated in the construction stage process and is based on the Menschke equation [26]: Et = βEtE28, (51) where βEt =  0.0468t− 0.00211t 2 if t < 1 day, βEt = ( 0.9506 + 32.89 t−6 )−0.5 if t ≥ 1 day. (52) 5.4. Results of the numerical analysis Four 3D analyses are carried out to determine the effect of the material models on the deformation of the ground during the excavation of the adit in stiff clay. Three analyses are done using the standard material models included in software PLAXIS (Mohr Coulomb, Modified Cam-Clay and Hardening Small Strain model) while the fourth analysis is done using the JPHS model. In the following chapter, the results of all analyses are discussed and compared with the real data from geotechnical monitoring. Figure 16 presents the surface settlement after the full excavation of the adit. The results show that the settlement trough calculated by the MC model has a totally unrealistic shape. The vertical displacement above the tunnel is lower than at the certain distance from the axis. This phenomenon is caused by a high value of K0 and it is evident that the MC model is not suitable for a numerical modelling of the overcon- solidated soils. The other models predict better shape of settlement though and differs mainly in its width 67 Jan Vesely Acta Polytechnica Figure 16. The settlement trough after full profile excavation, monitoring data from [11]. Figure 17. Development of the settlement in the longitudinal direction, monitoring data from [15]. and depth. The closest results to the real monitored data are calculated by the JPHS model. If we look at the absolute values of the surface settlement, we can clearly see that the JPHS model almost predicts the same values (difference around 5%), while the HSS model underestimates the vertical displacement by more than 25%. Also, the width of the settlement trough calculated by the JPHS model is almost two times more accurate than the one predicted by the Modified Cam-Clay model. Another interesting results can be observed when the development of the settlement in the longitudinal direction is compared (see Figure 17). The development of the settlement in the longitu- dinal direction calculated by the MC and Modified Cam-Clay model has an unrealistic shape. The models are predicting that 80% of the deformation will occur before the adit face pass through the measured profile. Figure 18. Horizontal displacement after full profile excavation, monitoring data from [15]. Figure 19. Deformations of the adit lining, monitor- ing data from [15]. One of the reasons for this behaviour is that the MC and Modified Cam-Clay model does not contain the non-linear small-strain behaviour of soils. However, this feature is implemented in the HSS and the JPHS model and especially the JPHS model fits very well with the monitoring data. The similar results can also be observed for the horizontal displacement as presented in Figure 18. The JPHS model is closest to the real measured data. The difference is higher (approx. 50%), but one of the reasons can be the anisotropy of Brno Clay (mentined by [27]), which is not implemented in the model. The JPHS model is also the only model that correctly simulates the increase of the horizontal displacement near to the ground surface. The good match is also indicated when comparing 68 vol. 57 no. 1/2017 Numerical Modelling of Soil Behaviour the deformations of the adit lining. The differences between the material models are, in this case, not so significant but the JPHS model is still closest to the reality (see Figure 19). 6. Conclusions The presented results show that the choice of the ma- terial model has a significant influence on the correct modelling of the soil behaviour. The commonly used Mohr-Coulomb is not appropriate for the numerical analysis in the overconsolidated clay and it is neces- sary to use advanced material models. The presented JPHS model demonstrates that is able to quite pre- cisely simulate behaviour of the soils and with proper calibration of input parameters, it is possible to pre- dict the correct displacements of the ground during excavation. List of symbols αDP Drucker-Prager coefficient [–] δ Hard./soft. constant based on the triaxial test [–] �a Axial strain [–] �dev Deviatoric strain invariant [–] �hardening Hardening constant based on the triaxial test [–] �n Strain tensor [–] �softening Softening constant based on the triaxial test [–] �v Volumetric strain [–] d�p Plastic strain scalar [–] κ Hardening/softening variable [–] λ Plastic multiplier [–] φ′ Actual friction angle [rad] φ′peak Peak friction angle [rad] φ′res Residual friction angle [rad] ψ Dilatation angle [rad] σm Mean stress pressure [kPa] σn Stress tensor [kPa] θ Loge angle [rad] cu Undrained shear strength [kPa] kcap Cap coefficient [–] kunload Unload coefficient [–] p Hydrostatic stress [kPa] pc Pre-consolidation stress [kPa] q Deviatoric stress [kPa] rc Distance from the hydrostatic axis to the failure surface at the compressive meridian [m] rt Distance from the hydrostatic axis to the failure surface at the tension meridian [m] A,B,C,α,γ,�min,�max Jardine material model constants based on the triaxial test [–] D Stiffness matrix [–] Et Young’s modulus of sprayed concrete in time [kPa] E28 Young’s modulus of sprayed concrete after 28 days [kPa] Eunload Unloading Young’s modulus [kPa] Eu Secant Young’s modulus [kPa] Eut Tangent Young’s modulus [kPa] F Yield function [–] K Bulk modulus [kPa] OCR Over consolidation ratio [–] Q Potential function [–] Acknowledgements This work was supported by the Grant Agency of the Czech Technical University in Prague, grants SGS15/045/OHK1/1T/11 and SGS16/051/OHK1/1T/11. References [1] G.R.Dasari, C.G.Rawlings, M.D.Bolton. Numerical modelling of a natm tunnel construction in london clay. In Geotechnical Aspects of Underground Construction in Soft Ground, Rotterdam, pp. 491–496. 1996. [2] P.A.Vermeer, P. Bonnier, S.C.Moller. On a smart use of 3d-fem in tunnelling. In Proceedings of the 8th International Symposium on Numerical Models in Geomechanics - NUMOG VIII, Rome, Italy, pp. 361–366. 1996. doi:10.1201/9781439833797-c52. [3] M.Dolezalova. Approaches to numerical modelling of ground movements due to shallow tunnelling. In Proc. 2nd Int. Conference on Soil Structure Interaction in Urban Civil Engineering, ETH Zurich, pp. 365–376. 2002. [4] C. Ng, G. Lee. Three-dimensional ground settlements and stress-transfer mechanisms due to open-face tunnelling. Canadian Geotechnical Journal 42(4):1015–2029, 2005. doi:10.1139/t05-025. [5] L.Vydrova, J.Vesely. Optimization of the numerical modeling utilization for the design of undeground structures. In Geotechnical Engineering: New horizons. 2011. ISBN 978-1-60750-807-6. 2011. [6] S. Maras-Dragojevic. Analysis of ground settlement caused by tunnel construction. Gradevinar 64(7):573–581, 2012. [7] A. Lambrughi, L. Rodriguez, R. Castellanza. Development and validation of a 3d numerical model for tbm-epb mechanised excavations. Computers and Geotechnics 40:97–113, 2012. doi:http://dx.doi.org/10.1016/j.compgeo.2011.10.004. [8] T.Janda, M.Sejnoha, J.Sejnoha. Modeling of soil structure interaction during tunnel excavation: An engineering approach. Advances in Engineering Software 62-63:51–60, 2013. doi:http://dx.doi.org/10.1016/j.advengsoft.2013.04.011. [9] J.Bartak, J.Pruska, M.Hilar. Probability analysis of the effect of input parameters on the mrazovka tunnel deformations modelling. Tunel 11(3):27–33, 2002. [10] T.Svoboda, , M.Hilar. Probabilistic analyses of tunnel loads using variance reduction. Proceedings of the Institution of Civil Engineers - Geotechnical Engineering 168(4):348–357, 2015. doi:http://dx.doi.org/10.1680/geng.14.00062. [11] D.Masin. 3d modeling of an natm tunnel in high k0 clay using two different constitutive models. Journal of Geotechnical and Geoenvironmental Engineering 135:1326–1335, 2009. doi:10.1061/(ASCE)GT.1943-5606.0000017. 69 http://dx.doi.org/10.1201/9781439833797-c52 http://dx.doi.org/10.1139/t05-025 http://dx.doi.org/http://dx.doi.org/10.1016/j.compgeo.2011.10.004 http://dx.doi.org/http://dx.doi.org/10.1016/j.advengsoft.2013.04.011 http://dx.doi.org/http://dx.doi.org/10.1680/geng.14.00062 http://dx.doi.org/10.1061/(ASCE)GT.1943-5606.0000017 Jan Vesely Acta Polytechnica [12] D.Masin, I.Herle. Numerical analyses of a tunnel in london clay using different constitutive models. In Geotechnical aspects of underground construction in soft ground, At Amsterdam, pp. 795–801. 2015. doi:10.1201/NOE0415391245.ch81. [13] T. Benz. Small-Strain Stiffness of Soils and its Numerical Consequences. Ph.D. thesis, University of Stuttgart, 2010. [14] R. Jardine, D. Potts, A. Fourie, J. Burland. Studie of the influence of non-linear stress-strain characteristics in soil-structure interaction. Geotechnique 36(3):377–396, 1986. doi:10.1680/geot.1986.36.3.377. [15] T. Svoboda. Numericky model NRTM tunelu v tuhem jilu. Ph.D. thesis, Charles University in Prague, 2010. [16] A. Gasparre. Advanced Laboratory Characterisation of London Clay. Ph.D. thesis, Imperial College London, 2005. [17] K.J.Willam, E.P.Warnke. Constitutive model for triaxial behavior of concrete. In International Association for Bridges and Structural Engineering Proceedings. Bergamo, Italy, vol. 19. 1975. [18] A. Prashant. Three-Dimensional Mechanical Behavior of Kaolin Clay with Controlled Microfabric Using True Triaxial Testing. Ph.D. thesis, University of Tennessee, 2004. [19] D. Drucker, W. Prager. Soil mechanics and plastic analysis or limit design. Quarterly of applied mathematics 10(2):157–165, 1952. [20] M. Jirasek, Z. Bazant. Inelastic Analysis of Structures. London: J. Wiley & Sons, 2002. [21] B. Jones, A. Thomas, Y. Hsu, M. Hilar. Evaluation of innovative sprayed-concrete-lined tunnelling. Geotechnical Engineering 161(GE3):137–149, 2008. doi:10.1680/geng.2008.161.3.137. [22] B.Jeremic, S. Sture. Implicit integrations in elasto-plastic geotechnics. International Journal for Mechanics of Cohesive-Frictional Materials and Structures 2:165–183, 1997. doi:10.1002/(SICI)1099- 1484(199704)2:2<165::AID-CFM31>3.0.CO;2-3. [23] J. Vesely. The Use of Advanced Material Model for Numerical Modelling of Underground Structures in Clay. Ph.D. thesis, Czech Technical University in Prague, 2017. [24] P. Mayne, F. Kulhawy. Soil mechanics and plastic analysis or limit design. Journal of Geotechnical Engineering, ASCE 108(GT6):851–872, 1986. [25] J. Rott. Homogenisation and modification of composite steel-concrete lining, with the modulus of elasticity of sprayed concrete growing with time. Tunel 23(3):53–60, 2014. [26] A. Thomas. Sprayed concrete lined tunnels: An introduction. Taylor & Francis, 2009. [27] J.Rott, D. Masin, J.Bohac, et al. Evaluation of k0 in stiff clay by back-analysis of convergence measurements from unsupported cylindrical cavity. Acta Geotechnica 10:719–733, 2015. doi:10.1007/s11440-015-0395-7. 70 http://dx.doi.org/10.1201/NOE0415391245.ch81 http://dx.doi.org/10.1680/geot.1986.36.3.377 http://dx.doi.org/10.1680/geng.2008.161.3.137 http://dx.doi.org/10.1002/(SICI)1099-1484(199704)2:2<165::AID-CFM31>3.0.CO;2-3 http://dx.doi.org/10.1002/(SICI)1099-1484(199704)2:2<165::AID-CFM31>3.0.CO;2-3 http://dx.doi.org/10.1007/s11440-015-0395-7 Acta Polytechnica 57(1):58–70, 2017 1 Introduction 2 Theoretical background of the model 2.1 Yield failure criterion 2.2 Plastic flow 2.3 Non-linear hardening/softening 2.4 Cap yield surface 2.5 Non-linear elasticity 3 Implementation of the model 4 Validation of the model 4.1 The triaxial test 4.2 Oedometer test 5 Practical use of the model — numerical modelling of tunnel excavation 5.1 Model discretization and boundary conditions 5.2 Ground properties and primary stress conditions 5.3 Primary lining and modelling stages 5.4 Results of the numerical analysis 6 Conclusions List of symbols Acknowledgements References