# A Time-Dependent Double Hardening Model for Soft Rock

Authors: Mohammed Saiful Alam Siddiquee; Amin Hamdi
Journal: Advances in Civil Engineering
 (2019)

---

## Abstract

The time-dependent behavior of soft rock is very important in tunnel construction through soft rock. There is always a chance of creep failure due to sustained loading along the crest of the tunnel. In this research, a three-component elasto-visco-plastic framework is used to develop a time-dependent double hardening model to predict the behavior of soft rock both in compression and shear. Due to the limitation of time-dependent single yield hardening model in predicting the behavior of soft rock in compression direction (volumetric deformation) of loading, another time-dependent yield surface is added to the compression direction. The intersection of two-yield surface usually gives rise to singularity phenomenon, which is avoided by using Koiter’s generalization principle. A new formulation of the stress-integration of incremental stress-strain equations is proposed. A series of triaxial compression tests were carried out at different strain rates to calculate the model parameters. Then, the model is used to simulate the various behavioral features of soft rock.


---

## Body

## 1. Introduction

The loading rate effects due to material viscosity on the stress-strain behavior of sand (not due to delayed dissipation of excess pore water) are often very important in geotechnical engineering practice. A number of researchers (e.g., [1–10]) reported significant loading rate effects observed in laboratory stress-strain tests on sand under drained conditions; i.e., effects of strain rate and its change on the stress-strain relation, creep deformation, and stress relaxation during otherwise monotonic loading (ML) at a constant strain rate.Within the framework of the general nonlinear three-component model (Figure1), Benedetto et al. [5] and Tatsuoka et al. [7] proposed a set of stress-strain models to simulate the effects of material viscosity on the stress-strain behaviour of the geomaterial (i.e., clay, sand, gravel, and sedimentary soft rock). They showed that the viscous property of clean sand (i.e., uniform sand) is different from that of clay in that the viscous effect decays with an increase in the irreversible strain (temporary effect of strain rate and acceleration (TESRA)) and proposed a specific model to describe the above (i.e., the TESRA model explained below).Figure 1
General nonlinear three-component model (Benedetto [11]).In [5] and [7], it has been showed that at least three different functional forms of the viscous component, σv, are necessary to describe the different viscous properties of the geomaterial. Firstly, a stress-strain model called the “new isotach” was proposed to describe the loading rate effects of clay-like materials, for which, for primary ML along a fixed stress path, the current value of σv is a nonlinear function of instantaneous value of ε˙ir while it is always proportional to the instantaneous value of σf as follows:(1)σvεir,ε˙ir=σfεir⋅gvε˙ir,σ=σfεir⋅1+gvε˙ir,where gvε˙ir is the viscosity function, which is always zero or positive and given as follows for any strain (εir) or stress (σf) history (with or without cyclic loading):(2)gvε˙ir=α⋅1−exp1−ε˙ir/ε˙rir+1m≥0,where ε˙ir is the absolute value of ε˙ir and α, ε˙rir, and m are positive material constants. According to this model, as far as ML continues along a fixed stress path, the viscous stress component, σv, is a unique function of instantaneous values of εir and ε˙ir, independent of the previous loading history. The term “new” of the model name comes from that, with the original isotach model (Tatsuoka et al. [10]), the stress σ (therefore σv) is a function of instantaneous strain rate, ε˙=∂ε/∂t, not ε˙ir, while, with the new isotach model, σv is a function of ε˙ir. This difference results into significant differences in the model behaviour, in particular during stress relaxation with ε˙=0 and immediately after a step change in ε˙ while doing an ML at a constant strain rate ε˙ (Tatsuoka et al. [10]).The above family of three-component material models can be solved by using three techniques, e.g., (1) Perzyna [11–14] method, (2) Duvaut–Lions [15–17] method, and (3) consistency method (Wang et al. [18]). The first two methods are based on the idea that stress state can stay outside the yield surface and directly define the plastic relaxation equations in the stress space. But in the consistency method, viscoplastic regularization is attained by introducing a rate-dependent yield surface and the consistency condition is fulfilled. In this proposal, the consistency method is used to integrate the incremental visco-elasto-plastic equations.The solution of two or more yield surface simultaneously is adopted from De Borst [19], where it has been showed that singularity in multiple yield surfaces can be avoided using the generalization of Koiter’s principle of energy dissipation.The TESRA model will be implemented into an existing nonlinear elastoplastic FE code ([20]), which will be developed as an user model written in FORTRAN 90 that was highly optimized for very fast and accurate computation of highly nonlinear equations emerging from material nonlinearity. Due to the inherent features of the commercially available nonlinear finite element program “Plaxis,” the model will be added to the program as DLL file.At the first level of integration, incremental elasto-visco-plastic equations will be solved inside the “return mapping algorithm” (Ortiz and Simo [20]), where the stress is returned to the growing yield surface while fulfilling the consistency condition (abiding by the flow rule). At each step of return mapping iteration, when necessary, the second level of integration is made, where the stress is returned to the inviscid yield surface with an incremental integration based on the TESRA model while calculating the viscous stress.


## 2. Element Test Results

In order to determine the model parameters, a series of triaxial tests Miyashita et al. [21] were conducted on the softrock material. A picture of a triaxial cell and a schematic figure of a specimen and deformation transducers used in this study are shown in Figure 2.Figure 2
(a) Picture of triaxial apparatus; (b) schematic figure of specimen and deformation transducers.
(a)
(b)


### 2.1. Test Materials

Undisturbed samples by block sampling taken from tunnel construction site were trimmed to rectangular specimen with the dimensions of 58 mm in width (W), 76 mm in depth (D), and 160 mm in height (H). The location of the samples was the connection tunnel between two express ways at a depth of about 45 m.


### 2.2. Test Preparation

Dental gypsum was put to the top and bottom end surfaces of specimen to reduce the effect of bedding error. The specimen on the UC test was covered by a membrane to retain water content, and saturation process was not applied. Specimens on CD tests were saturated by double vacuuming method and applied a back pressure of 200 kPa. The range of strain rate measured by LVDT was 0.0015%/min.–1.5%/min. Unconfined compression test and a series of drained triaxial compression tests were conducted. Small strain cyclic loading was applied at the end of isotropic consolidation and during shearing (5 cycles on each step). Some stages of creep loading were applied to CD_400, 1600, and 2700. Creep failure occurred at the deviator stress of 2545 kPa in the test CD_400.


### 2.3. Apparatus

Combination of linear roller bearings put perpendicular to each other below the pedestal made it possible not to interfere with the deformation or displacement of the specimen. Axial load was measured by using an inner load cell to avoid the influence of loading shaft friction. Global axial strain was measured by an external displacement transducer (herein denoted as EDT). The effect of bedding error was reduced by introducing thin layers of dental gypsum at the top and bottom ends of the specimen. In addition, local measurements of axial and horizontal strains were also conducted using vertical and horizontal local deformation transducers (herein denoted as LDTs; refer [8] for the details). Two vertical and two horizontal LDTs were placed on opposite sides of the specimen as shown in Figure 2(b). To evaluate global volumetric strain during testing by measuring the volume of drained water from the specimen, low capacity differential pressure transducer (herein denoted as LCDPT) was employed. The summation of vertical and horizontal strains obtained by LDTs was also evaluated as volumetric strain without the effect of bedding error and compliance associated with drainage system from the specimen to the LCDPT.During axial loading, the axial strain rate was changed stepwise many times by a factor of 10 in all the tests. The maximum value of the strain rate was 1000 times larger than the minimum value which was equal to 0.0015%/min. In the test CD_400, two stages of sustained loading were applied during monotonic loading keeping other test conditions unchanged.


### 2.4. Test Program (Loading Scheme/Stress Paths)

The loading scheme is given in Table1.Table 1
ρ
t0 (g/cm3)
B value (after IC)
σ
c′ (kPa)
Creep loading
Small strain cyclic loading
End of IC
During shearing
UC
1.884
—
—
N
Y
N
CD_100
1.896
0.95
100
N
Y
N
CD_400
1.894
0.90
400
Y
Y
Y
CD_800
1.896
0.86
800
N
Y
Y
CD_1600
1.875
0.61
1600
Y
Y
Y
CD_2700
1.897
0.79
2700
Y
Y
Y


## 2.1. Test Materials

Undisturbed samples by block sampling taken from tunnel construction site were trimmed to rectangular specimen with the dimensions of 58 mm in width (W), 76 mm in depth (D), and 160 mm in height (H). The location of the samples was the connection tunnel between two express ways at a depth of about 45 m.


## 2.2. Test Preparation

Dental gypsum was put to the top and bottom end surfaces of specimen to reduce the effect of bedding error. The specimen on the UC test was covered by a membrane to retain water content, and saturation process was not applied. Specimens on CD tests were saturated by double vacuuming method and applied a back pressure of 200 kPa. The range of strain rate measured by LVDT was 0.0015%/min.–1.5%/min. Unconfined compression test and a series of drained triaxial compression tests were conducted. Small strain cyclic loading was applied at the end of isotropic consolidation and during shearing (5 cycles on each step). Some stages of creep loading were applied to CD_400, 1600, and 2700. Creep failure occurred at the deviator stress of 2545 kPa in the test CD_400.


## 2.3. Apparatus

Combination of linear roller bearings put perpendicular to each other below the pedestal made it possible not to interfere with the deformation or displacement of the specimen. Axial load was measured by using an inner load cell to avoid the influence of loading shaft friction. Global axial strain was measured by an external displacement transducer (herein denoted as EDT). The effect of bedding error was reduced by introducing thin layers of dental gypsum at the top and bottom ends of the specimen. In addition, local measurements of axial and horizontal strains were also conducted using vertical and horizontal local deformation transducers (herein denoted as LDTs; refer [8] for the details). Two vertical and two horizontal LDTs were placed on opposite sides of the specimen as shown in Figure 2(b). To evaluate global volumetric strain during testing by measuring the volume of drained water from the specimen, low capacity differential pressure transducer (herein denoted as LCDPT) was employed. The summation of vertical and horizontal strains obtained by LDTs was also evaluated as volumetric strain without the effect of bedding error and compliance associated with drainage system from the specimen to the LCDPT.During axial loading, the axial strain rate was changed stepwise many times by a factor of 10 in all the tests. The maximum value of the strain rate was 1000 times larger than the minimum value which was equal to 0.0015%/min. In the test CD_400, two stages of sustained loading were applied during monotonic loading keeping other test conditions unchanged.


## 2.4. Test Program (Loading Scheme/Stress Paths)

The loading scheme is given in Table1.Table 1
ρ
t0 (g/cm3)
B value (after IC)
σ
c′ (kPa)
Creep loading
Small strain cyclic loading
End of IC
During shearing
UC
1.884
—
—
N
Y
N
CD_100
1.896
0.95
100
N
Y
N
CD_400
1.894
0.90
400
Y
Y
Y
CD_800
1.896
0.86
800
N
Y
Y
CD_1600
1.875
0.61
1600
Y
Y
Y
CD_2700
1.897
0.79
2700
Y
Y
Y


## 3. Viscoplastic Framework

### 3.1. Shear Yield Behaviour

This part of the modeling is achieved by frictional hardening with the apparent cohesion of the material. A strain-rate-dependent phenomenological model has been developed to describe the behavior of soft rock. As soft rock involves a significant amount of cohesion, the model is developed within the framework of frictional hardening but including the cohesion.From Figure3, the definition of normal stresses can be redefined as follows:(3)σij′=σij+Ccotφ.Figure 3
Schematic diagram of failure criterion of a cohesive frictional hardening model.
(a)
(b)Using the extendedx-axis as equivalent normal stress (as in Figure 3), the yield function and plastic potential of the elasto-visco-plastic framework can be easily established in usual way.This is essentially an elastoplastic isotropic strain hardened analysis with nonassociated flow rule:(4)f1=αθI1+J2=0,where(4a)αθ=sinϕmob3cosθ−3cosθsinϕmob,where I1 is the 1st stress invariant, J2 is the 2nd invariant of deviatoric stress, K is the cohesion terms (used only in Figure 3(b)), α(θ) is the Lode angle function, and sinϕmob is the mobilized angle of internal friction for the modified normal stress definition. In this way, sinϕmob is defined as follows:(4b)sinϕmob=σ1′−σ3′σ1′+σ3′=σ1+Ccotφ−σ3+Ccotφσ1+Ccotφ+σ3+Ccotφ.The original Rowe’s stress-dilatancy relation was modified to define the plastic potential as stated as follows:(5)g1=α′θI1+J2=0,where(6)α′θ=sinψ3cosθ−3cosθsinψ,sinψ=900σmsinϕmob−sinϕr1−sinϕmobsinϕr,And, sinψ has a lower limit as follows: if sinψ<0.105, sinψ=0.105.This dilatancy cutoff is very important for large strain computation in the softening regions, whereϕmob is the mobilized angle of friction, sinϕr is the residual angle of internal friction, ψ is the mobilized angle of dilatancy, and σm is the mean pressure in kPa.The evolution of the yield function is modeled as described in:(7)κ=∫dεP,where dεp = [2 (deP11)2 + 2 (deP22)2 + 2 (deP33)2 + (dγP11)2 +(dγP23)2 + (dγP31)2]1/2dePij ∼ dγPij are the deviatoric components of plastic strain increments and κ = plasticity parameter (internal variable).The evolution equation is defined by the following sets of equations:(8)sinϕmob=σ1′−σ3′σ1′+σ3′=R′κ−1R′κ+1.Ifκ<κf then,(9)R′κ=Rmax−12κ×κfκ+κfbcont+1.0.Ifκ>κf then,(10)R′κ=Rr+Rmax−Rr×Expκ−κf∗κ−κfccont∗ccont.


### 3.2. Compression Yield Surface

The compressional yield surface is a yield surface controlling the yielding in the direction of mean pressure (volumetric strain). It is given by(11)f2=I1+I0=0,where I1 is the 1st stress invariant of deviatoric stresses and I0 is the evolution parameter in the compression direction given by(12)I0=δ⋅εvp,where δ is the material constant signifying the rate of movement of I0 in the compression direction and εvp is the plastic volumetric strain arising out of model response.


### 3.3. Viscous Model

As stated in several papers by Prof. Tatsuoka [7, 10, 20, 22], the following equations are used for the implementation of viscous behavior (the New Isotach model) within the framework of “frictional hardening model with cohesion.” Following the fundamental formulation of the new isotach model (equations (1) and (2)), the fictional hardening model with cohesion can be formulated as follows.Equations (9) and (10) are recast into the following equation in order to follow the modified framework of the “frictional hardening model with cohesion:”(13)Rv′=R′κ⋅1+gvε˙ir,where R′κ is described in equation (10) and gvε˙ir is explained in equations (1) and (2).The same viscous formulation is used here for the second yield surface in compression direction. Equations (11) and (12) are recast into equation (14). So, the viscous formulation for the compression direction can be given as(14)δv=δ⋅1+gvε˙ir.


### 3.4. Elastic Behavior

The directional stress dependent cross-anisotropic model is used in this analysis. Here in this analysis, the following equations of elastic modulus and Poisson’s ratio are used to model the elastic behavior:(15)Ev=E1σvm,Eh=E1σhm,where E1 = “Ev or Eh at σv or σh = 1.0 kN/m2.”E
v and Eh are isotropic (Ev = Eh = Ep), when σv = σh, as supported by the data of Tatsuoka et al., [7]. Poisson’s ratio is determined according to equation (16). Using Young’s modulus and Poisson’s ratio, a cross-anisotropic Hook’s law is used to determine the elastic response of the model:(16)νhv=ν0σhσvm/2,νvh=ν0σvσhm/2.


## 3.1. Shear Yield Behaviour

This part of the modeling is achieved by frictional hardening with the apparent cohesion of the material. A strain-rate-dependent phenomenological model has been developed to describe the behavior of soft rock. As soft rock involves a significant amount of cohesion, the model is developed within the framework of frictional hardening but including the cohesion.From Figure3, the definition of normal stresses can be redefined as follows:(3)σij′=σij+Ccotφ.Figure 3
Schematic diagram of failure criterion of a cohesive frictional hardening model.
(a)
(b)Using the extendedx-axis as equivalent normal stress (as in Figure 3), the yield function and plastic potential of the elasto-visco-plastic framework can be easily established in usual way.This is essentially an elastoplastic isotropic strain hardened analysis with nonassociated flow rule:(4)f1=αθI1+J2=0,where(4a)αθ=sinϕmob3cosθ−3cosθsinϕmob,where I1 is the 1st stress invariant, J2 is the 2nd invariant of deviatoric stress, K is the cohesion terms (used only in Figure 3(b)), α(θ) is the Lode angle function, and sinϕmob is the mobilized angle of internal friction for the modified normal stress definition. In this way, sinϕmob is defined as follows:(4b)sinϕmob=σ1′−σ3′σ1′+σ3′=σ1+Ccotφ−σ3+Ccotφσ1+Ccotφ+σ3+Ccotφ.The original Rowe’s stress-dilatancy relation was modified to define the plastic potential as stated as follows:(5)g1=α′θI1+J2=0,where(6)α′θ=sinψ3cosθ−3cosθsinψ,sinψ=900σmsinϕmob−sinϕr1−sinϕmobsinϕr,And, sinψ has a lower limit as follows: if sinψ<0.105, sinψ=0.105.This dilatancy cutoff is very important for large strain computation in the softening regions, whereϕmob is the mobilized angle of friction, sinϕr is the residual angle of internal friction, ψ is the mobilized angle of dilatancy, and σm is the mean pressure in kPa.The evolution of the yield function is modeled as described in:(7)κ=∫dεP,where dεp = [2 (deP11)2 + 2 (deP22)2 + 2 (deP33)2 + (dγP11)2 +(dγP23)2 + (dγP31)2]1/2dePij ∼ dγPij are the deviatoric components of plastic strain increments and κ = plasticity parameter (internal variable).The evolution equation is defined by the following sets of equations:(8)sinϕmob=σ1′−σ3′σ1′+σ3′=R′κ−1R′κ+1.Ifκ<κf then,(9)R′κ=Rmax−12κ×κfκ+κfbcont+1.0.Ifκ>κf then,(10)R′κ=Rr+Rmax−Rr×Expκ−κf∗κ−κfccont∗ccont.


## 3.2. Compression Yield Surface

The compressional yield surface is a yield surface controlling the yielding in the direction of mean pressure (volumetric strain). It is given by(11)f2=I1+I0=0,where I1 is the 1st stress invariant of deviatoric stresses and I0 is the evolution parameter in the compression direction given by(12)I0=δ⋅εvp,where δ is the material constant signifying the rate of movement of I0 in the compression direction and εvp is the plastic volumetric strain arising out of model response.


## 3.3. Viscous Model

As stated in several papers by Prof. Tatsuoka [7, 10, 20, 22], the following equations are used for the implementation of viscous behavior (the New Isotach model) within the framework of “frictional hardening model with cohesion.” Following the fundamental formulation of the new isotach model (equations (1) and (2)), the fictional hardening model with cohesion can be formulated as follows.Equations (9) and (10) are recast into the following equation in order to follow the modified framework of the “frictional hardening model with cohesion:”(13)Rv′=R′κ⋅1+gvε˙ir,where R′κ is described in equation (10) and gvε˙ir is explained in equations (1) and (2).The same viscous formulation is used here for the second yield surface in compression direction. Equations (11) and (12) are recast into equation (14). So, the viscous formulation for the compression direction can be given as(14)δv=δ⋅1+gvε˙ir.


## 3.4. Elastic Behavior

The directional stress dependent cross-anisotropic model is used in this analysis. Here in this analysis, the following equations of elastic modulus and Poisson’s ratio are used to model the elastic behavior:(15)Ev=E1σvm,Eh=E1σhm,where E1 = “Ev or Eh at σv or σh = 1.0 kN/m2.”E
v and Eh are isotropic (Ev = Eh = Ep), when σv = σh, as supported by the data of Tatsuoka et al., [7]. Poisson’s ratio is determined according to equation (16). Using Young’s modulus and Poisson’s ratio, a cross-anisotropic Hook’s law is used to determine the elastic response of the model:(16)νhv=ν0σhσvm/2,νvh=ν0σvσhm/2.


## 4. Stress Integration

The stress integration is carried out using the return mapping algorithm (Simo [16]) along with the consistent viscoplastic flow theory (Wang et al., [18]). A consistent double hardening formulation is developed depending on different cases of loading in various directions. There are obviously two cases. The first case is the separate yielding of shear yield surface and compression yield surface independent of each other. The second case is the simultaneous yielding of both surfaces. Simultaneous yielding of both yield surfaces happens at the corner of the intersection of the two-yield surfaces as shown in Figure 4. In the following subsections, the stress integration steps are detailed according to the return mapping algorithm. It is noteworthy that here stress integration takes care of the rate of plastic loading as well [23, 24].Figure 4
Schematic diagram of two interacting yield surfaces.


### 4.1. Case I: Separate Yielding

Case I describes the stress integration about the yielding in shear direction due to the increase in the slope of shear yield surface and the yielding in compression direction. The evolution of the plastic flow in both directions is determined by the discretized consistency condition of the of the rate-dependent yield surface, which can be expanded using Tailor’s series as follows:(17)f˙1σn+1,κn+1,κ˙n+1=0.0,f1+∂f1∂σδσ+∂f1∂κδk+∂f1∂κ˙δk˙=0.0,(18)f˙2σn+1,μn+1,μ˙n+1=0.0,f2+∂f2∂σδσ+∂f2∂μδμ+∂f2∂μ˙δμ˙=0.0.From the additive decomposition of strain in a framework of small strain theory, the following equation can be written:(19)σ=De:ε−De:εp,σ=De:ε−De:λ1∂g1∂σ.Now, taking differential variation of the following equation (19), we can get the value of δσ:(20)δσ=De:δε−De:δλ1∂g1∂σ+λ1∂2g1∂σ∂σδσ+λ1∂2g1∂σ∂κδκ+λ1∂2g1∂σ∂κ˙δκ˙.During the incremental elastic loading step, the elastoplastic integration of the incremental elastoplastic equation through return mapping total strain (ε) does not change. So the following equation (21) can be written:(21)δε=0.0.Putting equation (21) into equation (20), we get equations (22) and (24):(22)δσ+λ1De:∂2g1∂σ∂σδσ=−De:δλ1∂g1∂σ+λ1∂2g1∂σ∂κδκ+λ1∂2g1∂σ∂κ˙δκ˙,δσI+λ1De:∂2g1∂σ∂σ=−De:δλ1∂g1∂σ+λ1∂2g1∂σ∂κδκ+λ1∂2g1∂σ∂κ˙δκ˙,δσ=−De:I+λ1De:∂2g1∂σ∂σ−1δλ1∂g1∂σ+λ1∂2g1∂σ∂κδκ+λ1∂2g1∂σ∂κ˙δκ˙,δσ=−De−1+λ1∂2g1∂σ∂σ−1δλ1∂g1∂σ+λ1∂2g1∂σ∂κδκ+λ1∂2g1∂σ∂κ˙δκ˙.Similarly another equation (23) can be written also for the plastic potential in the compression direction, g2:(23)δσ=−De−1+λ1∂2g2∂σ∂σ−1δλ1∂g2∂σ+λ1∂2g2∂σ∂μδμ+λ1∂2g2∂σ∂μ˙δμ˙.Now, substituting the value ofδσ from equations (22) and (23) into lower part of equations (17) and (18), respectively, we get(24)f1+∂f1∂σδσ+∂f1∂κδk+∂f1∂κ˙δk˙=0.0,f1−∂f1∂σDe−1+λ1∂2g1∂σ∂σ−1δλ1∂g1∂σ+λ1∂2g1∂σ∂κδκ+λ1∂2g1∂σ∂κ˙δκ˙+∂f1∂κδk+∂f1∂κ˙δk˙=0.0,f2+∂f2∂σδσ+∂f2∂κδμ+∂f2∂κ˙δμ˙=0.0,f2−∂f2∂σDe−1+λ2∂2g2∂σ∂σ−1δλ2∂g2∂σ+λ2∂2g2∂σ∂μδμ+λ2∂2g2∂σ∂μ˙δμ˙+∂f2∂μδμ+∂f2∂μ˙δμ˙=0.0.Then, defining hardening modulusH as H1=−∂f1/∂κ;H2=−∂f2/∂μ and rate of hardening modulus Hr as H1r=−∂f1/∂κ˙;H2r=−∂f2/∂μ˙ and it is well-known that for plane strain case that δλ=δκandδλ˙=δκ˙, equation (24) can be rewritten as follows:(25)f1−∂f1∂σDe−1+λ1∂2g1∂σ∂σ−1δλ1∂g1∂σ+λ1∂2g1∂σ∂κδκ+λ1∂2g1∂σ∂κ˙δκ˙+∂f1∂κδk+∂f1∂κ˙δk˙=0,∂f1∂σDe−1+λ1∂2g1∂σ∂σ−1δλ1∂g1∂σ+λ1∂2g1∂σ∂κδλ+λ1∂2g1∂σ∂κ˙δκ˙+H1δλ+H1rδk˙=f1,δλ1=f1H1+H1r/Δt+∂f1/∂σDe−1+λ1∂2g1/∂σ∂σ−1∂g1/∂σ+λ1∂2g1/∂σ∂κ+λ1/Δt∂2g1/∂σ∂κ˙.Similarly, for the compression zone, the following equation can be written as follows:(26)δλ2=f2H2+H2r/Δt+∂f1/∂σDe−1+λ2∂2g1/∂σ∂σ−1∂g2/∂σ+λ2∂2g2/∂σ∂κ+λ2/Δt∂2g2/∂σ∂κ˙.Equation (11) can be rewritten in the following conventional way:(27)δλ1=f1H1+H1r/Δt+∂f1/∂σDe−1+λ1∂2g1/∂σ∂σ−1∂g1/∂σ+λ1∂2g1/∂σ∂κ+λ1/Δt∂2g1/∂σ∂κ˙,δλ1=f1H1+H1r/Δt+n1tD¯1em¯1,where(28)n1t=∂f1∂σ,D¯1e=D1e−1+λ1∂2g1∂σ∂σ−1,m¯1=∂g1∂σ+λ1∂2g1∂σ∂κ+λ1Δt∂2g1∂σ∂κ˙.Similarly, equation (26) for the compression zone can be also rewritten as(29)δλ2=f2H2+H2r/Δt+∂f2/∂σDe−1+λ2∂2g2/∂σ∂σ−1∂g2/∂σ+λ2∂2g2/∂σ∂κ+λ2/Δt∂2g2/∂σ∂κ˙,δλ2=f2H2+H2r/Δt+n2tD¯2em¯2,where(30)n2t=∂f2∂σ,D¯2e=D2e−1+λ2∂2g2∂σ∂σ−1,m¯2=∂g2∂σ+λ2∂2g2∂σ∂κ+λ2Δt∂2g2∂σ∂κ˙.Once, the increment of plasticity parameter,δλ is known for shear and compression direction, corrected stress can be calculated, and plasticity parameter can be updated as follows:(31)σn+1=σtrial−δλ1De∂g1∂σ,κn+1=κn+δλ1,σn+1=σtrial−δλ2De∂g1∂σ,μn+1=μn+δλ2.


### 4.2. Case II: Simultaneous Yielding

Case II describes the stress integration when both the yield surfaces are yielding at the same time. During the general loading, it is possible that both the yield functions are acting simultaneously. It represents a stress state lying inside the corner, bounded by the two-yield functions (Figure4). Following the same arguments as case I, the following equations can be written:(32)f˙1σn+1,κn+1,κ˙n+1=0.0,f1+∂f1∂σδσ+∂f1∂κδk+∂f1∂κ˙δk˙=0.0,(33)f˙2σn+1,μn+1,μ˙n+1=0.0,f2+∂f2∂σδσ+∂f2∂μδμ+∂f2∂μ˙δμ˙=0.0.From the additive decomposition of strain in a framework of small strain theory and Koiter's generalization of plastic strain generation, the following equation can be written:(34)σ=De:ε−De:εp,σ=De:ε−De:λ1∂g1∂σ−De:λ2∂g2∂σ.Now, taking differential variation of the following equation (34), the value of δσ can be deduced as follows:(35)δσ=De:δε−De:δλ1∂g1∂σ+λ1∂2g1∂σ∂σδσ+λ1∂2g1∂σ∂κδκ+λ1∂2g1∂σ∂κ˙δκ˙−De:δλ2∂g2∂σ+λ2∂2g2∂σ∂σδσ+λ2∂2g2∂σ∂μδμ+λ2∂2g1∂σ∂μ˙δμ˙.Usingδε=0.0 into equation (35), the following equation (36) can be written:(36)δσ+λ1De:∂2g1∂σ∂σδσ+λ2De:∂2g2∂σ∂σδσ=−De:δλ1∂g1∂σ+λ1∂2g1∂σ∂κδκ+λ1∂2g1∂σ∂κ˙δκ˙−De:δλ2∂g2∂σ+λ2∂2g2∂σ∂μδμ+λ1∂2g2∂σ∂μ˙δμ˙,δσ=−Rδλ1∂g1∂σ+λ1∂2g1∂σ∂κδκ+λ1∂2g1∂σ∂κ˙δκ˙+δλ2∂g2∂σ+λ2∂2g2∂σ∂μδμ+λ1∂2g2∂σ∂μ˙δμ˙,where R=De−1+λ1∂2g1/∂σ∂σ++λ2∂2g2/∂σ∂σ−1.Now substituting the value ofδσ from equation (36) into equation (32), the following equation can be derived:(37)f1−∂f1∂σRδλ1∂g1∂σ+λ1∂2g1∂σ∂κδκ+λ1∂2g1∂σ∂κ˙δκ˙+δλ2∂g2∂σ+λ2∂2g2∂σ∂μδμ+λ2∂2g2∂σ∂μ˙δμ˙+∂f1∂κδk+∂f1∂κ˙δk˙=0.0,∂f1∂κ+∂f1Δt∂κ˙+∂f1∂σR∂g1∂σ+∂f1∂σRλ1∂2g1∂σ∂κ+∂f1∂σRλ1Δt∂2g1∂σ∂κ˙δλ1+∂f1∂σR∂g2∂σ+∂f1∂σRλ2∂2g2∂σ∂μ+∂f1∂σRλ2Δt∂2g2∂σ∂μ˙δλ2=f1.Again substituting the value ofδσ from equation (36) into equation (33), the following equation can be derived:(38)f2−∂f2∂σRδλ1∂g1∂σ+λ1∂2g1∂σ∂κδκ+λ1∂2g1∂σ∂κ˙δκ˙+δλ2∂g2∂σ+λ2∂2g2∂σ∂μδμ+λ2∂2g2∂σ∂μ˙δμ˙+∂f2∂μδμ+∂f2∂μ˙δμ˙=0.0,∂f2∂μ+∂f2Δt∂μ˙+∂f2∂σR∂g2∂σ+∂f2∂σRλ2∂2g2∂σ∂μ+∂f2∂σRλ2Δt∂2g2∂σ∂μ˙δλ2+∂f2∂σR∂g1∂σ+∂f2∂σRλ1∂2g1∂σ∂κ+∂f2∂σRλ1Δt∂2g1∂σ∂κ˙δλ1=f2.Abbreviating and rearranging equations (37) and (38), the following equation (39) can be written:(39)β1δλ1+β2δλ2=f1,β3δλ1+β4δλ2=f2,where(40)β1=∂f1∂κ+∂f1Δt∂κ˙+∂f1∂σR∂g1∂σ+∂f1∂σRλ1∂2g1∂σ∂κ+∂f1∂σRλ1Δt∂2g1∂σ∂κ˙,β2=∂f1∂σR∂g2∂σ+∂f1∂σRλ2∂2g2∂σ∂μ+∂f1∂σRλ2Δt∂2g2∂σ∂μ˙,β3=∂f2∂σR∂g1∂σ+∂f2∂σRλ1∂2g1∂σ∂κ+∂f2∂σRλ1Δt∂2g1∂σ∂κ˙,β4=∂f2∂μ+∂f2Δt∂μ˙+∂f2∂σR∂g2∂σ+∂f2∂σRλ2∂2g2∂σ∂μ+∂f2∂σRλ2Δt∂2g2∂σ∂μ˙.Now the systems of simultaneous equations (39) can be solved, and δλ1 and δλ2 are determined as follows:(41)δλ1=β4f1−β2f2β1β4−β2β3,δλ2=β1f2−β3f1β1β4−β2β3,Once, the increment of plasticity parameters,δλ1 and δλ2, are known, the corrected stress can be calculated and plasticity parameter can be updated as follows:(42)σn+1=σtrial−δλ1De∂g1∂σ−δλ2De∂g2∂σ,κn+1=κn+δλ1,μn+1=μn+δλ2.


## 4.1. Case I: Separate Yielding

Case I describes the stress integration about the yielding in shear direction due to the increase in the slope of shear yield surface and the yielding in compression direction. The evolution of the plastic flow in both directions is determined by the discretized consistency condition of the of the rate-dependent yield surface, which can be expanded using Tailor’s series as follows:(17)f˙1σn+1,κn+1,κ˙n+1=0.0,f1+∂f1∂σδσ+∂f1∂κδk+∂f1∂κ˙δk˙=0.0,(18)f˙2σn+1,μn+1,μ˙n+1=0.0,f2+∂f2∂σδσ+∂f2∂μδμ+∂f2∂μ˙δμ˙=0.0.From the additive decomposition of strain in a framework of small strain theory, the following equation can be written:(19)σ=De:ε−De:εp,σ=De:ε−De:λ1∂g1∂σ.Now, taking differential variation of the following equation (19), we can get the value of δσ:(20)δσ=De:δε−De:δλ1∂g1∂σ+λ1∂2g1∂σ∂σδσ+λ1∂2g1∂σ∂κδκ+λ1∂2g1∂σ∂κ˙δκ˙.During the incremental elastic loading step, the elastoplastic integration of the incremental elastoplastic equation through return mapping total strain (ε) does not change. So the following equation (21) can be written:(21)δε=0.0.Putting equation (21) into equation (20), we get equations (22) and (24):(22)δσ+λ1De:∂2g1∂σ∂σδσ=−De:δλ1∂g1∂σ+λ1∂2g1∂σ∂κδκ+λ1∂2g1∂σ∂κ˙δκ˙,δσI+λ1De:∂2g1∂σ∂σ=−De:δλ1∂g1∂σ+λ1∂2g1∂σ∂κδκ+λ1∂2g1∂σ∂κ˙δκ˙,δσ=−De:I+λ1De:∂2g1∂σ∂σ−1δλ1∂g1∂σ+λ1∂2g1∂σ∂κδκ+λ1∂2g1∂σ∂κ˙δκ˙,δσ=−De−1+λ1∂2g1∂σ∂σ−1δλ1∂g1∂σ+λ1∂2g1∂σ∂κδκ+λ1∂2g1∂σ∂κ˙δκ˙.Similarly another equation (23) can be written also for the plastic potential in the compression direction, g2:(23)δσ=−De−1+λ1∂2g2∂σ∂σ−1δλ1∂g2∂σ+λ1∂2g2∂σ∂μδμ+λ1∂2g2∂σ∂μ˙δμ˙.Now, substituting the value ofδσ from equations (22) and (23) into lower part of equations (17) and (18), respectively, we get(24)f1+∂f1∂σδσ+∂f1∂κδk+∂f1∂κ˙δk˙=0.0,f1−∂f1∂σDe−1+λ1∂2g1∂σ∂σ−1δλ1∂g1∂σ+λ1∂2g1∂σ∂κδκ+λ1∂2g1∂σ∂κ˙δκ˙+∂f1∂κδk+∂f1∂κ˙δk˙=0.0,f2+∂f2∂σδσ+∂f2∂κδμ+∂f2∂κ˙δμ˙=0.0,f2−∂f2∂σDe−1+λ2∂2g2∂σ∂σ−1δλ2∂g2∂σ+λ2∂2g2∂σ∂μδμ+λ2∂2g2∂σ∂μ˙δμ˙+∂f2∂μδμ+∂f2∂μ˙δμ˙=0.0.Then, defining hardening modulusH as H1=−∂f1/∂κ;H2=−∂f2/∂μ and rate of hardening modulus Hr as H1r=−∂f1/∂κ˙;H2r=−∂f2/∂μ˙ and it is well-known that for plane strain case that δλ=δκandδλ˙=δκ˙, equation (24) can be rewritten as follows:(25)f1−∂f1∂σDe−1+λ1∂2g1∂σ∂σ−1δλ1∂g1∂σ+λ1∂2g1∂σ∂κδκ+λ1∂2g1∂σ∂κ˙δκ˙+∂f1∂κδk+∂f1∂κ˙δk˙=0,∂f1∂σDe−1+λ1∂2g1∂σ∂σ−1δλ1∂g1∂σ+λ1∂2g1∂σ∂κδλ+λ1∂2g1∂σ∂κ˙δκ˙+H1δλ+H1rδk˙=f1,δλ1=f1H1+H1r/Δt+∂f1/∂σDe−1+λ1∂2g1/∂σ∂σ−1∂g1/∂σ+λ1∂2g1/∂σ∂κ+λ1/Δt∂2g1/∂σ∂κ˙.Similarly, for the compression zone, the following equation can be written as follows:(26)δλ2=f2H2+H2r/Δt+∂f1/∂σDe−1+λ2∂2g1/∂σ∂σ−1∂g2/∂σ+λ2∂2g2/∂σ∂κ+λ2/Δt∂2g2/∂σ∂κ˙.Equation (11) can be rewritten in the following conventional way:(27)δλ1=f1H1+H1r/Δt+∂f1/∂σDe−1+λ1∂2g1/∂σ∂σ−1∂g1/∂σ+λ1∂2g1/∂σ∂κ+λ1/Δt∂2g1/∂σ∂κ˙,δλ1=f1H1+H1r/Δt+n1tD¯1em¯1,where(28)n1t=∂f1∂σ,D¯1e=D1e−1+λ1∂2g1∂σ∂σ−1,m¯1=∂g1∂σ+λ1∂2g1∂σ∂κ+λ1Δt∂2g1∂σ∂κ˙.Similarly, equation (26) for the compression zone can be also rewritten as(29)δλ2=f2H2+H2r/Δt+∂f2/∂σDe−1+λ2∂2g2/∂σ∂σ−1∂g2/∂σ+λ2∂2g2/∂σ∂κ+λ2/Δt∂2g2/∂σ∂κ˙,δλ2=f2H2+H2r/Δt+n2tD¯2em¯2,where(30)n2t=∂f2∂σ,D¯2e=D2e−1+λ2∂2g2∂σ∂σ−1,m¯2=∂g2∂σ+λ2∂2g2∂σ∂κ+λ2Δt∂2g2∂σ∂κ˙.Once, the increment of plasticity parameter,δλ is known for shear and compression direction, corrected stress can be calculated, and plasticity parameter can be updated as follows:(31)σn+1=σtrial−δλ1De∂g1∂σ,κn+1=κn+δλ1,σn+1=σtrial−δλ2De∂g1∂σ,μn+1=μn+δλ2.


## 4.2. Case II: Simultaneous Yielding

Case II describes the stress integration when both the yield surfaces are yielding at the same time. During the general loading, it is possible that both the yield functions are acting simultaneously. It represents a stress state lying inside the corner, bounded by the two-yield functions (Figure4). Following the same arguments as case I, the following equations can be written:(32)f˙1σn+1,κn+1,κ˙n+1=0.0,f1+∂f1∂σδσ+∂f1∂κδk+∂f1∂κ˙δk˙=0.0,(33)f˙2σn+1,μn+1,μ˙n+1=0.0,f2+∂f2∂σδσ+∂f2∂μδμ+∂f2∂μ˙δμ˙=0.0.From the additive decomposition of strain in a framework of small strain theory and Koiter's generalization of plastic strain generation, the following equation can be written:(34)σ=De:ε−De:εp,σ=De:ε−De:λ1∂g1∂σ−De:λ2∂g2∂σ.Now, taking differential variation of the following equation (34), the value of δσ can be deduced as follows:(35)δσ=De:δε−De:δλ1∂g1∂σ+λ1∂2g1∂σ∂σδσ+λ1∂2g1∂σ∂κδκ+λ1∂2g1∂σ∂κ˙δκ˙−De:δλ2∂g2∂σ+λ2∂2g2∂σ∂σδσ+λ2∂2g2∂σ∂μδμ+λ2∂2g1∂σ∂μ˙δμ˙.Usingδε=0.0 into equation (35), the following equation (36) can be written:(36)δσ+λ1De:∂2g1∂σ∂σδσ+λ2De:∂2g2∂σ∂σδσ=−De:δλ1∂g1∂σ+λ1∂2g1∂σ∂κδκ+λ1∂2g1∂σ∂κ˙δκ˙−De:δλ2∂g2∂σ+λ2∂2g2∂σ∂μδμ+λ1∂2g2∂σ∂μ˙δμ˙,δσ=−Rδλ1∂g1∂σ+λ1∂2g1∂σ∂κδκ+λ1∂2g1∂σ∂κ˙δκ˙+δλ2∂g2∂σ+λ2∂2g2∂σ∂μδμ+λ1∂2g2∂σ∂μ˙δμ˙,where R=De−1+λ1∂2g1/∂σ∂σ++λ2∂2g2/∂σ∂σ−1.Now substituting the value ofδσ from equation (36) into equation (32), the following equation can be derived:(37)f1−∂f1∂σRδλ1∂g1∂σ+λ1∂2g1∂σ∂κδκ+λ1∂2g1∂σ∂κ˙δκ˙+δλ2∂g2∂σ+λ2∂2g2∂σ∂μδμ+λ2∂2g2∂σ∂μ˙δμ˙+∂f1∂κδk+∂f1∂κ˙δk˙=0.0,∂f1∂κ+∂f1Δt∂κ˙+∂f1∂σR∂g1∂σ+∂f1∂σRλ1∂2g1∂σ∂κ+∂f1∂σRλ1Δt∂2g1∂σ∂κ˙δλ1+∂f1∂σR∂g2∂σ+∂f1∂σRλ2∂2g2∂σ∂μ+∂f1∂σRλ2Δt∂2g2∂σ∂μ˙δλ2=f1.Again substituting the value ofδσ from equation (36) into equation (33), the following equation can be derived:(38)f2−∂f2∂σRδλ1∂g1∂σ+λ1∂2g1∂σ∂κδκ+λ1∂2g1∂σ∂κ˙δκ˙+δλ2∂g2∂σ+λ2∂2g2∂σ∂μδμ+λ2∂2g2∂σ∂μ˙δμ˙+∂f2∂μδμ+∂f2∂μ˙δμ˙=0.0,∂f2∂μ+∂f2Δt∂μ˙+∂f2∂σR∂g2∂σ+∂f2∂σRλ2∂2g2∂σ∂μ+∂f2∂σRλ2Δt∂2g2∂σ∂μ˙δλ2+∂f2∂σR∂g1∂σ+∂f2∂σRλ1∂2g1∂σ∂κ+∂f2∂σRλ1Δt∂2g1∂σ∂κ˙δλ1=f2.Abbreviating and rearranging equations (37) and (38), the following equation (39) can be written:(39)β1δλ1+β2δλ2=f1,β3δλ1+β4δλ2=f2,where(40)β1=∂f1∂κ+∂f1Δt∂κ˙+∂f1∂σR∂g1∂σ+∂f1∂σRλ1∂2g1∂σ∂κ+∂f1∂σRλ1Δt∂2g1∂σ∂κ˙,β2=∂f1∂σR∂g2∂σ+∂f1∂σRλ2∂2g2∂σ∂μ+∂f1∂σRλ2Δt∂2g2∂σ∂μ˙,β3=∂f2∂σR∂g1∂σ+∂f2∂σRλ1∂2g1∂σ∂κ+∂f2∂σRλ1Δt∂2g1∂σ∂κ˙,β4=∂f2∂μ+∂f2Δt∂μ˙+∂f2∂σR∂g2∂σ+∂f2∂σRλ2∂2g2∂σ∂μ+∂f2∂σRλ2Δt∂2g2∂σ∂μ˙.Now the systems of simultaneous equations (39) can be solved, and δλ1 and δλ2 are determined as follows:(41)δλ1=β4f1−β2f2β1β4−β2β3,δλ2=β1f2−β3f1β1β4−β2β3,Once, the increment of plasticity parameters,δλ1 and δλ2, are known, the corrected stress can be calculated and plasticity parameter can be updated as follows:(42)σn+1=σtrial−δλ1De∂g1∂σ−δλ2De∂g2∂σ,κn+1=κn+δλ1,μn+1=μn+δλ2.


## 5. Consistent Modulus Derivation

In order to obtain the quadratic convergence in solving the nonlinear systems of simultaneous equation through finite element method, consistent or algorithmic modulus is an essential element. Here, the derivation of the consistent modulus has been divided into two parts as before ((1) case I: separate yielding; (2) case-II: simultaneous yielding)


### 5.1. Case I: Separate Yielding

#### 5.1.1. Shear Yielding

(43)
δ
σ
+
λ
1
D
e
:
∂
2
g
1
∂
σ
∂
σ
δ
σ
=
D
e
:
δ
ε
−
D
e
:
δ
λ
1
∂
g
1
∂
σ
+
λ
1
∂
2
g
1
∂
σ
∂
κ
δ
κ
+
λ
1
∂
2
g
1
∂
σ
∂
κ
˙
δ
κ
˙
,
δ
σ
I
+
λ
1
D
e
:
∂
2
g
1
∂
σ
∂
σ
=
D
e
:
δ
ε
−
D
e
:
δ
λ
1
∂
g
1
∂
σ
+
λ
1
∂
2
g
1
∂
σ
∂
κ
δ
κ
+
λ
1
∂
2
g
1
∂
σ
∂
κ
˙
δ
κ
˙
,
δ
σ
=
I
+
λ
1
D
e
:
∂
2
g
1
∂
σ
∂
σ
−
1
D
e
:
δ
ε
−
D
e
:
δ
λ
1
∂
g
1
∂
σ
+
λ
1
∂
2
g
1
∂
σ
∂
κ
δ
κ
+
λ
1
∂
2
g
1
∂
σ
∂
κ
˙
δ
κ
˙
,
δ
σ
=
D
e
−
1
+
λ
1
∂
2
g
1
∂
σ
∂
σ
−
1
δ
ε
−
δ
λ
1
∂
g
1
∂
σ
+
λ
1
∂
2
g
1
∂
σ
∂
κ
δ
κ
+
λ
1
∂
2
g
1
∂
σ
∂
κ
˙
δ
κ
˙
,
δ
σ
=
D
¯
1
e
δ
ε
−
δ
λ
1
∂
g
1
∂
σ
+
λ
1
∂
2
g
1
∂
σ
∂
κ
δ
κ
+
λ
1
∂
2
g
1
∂
σ
∂
κ
˙
δ
κ
˙
.Now, substituting the value ofδσ from equation (43) into lower part of equation (17), we get(44)f1+∂f1∂σδσ+∂f1∂κδk+∂f1∂κ˙δk˙=0.0,f1+∂f1∂σD¯1eδε−δλ1∂g1∂σ+λ1∂2g1∂σ∂κδκ+λ1∂2g1∂σ∂κ˙δκ˙+∂f1∂κδk+∂f1∂κ˙δk˙=0.0.Then, defining hardening modulusH as H=−∂f1/∂κ and rate hardening modulus Hr as Hr=−∂f1/∂κ˙ and as we know that for plane strain case that δλ=δκandδλ˙=δκ˙, we can rewrite equation (44) as follows:(45)f1+∂f1∂σD¯1eδε−δλ1∂g1∂σ+λ1∂2g1∂σ∂κδκ+λ1∂2g1∂σ∂κ˙δκ˙+∂f1∂κδk+∂f1∂κ˙δk˙=0.0,∂f1∂σD¯1eδλ1∂g1∂σ+λ1∂2g1∂σ∂κδκ+λ1∂2g1∂σ∂κ˙δκ˙+Hδλ+Hrδk˙=f1+∂f1∂σD¯1eδε,∂f1∂σD¯1eδλ∂g1∂σ+λ1∂2g1∂σ∂κδλ+λ1Δt∂2g1∂σ∂κ˙δλ+Hδλ+HrΔtδλ=f1+∂f1∂σD¯1eδε,δλ=f1+∂f1/∂σD¯1eδεH+Hr/Δt+∂f1/∂σD¯1e∂g1/∂σ+λ1∂2g1/∂σ∂κ+λ1/Δt∂2g1/∂σ∂κ˙.Equation (45) can be rewritten in the following conventional way:(46)δλ=f1+∂f1/∂σD¯1eδεH+Hr/Δt+∂f1/∂σD¯1e∂g1/∂σ+λ1∂2g1/∂σ∂κ+λ1/Δt∂2g1/∂σ∂κ˙,δλ=f1+nTD¯eδεH+Hr/Δt+nTD¯1em¯,where(47)n1T=∂f1∂σ,D¯1e=De−1+λ∂2g1∂σ∂σ−1,m¯1=∂g1∂σ+λ1∂2g1∂σ∂κ+λ1Δt∂2g1∂σ∂κ˙.Now, substituting the value ofδλ from equation (46) to the last step of equation (43), we get(48)δσ=D¯1eδε−D¯1em¯δλ,δσ=D¯1eδε−D¯1em¯f+nTD¯1eδεH+Hr/Δt+ntD¯1em¯.At the end of the load stepf≈0.0, then equation (48) will be as follows:(49)δσ=D¯1eδε−D¯1em¯nTD¯1eδεH+Hr/Δt+nTD¯1em¯,δσ=D¯1e−D¯1em¯nTD¯1eH+Hr/Δt+nTD¯1em¯δε,δσ=D¯1vpδε,where consistent tangent moduli, D¯1vp is defined as follows:(50)D¯1vp=D¯1e−D¯1em¯nTD¯1eH+Hr/Δt+nTD¯1em¯.


#### 5.1.2. Compression Yield

In a similar way as the shearing yield, the consistent modulus for the compression yield can be derived as(51)D¯2vp=D¯2e−D¯2em¯2n2TD¯2eHc+Hrc/Δt+n2TD¯2em¯2,where(52)n2T=∂f2∂σ,D¯2e=De−1+λ∂2g2∂σ∂σ−1,m¯2=∂g2∂σ+λ2∂2g2∂σ∂μ+λ2Δt∂2g2∂σ∂μ˙,Hc=−∂f2∂μ,Hrc=−∂f2∂μ˙.


### 5.2. Case II: Simultaneous Yielding

Now from equation (35) of stress integration by return mapping scheme, we get(53)δσ=De:δε−De:δλ1∂g1∂σ+λ1∂2g1∂σ∂σδσ+λ1∂2g1∂σ∂κδκ+λ1∂2g1∂σ∂κ˙δκ˙−De:δλ2∂g2∂σ+λ2∂2g2∂σ∂σδσ+λ2∂2g2∂σ∂μδμ+λ2∂2g1∂σ∂μ˙δμ˙.Now considering all the equations after the equilibrium iterations, strain increment does not vanish. Equation (53) can be rearranged into the following form:(54)δσ+λ1De:∂2g1∂σ∂σδσ+λ2De:∂2g2∂σ∂σδσ=De:δε−De:δλ1∂g1∂σ+λ1∂2g1∂σ∂κδκ+λ1∂2g1∂σ∂κ˙δκ˙−De:δλ2∂g2∂σ+λ2∂2g2∂σ∂μδμ+λ2∂2g2∂σ∂μ˙δμ˙,δσ=De−1+λ1∂2g1∂σ∂σ++λ2∂2g2∂σ∂σ−1δε−δλ1∂g1∂σ+λ1∂2g1∂σ∂κδκ+λ1∂2g1∂σ∂κ˙δκ˙−δλ2∂g2∂σ+λ2∂2g2∂σ∂μδμ+λ2∂2g2∂σ∂μ˙δμ˙,δσ=Rδε−δλ1∂g1∂σ+λ1∂2g1∂σ∂κδκ+λ1∂2g1∂σ∂κ˙δκ˙−δλ2∂g2∂σ+λ2∂2g2∂σ∂μδμ+λ2∂2g2∂σ∂μ˙δμ˙,where R=De−1+λ1∂2g1/∂σ∂σ++λ2∂2g2/∂σ∂σ−1.Substituting the value ofδσ in equation (54) into equations (32) and (33), the following equations can be derived:(55)f1+∂f1∂σRδε−δλ1∂g1∂σ+λ1∂2g1∂σ∂κδκ+λ1∂2g1∂σ∂κ˙δκ˙−δλ2∂g2∂σ+λ2∂2g2∂σ∂μδμ+λ2∂2g2∂σ∂μ˙δμ˙+∂f1∂κδk+∂f1∂κ˙δk˙=0.0,∂f1∂κ+∂f1Δt∂κ˙+∂f1∂σR∂g1∂σ+∂f1∂σRλ1∂2g1∂σ∂κ+∂f1∂σRλ1Δt∂2g1∂σ∂κ˙δλ1+∂f1∂σR∂g2∂σ+∂f1∂σRλ2∂2g2∂σ∂μ+∂f1∂σRλ2Δt∂2g2∂σ∂μ˙δλ2=f1+∂f1∂σRδε,f2+∂f2∂σRδε−δλ1∂g1∂σ+λ1∂2g1∂σ∂κδκ+λ1∂2g1∂σ∂κ˙δκ˙−δλ2∂g2∂σ+λ2∂2g2∂σ∂μδμ+λ2∂2g2∂σ∂μ˙δμ˙+∂f2∂μδμ+∂f2∂μ˙δμ˙=0.0,∂f2∂μ+∂f2Δt∂μ˙+∂f2∂σR∂g2∂σ+∂f2∂σRλ2∂2g2∂σ∂μ+∂f2∂σRλ2Δt∂2g2∂σ∂μ˙δλ2+∂f2∂σR∂g1∂σ+∂f2∂σRλ1∂2g1∂σ∂κ+∂f2∂σRλ1Δt∂2g1∂σ∂κ˙δλ1=f2+∂f2∂σRδε.Using the same abbreviations of equation (39), the following equations can be derived:(56)β1δλ1+β2δλ2=f1+∂f1∂σRδε,β3δλ1+β4δλ2=f2+∂f2∂σRδε.Solving the 2 × 2 systems of simultaneous equations, following (δλ1,δλ2) variables can be obtained:(57)δλ1=β4f1+∂f1/∂σRdε−β2f2+∂f2/∂σRdεβ1β4−β2β3,δλ2=β1f2+∂f2/∂σRdε−β3f1+∂f1/∂σRdεβ1β4−β2β3.Now at the end of a time step, when convergence happens,f1=f2≈0.0, so equation (57) turns out to be as follows:(58)δλ1=β4∂f1/∂σRdε−β2∂f2/∂σRdεβ1β4−β2β3,δλ2=β1∂f2/∂σRdε−β3∂f1/∂σRdεβ1β4−β2β3.Now inputting the value ofδλ1 and δλ2 from equation (58) into equation (54), we find the consistent modulus (for plane strain case, δλ1=δκandδλ˙1=δκ˙ and δλ2=μandδλ˙2=δμ˙)(59)δσ=R−β4∂f1/∂σRR∂g1/∂σ−β2∂f2/∂σRR∂g1/∂σβ1β4−β2β3+λ1β4∂f1/∂σRR∂2g1/∂σ∂κ−β2∂f2/∂σRR∂2g1/∂σ∂κβ1β4−β2β3+λ1Δtβ4∂f1/∂σRR∂2g1/∂σ∂κ˙−β2∂f2/∂σRR∂2g1/∂σ∂κ˙β1β4−β2β3−β1∂f2/∂σRR∂g2/∂σ−β3∂f1/∂σRR∂g2/∂σβ1β4−β2β3+λ2β1∂f2/∂σRR∂2g2/∂σ∂μ−β3∂f1/∂σRR∂2g2/∂σ∂μβ1β4−β2β3+λ2Δtβ1∂f2/∂σRR∂2g2/∂σ∂μ˙−β3∂f1/∂σRR∂2g2/∂σ∂μ˙β1β4−β2β3dε.Details of the compression yield/plastic potential function is as follows:(60)f2=g2=−I−I0σ,μ,μ˙,I0=αμ1+gc,∂g2∂σ=∂f2∂σ=−1,−1,−1,0.All the 2nd derivative of the above equation is zero or null:(61)∂f2∂μ=−α1+gc,∂f2∂μ˙=−αμ∂gc∂μ˙.


## 5.1. Case I: Separate Yielding

### 5.1.1. Shear Yielding

(43)
δ
σ
+
λ
1
D
e
:
∂
2
g
1
∂
σ
∂
σ
δ
σ
=
D
e
:
δ
ε
−
D
e
:
δ
λ
1
∂
g
1
∂
σ
+
λ
1
∂
2
g
1
∂
σ
∂
κ
δ
κ
+
λ
1
∂
2
g
1
∂
σ
∂
κ
˙
δ
κ
˙
,
δ
σ
I
+
λ
1
D
e
:
∂
2
g
1
∂
σ
∂
σ
=
D
e
:
δ
ε
−
D
e
:
δ
λ
1
∂
g
1
∂
σ
+
λ
1
∂
2
g
1
∂
σ
∂
κ
δ
κ
+
λ
1
∂
2
g
1
∂
σ
∂
κ
˙
δ
κ
˙
,
δ
σ
=
I
+
λ
1
D
e
:
∂
2
g
1
∂
σ
∂
σ
−
1
D
e
:
δ
ε
−
D
e
:
δ
λ
1
∂
g
1
∂
σ
+
λ
1
∂
2
g
1
∂
σ
∂
κ
δ
κ
+
λ
1
∂
2
g
1
∂
σ
∂
κ
˙
δ
κ
˙
,
δ
σ
=
D
e
−
1
+
λ
1
∂
2
g
1
∂
σ
∂
σ
−
1
δ
ε
−
δ
λ
1
∂
g
1
∂
σ
+
λ
1
∂
2
g
1
∂
σ
∂
κ
δ
κ
+
λ
1
∂
2
g
1
∂
σ
∂
κ
˙
δ
κ
˙
,
δ
σ
=
D
¯
1
e
δ
ε
−
δ
λ
1
∂
g
1
∂
σ
+
λ
1
∂
2
g
1
∂
σ
∂
κ
δ
κ
+
λ
1
∂
2
g
1
∂
σ
∂
κ
˙
δ
κ
˙
.Now, substituting the value ofδσ from equation (43) into lower part of equation (17), we get(44)f1+∂f1∂σδσ+∂f1∂κδk+∂f1∂κ˙δk˙=0.0,f1+∂f1∂σD¯1eδε−δλ1∂g1∂σ+λ1∂2g1∂σ∂κδκ+λ1∂2g1∂σ∂κ˙δκ˙+∂f1∂κδk+∂f1∂κ˙δk˙=0.0.Then, defining hardening modulusH as H=−∂f1/∂κ and rate hardening modulus Hr as Hr=−∂f1/∂κ˙ and as we know that for plane strain case that δλ=δκandδλ˙=δκ˙, we can rewrite equation (44) as follows:(45)f1+∂f1∂σD¯1eδε−δλ1∂g1∂σ+λ1∂2g1∂σ∂κδκ+λ1∂2g1∂σ∂κ˙δκ˙+∂f1∂κδk+∂f1∂κ˙δk˙=0.0,∂f1∂σD¯1eδλ1∂g1∂σ+λ1∂2g1∂σ∂κδκ+λ1∂2g1∂σ∂κ˙δκ˙+Hδλ+Hrδk˙=f1+∂f1∂σD¯1eδε,∂f1∂σD¯1eδλ∂g1∂σ+λ1∂2g1∂σ∂κδλ+λ1Δt∂2g1∂σ∂κ˙δλ+Hδλ+HrΔtδλ=f1+∂f1∂σD¯1eδε,δλ=f1+∂f1/∂σD¯1eδεH+Hr/Δt+∂f1/∂σD¯1e∂g1/∂σ+λ1∂2g1/∂σ∂κ+λ1/Δt∂2g1/∂σ∂κ˙.Equation (45) can be rewritten in the following conventional way:(46)δλ=f1+∂f1/∂σD¯1eδεH+Hr/Δt+∂f1/∂σD¯1e∂g1/∂σ+λ1∂2g1/∂σ∂κ+λ1/Δt∂2g1/∂σ∂κ˙,δλ=f1+nTD¯eδεH+Hr/Δt+nTD¯1em¯,where(47)n1T=∂f1∂σ,D¯1e=De−1+λ∂2g1∂σ∂σ−1,m¯1=∂g1∂σ+λ1∂2g1∂σ∂κ+λ1Δt∂2g1∂σ∂κ˙.Now, substituting the value ofδλ from equation (46) to the last step of equation (43), we get(48)δσ=D¯1eδε−D¯1em¯δλ,δσ=D¯1eδε−D¯1em¯f+nTD¯1eδεH+Hr/Δt+ntD¯1em¯.At the end of the load stepf≈0.0, then equation (48) will be as follows:(49)δσ=D¯1eδε−D¯1em¯nTD¯1eδεH+Hr/Δt+nTD¯1em¯,δσ=D¯1e−D¯1em¯nTD¯1eH+Hr/Δt+nTD¯1em¯δε,δσ=D¯1vpδε,where consistent tangent moduli, D¯1vp is defined as follows:(50)D¯1vp=D¯1e−D¯1em¯nTD¯1eH+Hr/Δt+nTD¯1em¯.


### 5.1.2. Compression Yield

In a similar way as the shearing yield, the consistent modulus for the compression yield can be derived as(51)D¯2vp=D¯2e−D¯2em¯2n2TD¯2eHc+Hrc/Δt+n2TD¯2em¯2,where(52)n2T=∂f2∂σ,D¯2e=De−1+λ∂2g2∂σ∂σ−1,m¯2=∂g2∂σ+λ2∂2g2∂σ∂μ+λ2Δt∂2g2∂σ∂μ˙,Hc=−∂f2∂μ,Hrc=−∂f2∂μ˙.


## 5.1.1. Shear Yielding

(43)
δ
σ
+
λ
1
D
e
:
∂
2
g
1
∂
σ
∂
σ
δ
σ
=
D
e
:
δ
ε
−
D
e
:
δ
λ
1
∂
g
1
∂
σ
+
λ
1
∂
2
g
1
∂
σ
∂
κ
δ
κ
+
λ
1
∂
2
g
1
∂
σ
∂
κ
˙
δ
κ
˙
,
δ
σ
I
+
λ
1
D
e
:
∂
2
g
1
∂
σ
∂
σ
=
D
e
:
δ
ε
−
D
e
:
δ
λ
1
∂
g
1
∂
σ
+
λ
1
∂
2
g
1
∂
σ
∂
κ
δ
κ
+
λ
1
∂
2
g
1
∂
σ
∂
κ
˙
δ
κ
˙
,
δ
σ
=
I
+
λ
1
D
e
:
∂
2
g
1
∂
σ
∂
σ
−
1
D
e
:
δ
ε
−
D
e
:
δ
λ
1
∂
g
1
∂
σ
+
λ
1
∂
2
g
1
∂
σ
∂
κ
δ
κ
+
λ
1
∂
2
g
1
∂
σ
∂
κ
˙
δ
κ
˙
,
δ
σ
=
D
e
−
1
+
λ
1
∂
2
g
1
∂
σ
∂
σ
−
1
δ
ε
−
δ
λ
1
∂
g
1
∂
σ
+
λ
1
∂
2
g
1
∂
σ
∂
κ
δ
κ
+
λ
1
∂
2
g
1
∂
σ
∂
κ
˙
δ
κ
˙
,
δ
σ
=
D
¯
1
e
δ
ε
−
δ
λ
1
∂
g
1
∂
σ
+
λ
1
∂
2
g
1
∂
σ
∂
κ
δ
κ
+
λ
1
∂
2
g
1
∂
σ
∂
κ
˙
δ
κ
˙
.Now, substituting the value ofδσ from equation (43) into lower part of equation (17), we get(44)f1+∂f1∂σδσ+∂f1∂κδk+∂f1∂κ˙δk˙=0.0,f1+∂f1∂σD¯1eδε−δλ1∂g1∂σ+λ1∂2g1∂σ∂κδκ+λ1∂2g1∂σ∂κ˙δκ˙+∂f1∂κδk+∂f1∂κ˙δk˙=0.0.Then, defining hardening modulusH as H=−∂f1/∂κ and rate hardening modulus Hr as Hr=−∂f1/∂κ˙ and as we know that for plane strain case that δλ=δκandδλ˙=δκ˙, we can rewrite equation (44) as follows:(45)f1+∂f1∂σD¯1eδε−δλ1∂g1∂σ+λ1∂2g1∂σ∂κδκ+λ1∂2g1∂σ∂κ˙δκ˙+∂f1∂κδk+∂f1∂κ˙δk˙=0.0,∂f1∂σD¯1eδλ1∂g1∂σ+λ1∂2g1∂σ∂κδκ+λ1∂2g1∂σ∂κ˙δκ˙+Hδλ+Hrδk˙=f1+∂f1∂σD¯1eδε,∂f1∂σD¯1eδλ∂g1∂σ+λ1∂2g1∂σ∂κδλ+λ1Δt∂2g1∂σ∂κ˙δλ+Hδλ+HrΔtδλ=f1+∂f1∂σD¯1eδε,δλ=f1+∂f1/∂σD¯1eδεH+Hr/Δt+∂f1/∂σD¯1e∂g1/∂σ+λ1∂2g1/∂σ∂κ+λ1/Δt∂2g1/∂σ∂κ˙.Equation (45) can be rewritten in the following conventional way:(46)δλ=f1+∂f1/∂σD¯1eδεH+Hr/Δt+∂f1/∂σD¯1e∂g1/∂σ+λ1∂2g1/∂σ∂κ+λ1/Δt∂2g1/∂σ∂κ˙,δλ=f1+nTD¯eδεH+Hr/Δt+nTD¯1em¯,where(47)n1T=∂f1∂σ,D¯1e=De−1+λ∂2g1∂σ∂σ−1,m¯1=∂g1∂σ+λ1∂2g1∂σ∂κ+λ1Δt∂2g1∂σ∂κ˙.Now, substituting the value ofδλ from equation (46) to the last step of equation (43), we get(48)δσ=D¯1eδε−D¯1em¯δλ,δσ=D¯1eδε−D¯1em¯f+nTD¯1eδεH+Hr/Δt+ntD¯1em¯.At the end of the load stepf≈0.0, then equation (48) will be as follows:(49)δσ=D¯1eδε−D¯1em¯nTD¯1eδεH+Hr/Δt+nTD¯1em¯,δσ=D¯1e−D¯1em¯nTD¯1eH+Hr/Δt+nTD¯1em¯δε,δσ=D¯1vpδε,where consistent tangent moduli, D¯1vp is defined as follows:(50)D¯1vp=D¯1e−D¯1em¯nTD¯1eH+Hr/Δt+nTD¯1em¯.


## 5.1.2. Compression Yield

In a similar way as the shearing yield, the consistent modulus for the compression yield can be derived as(51)D¯2vp=D¯2e−D¯2em¯2n2TD¯2eHc+Hrc/Δt+n2TD¯2em¯2,where(52)n2T=∂f2∂σ,D¯2e=De−1+λ∂2g2∂σ∂σ−1,m¯2=∂g2∂σ+λ2∂2g2∂σ∂μ+λ2Δt∂2g2∂σ∂μ˙,Hc=−∂f2∂μ,Hrc=−∂f2∂μ˙.


## 5.2. Case II: Simultaneous Yielding

Now from equation (35) of stress integration by return mapping scheme, we get(53)δσ=De:δε−De:δλ1∂g1∂σ+λ1∂2g1∂σ∂σδσ+λ1∂2g1∂σ∂κδκ+λ1∂2g1∂σ∂κ˙δκ˙−De:δλ2∂g2∂σ+λ2∂2g2∂σ∂σδσ+λ2∂2g2∂σ∂μδμ+λ2∂2g1∂σ∂μ˙δμ˙.Now considering all the equations after the equilibrium iterations, strain increment does not vanish. Equation (53) can be rearranged into the following form:(54)δσ+λ1De:∂2g1∂σ∂σδσ+λ2De:∂2g2∂σ∂σδσ=De:δε−De:δλ1∂g1∂σ+λ1∂2g1∂σ∂κδκ+λ1∂2g1∂σ∂κ˙δκ˙−De:δλ2∂g2∂σ+λ2∂2g2∂σ∂μδμ+λ2∂2g2∂σ∂μ˙δμ˙,δσ=De−1+λ1∂2g1∂σ∂σ++λ2∂2g2∂σ∂σ−1δε−δλ1∂g1∂σ+λ1∂2g1∂σ∂κδκ+λ1∂2g1∂σ∂κ˙δκ˙−δλ2∂g2∂σ+λ2∂2g2∂σ∂μδμ+λ2∂2g2∂σ∂μ˙δμ˙,δσ=Rδε−δλ1∂g1∂σ+λ1∂2g1∂σ∂κδκ+λ1∂2g1∂σ∂κ˙δκ˙−δλ2∂g2∂σ+λ2∂2g2∂σ∂μδμ+λ2∂2g2∂σ∂μ˙δμ˙,where R=De−1+λ1∂2g1/∂σ∂σ++λ2∂2g2/∂σ∂σ−1.Substituting the value ofδσ in equation (54) into equations (32) and (33), the following equations can be derived:(55)f1+∂f1∂σRδε−δλ1∂g1∂σ+λ1∂2g1∂σ∂κδκ+λ1∂2g1∂σ∂κ˙δκ˙−δλ2∂g2∂σ+λ2∂2g2∂σ∂μδμ+λ2∂2g2∂σ∂μ˙δμ˙+∂f1∂κδk+∂f1∂κ˙δk˙=0.0,∂f1∂κ+∂f1Δt∂κ˙+∂f1∂σR∂g1∂σ+∂f1∂σRλ1∂2g1∂σ∂κ+∂f1∂σRλ1Δt∂2g1∂σ∂κ˙δλ1+∂f1∂σR∂g2∂σ+∂f1∂σRλ2∂2g2∂σ∂μ+∂f1∂σRλ2Δt∂2g2∂σ∂μ˙δλ2=f1+∂f1∂σRδε,f2+∂f2∂σRδε−δλ1∂g1∂σ+λ1∂2g1∂σ∂κδκ+λ1∂2g1∂σ∂κ˙δκ˙−δλ2∂g2∂σ+λ2∂2g2∂σ∂μδμ+λ2∂2g2∂σ∂μ˙δμ˙+∂f2∂μδμ+∂f2∂μ˙δμ˙=0.0,∂f2∂μ+∂f2Δt∂μ˙+∂f2∂σR∂g2∂σ+∂f2∂σRλ2∂2g2∂σ∂μ+∂f2∂σRλ2Δt∂2g2∂σ∂μ˙δλ2+∂f2∂σR∂g1∂σ+∂f2∂σRλ1∂2g1∂σ∂κ+∂f2∂σRλ1Δt∂2g1∂σ∂κ˙δλ1=f2+∂f2∂σRδε.Using the same abbreviations of equation (39), the following equations can be derived:(56)β1δλ1+β2δλ2=f1+∂f1∂σRδε,β3δλ1+β4δλ2=f2+∂f2∂σRδε.Solving the 2 × 2 systems of simultaneous equations, following (δλ1,δλ2) variables can be obtained:(57)δλ1=β4f1+∂f1/∂σRdε−β2f2+∂f2/∂σRdεβ1β4−β2β3,δλ2=β1f2+∂f2/∂σRdε−β3f1+∂f1/∂σRdεβ1β4−β2β3.Now at the end of a time step, when convergence happens,f1=f2≈0.0, so equation (57) turns out to be as follows:(58)δλ1=β4∂f1/∂σRdε−β2∂f2/∂σRdεβ1β4−β2β3,δλ2=β1∂f2/∂σRdε−β3∂f1/∂σRdεβ1β4−β2β3.Now inputting the value ofδλ1 and δλ2 from equation (58) into equation (54), we find the consistent modulus (for plane strain case, δλ1=δκandδλ˙1=δκ˙ and δλ2=μandδλ˙2=δμ˙)(59)δσ=R−β4∂f1/∂σRR∂g1/∂σ−β2∂f2/∂σRR∂g1/∂σβ1β4−β2β3+λ1β4∂f1/∂σRR∂2g1/∂σ∂κ−β2∂f2/∂σRR∂2g1/∂σ∂κβ1β4−β2β3+λ1Δtβ4∂f1/∂σRR∂2g1/∂σ∂κ˙−β2∂f2/∂σRR∂2g1/∂σ∂κ˙β1β4−β2β3−β1∂f2/∂σRR∂g2/∂σ−β3∂f1/∂σRR∂g2/∂σβ1β4−β2β3+λ2β1∂f2/∂σRR∂2g2/∂σ∂μ−β3∂f1/∂σRR∂2g2/∂σ∂μβ1β4−β2β3+λ2Δtβ1∂f2/∂σRR∂2g2/∂σ∂μ˙−β3∂f1/∂σRR∂2g2/∂σ∂μ˙β1β4−β2β3dε.Details of the compression yield/plastic potential function is as follows:(60)f2=g2=−I−I0σ,μ,μ˙,I0=αμ1+gc,∂g2∂σ=∂f2∂σ=−1,−1,−1,0.All the 2nd derivative of the above equation is zero or null:(61)∂f2∂μ=−α1+gc,∂f2∂μ˙=−αμ∂gc∂μ˙.


## 6. Model Implementation and Calibration

The above numerical model is coded in FORTRAN-90 as a generic “usermat,” which works in all three major finite element packages “Abaqus,” “Plaxis 2D,” and “Midas GTS.” The results are found to be same for all three cases. Here in this study, only the results from Plaxis2D are reported. Model calibration is usually a trial and error process. After a comprehensive trial and error process, the following basic elastoplastic parameters (Table2) are found to calibrate the model and simulate the existing data from the triaxial tests carried out in Section 2. The simulations are carried out in 8-element plane strain idealization.Table 2
Model parameters.
Parameters of softrock model
Type of parameter
Typical values
Young’s modulusE (kPa)
Value
800000.0
Poisson’s ratio
Value
0.2
Elastic model
Type
1
Cyclic model
Type
2
Plastic model
Type
5
R
max
Compression
0.001
R
r
Compression
1.6
k
f
Compression
1.1
R
max
Extension
0.0005
R
r
Extension
1.15
k
f
Extension
1.01
Bcont
All
0.3
Ccont
All
0.1
Α
Isotach
1.15
M
Isotach
0.035
ε′ir-ref
Isotach
0.1e − 9
Cohesion (kN/m2)
Value
1230
Friction
Degree
12
Dilatancy parameter
Value
900
Dilatancy cutoff
Value
0.10
Unit weight,γt, (kN/m3)
Value
19It is interesting to note that, to make the stress-dilatancy rule to become more pressure dependent, a term of (900/σm) is added to equation (6). So here, 900 is the pressure dependency parameter for the calibration of the volumetric strain response of this material. It was necessary to adjust the volumetric plastic strain to the actual experimental value.


## 7. Results and Discussion

The result shown here in this section is a measure of the capability of this model. The effect of loading rate is seen in strength and deformation of soft rock. The results as depicted in Figures5–11 show these features in detail. Figure 5 shows the variation of deviatoric stress versus axial strain for different confining pressures, while keeping the loading strain constant (0.008%/minute). From Figure 5, it can be observed that slope of the ascending loading curve is almost same for all the different confining pressures. Just before the peak, the curves start to deviate and go for a reduced peak value with the decrease in confining stresses as seen in the experimental results. Figure 6 shows the variation of volumetric strain (%) versus axial strain (%) for various confining stresses for a single strain rate (0.008%/minute). Here, the material shows more contractive behavior as the confining stress increases. Similar behavior is observed in the experiment as well. Figures 7 and 8 show the behavior of this model in case of UC type and plane strain type analysis for various strain rates. For both the cases, strain rates are varied from 0.0015%/min to 1.5%, which is about 1000 times. In both the cases, confining stresses are kept constant to be 10 kPa. In both the cases, it can be seen that the ultimate strength of the material increases with the increase in rate of loading. The slope of the loading curve also steepens with the increase in the rate of loading although it was not visible in naked eye. In order to make it visible, Figures 7(b) and 8(b) are zoomed and the trend was clearly visible. Figure 9 shows the cyclic capability of the model and shows the variation of deviatoric stress versus maximum shear strain for different strain rates. The strain rates are varied from 0.0015%/min to 1.5%, which is about 1000 times. The confining pressure was kept constant at 10 kPa. The effect of strain rate is similar to the monotonic loading case in compression. It can also be observed from Figure 9 that effect of strain rate in extension is less than the compression as has been seen in the experiments. Figure 10 shows the variation of deviatoric stress versus axial strain for a full cycle of loading with stress relaxation in both loading and unloading parts of the response curve. The simulation was done at a constant confining stress of 10 kPa and strain rate of 0.0015%/minute. It can be observed from the figure that relaxation was perfectly carried out in both the compression and extension parts of the cyclic loading. Figure 11 shows the variation of deviatoric stress versus axial strain for a full cycle of loading achieved through load control (not displacement control). The simulation was done at a constant confining stress of 10 kPa and strain rate of 1.5%/minute. It can be observed that the cycles change directions from loading to unloading and unloading to loading in a more smooth kind of movement as observed in the experiments. The smooth rounded turn around of the deviatoric versus axial strain curve is due to the creep-effect at the peak of the curve. When the loading curve starts to move for the unloading in high rate of loading (1.5%/min in this case), load itself does not change much, only the rate starts to change to the lower value, thus lowering the deviatoric stress, until the rate of loading changes sign (from “+ve” loading rate to “−ve” loading rate). This rather complex mechanism of cyclic loading under high loading rate is completely captured by this model quite well.Figure 5
Deviatoric stress versus axial strain (%) for various confining stress for a single strain rate. Displacement control analysis is done.Figure 6
Volumetric strain (%) versus axial strain (%) for various confining stress for a single strain rate. Displacement control analysis is done.Figure 7
Deviatoric stress versus axial strain (%) for various strain rates under a single confining stress of 10 kPa. (a) Displacement control analysis is done. (b) Displacement control analysis is done in the compression direction (zoomed up figure of Figure 7(a)).
(a)
(b)Figure 8
Deviatoric stress versus axial strain (%) for various strain rates under a single confining stress of 10 kPa. (a) Displacement control analysis is done. (b) Displacement control analysis is done in the extension direction (zoomed up fig. of Figure 8(a)).
(a)
(b)Figure 9
Deviatoric stress versus maximum shear strain for various strain rates under a single confining stress of 10 kPa. Displacement control cyclic analysis is done in both compression and extension directions.Figure 10
Deviatoric stress versus axial strain for a single strain rate under a single confining stress of 10 kPa. Displacement control cyclic relaxation analysis is done in both compression and extension directions.Figure 11
Deviatoric stress versus axial strain for single strain rate under a single confining stress of 10 kPa. Load control cyclic analysis is done in both compression and extension directions.Figures12(a)–12(c) show the results of simulation of deviatoric stress versus axial strain for different strain rates and three different confining stresses 100 kN/m2, 400 kN/m2, and 800 kN/m2, respectively. It can be observed that the analytical solution by the double hardening model in red color fits quite well with the experimental data for all three confining stresses.Figure 12
Deviatoric stress versus axial strain of the simulation of the CD test by the double hardening model for multiple strain rates under confining stresses of (a) 100 kPa, (b) 400 kPa, and (c) 800 kPa.
(a)
(b)
(c)


## 8. Conclusions

In this paper, a time-dependent double hardening model has been formulated and a new algorithm of the stress-integration of incremental stress-strain equations is proposed. A series of triaxial compression tests were carried out at different strain rates to calculate the model parameters. The model showed the following capabilities:(1)
The model is pressure dependent(2)
The model is loading-rate dependent(3)
The model is capable of time-dependent cyclic loading(4)
The model is capable of producing creep and stress relaxation anywhere in the load/unload curve



---

*Source: 2927178-2019-10-07.xml*