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

#### 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 (Figure 1), 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).

In [5] and [7], it has been showed that at least three different functional forms of the viscous component, , 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 is a nonlinear function of instantaneous value of while it is always proportional to the instantaneous value of as follows:where is the viscosity function, which is always zero or positive and given as follows for any strain () or stress () history (with or without cyclic loading):where is the absolute value of and *α*, , and *m* are positive material constants. According to this model, as far as ML continues along a fixed stress path, the viscous stress component, , is a unique function of instantaneous values of and , 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 ) is a function of instantaneous strain rate, , not , while, with the new isotach model, is a function of . This difference results into significant differences in the model behaviour, in particular during stress relaxation with 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.

**(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 Table 1.

#### 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 Figure 3, the definition of normal stresses can be redefined as follows:

**(a)**

**(b)**

Using the extended *x*-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:wherewhere *I*_{1} is the 1st stress invariant, *J*_{2} is the 2nd invariant of deviatoric stress, *K* is the cohesion terms (used only in Figure 3(b)), *α*(*θ*) is the Lode angle function, and is the mobilized angle of internal friction for the modified normal stress definition. In this way, is defined as follows:

The original Rowe’s stress-dilatancy relation was modified to define the plastic potential as stated as follows:whereAnd, has a lower limit as follows: if , .

This dilatancy cutoff is very important for large strain computation in the softening regions, where is the mobilized angle of friction, is the residual angle of internal friction, is the mobilized angle of dilatancy, and is the mean pressure in kPa.

The evolution of the yield function is modeled as described in:where *dε*^{p} = [2 (*de*^{P}_{11})^{2} + 2 (*de*^{P}_{22})^{2} + 2 (*de*^{P}_{33})^{2} + (*dγ*^{P}_{11})^{2} +(*dγ*^{P}_{23})^{2} + (*dγ*^{P}_{31})^{2}]^{1/2}*de*^{P}_{ij} ∼ *dγ*^{P}_{ij} are the deviatoric components of plastic strain increments and *κ* = plasticity parameter (internal variable).

The evolution equation is defined by the following sets of equations:

If then,

If then,

##### 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 bywhere *I*_{1} is the 1st stress invariant of deviatoric stresses and *I*_{0} is the evolution parameter in the compression direction given bywhere is the material constant signifying the rate of movement of *I*_{0} in the compression direction and 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:”where is described in equation (10) and 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

##### 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:where *E*_{1} = “ or *E*_{h} at or *σ*_{h} = 1.0 kN/m^{2}.”

and *E*_{h} are isotropic ( = *E*_{h} = *E*_{p}), when = *σ*_{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:

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

##### 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:

From the additive decomposition of strain in a framework of small strain theory, the following equation can be written:

Now, taking differential variation of the following equation (19), we can get the value of :

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:

Putting equation (21) into equation (20), we get equations (22) and (24):

Similarly another equation (23) can be written also for the plastic potential in the compression direction, :

Now, substituting the value of from equations (22) and (23) into lower part of equations (17) and (18), respectively, we get

Then, defining hardening modulus *H* as and rate of hardening modulus as and it is well-known that for plane strain case that , equation (24) can be rewritten as follows:

Similarly, for the compression zone, the following equation can be written as follows:

Equation (11) can be rewritten in the following conventional way:where

Similarly, equation (26) for the compression zone can be also rewritten aswhere

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:

##### 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 (Figure 4). Following the same arguments as case I, the following equations can be written:

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:

Now, taking differential variation of the following equation (34), the value of can be deduced as follows:

Using into equation (35), the following equation (36) can be written:where .

Now substituting the value of from equation (36) into equation (32), the following equation can be derived:

Again substituting the value of from equation (36) into equation (33), the following equation can be derived:

Abbreviating and rearranging equations (37) and (38), the following equation (39) can be written:where

Now the systems of simultaneous equations (39) can be solved, and and are determined as follows:

Once, the increment of plasticity parameters, and , are known, the corrected stress can be calculated and plasticity parameter can be updated as follows:

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

Now, substituting the value of from equation (43) into lower part of equation (17), we get

Then, defining hardening modulus *H* as and rate hardening modulus as and as we know that for plane strain case that , we can rewrite equation (44) as follows:

Equation (45) can be rewritten in the following conventional way:where

Now, substituting the value of from equation (46) to the last step of equation (43), we get

At the end of the load step , then equation (48) will be as follows:where consistent tangent moduli, is defined as follows:

###### 5.1.2. Compression Yield

In a similar way as the shearing yield, the consistent modulus for the compression yield can be derived aswhere

##### 5.2. Case II: Simultaneous Yielding

Now from equation (35) of stress integration by return mapping scheme, we get

Now considering all the equations after the equilibrium iterations, strain increment does not vanish. Equation (53) can be rearranged into the following form:where .

Substituting the value of in equation (54) into equations (32) and (33), the following equations can be derived:

Using the same abbreviations of equation (39), the following equations can be derived:

Solving the 2 × 2 systems of simultaneous equations, following () variables can be obtained:

Now at the end of a time step, when convergence happens, , so equation (57) turns out to be as follows:

Now inputting the value of and from equation (58) into equation (54), we find the consistent modulus (for plane strain case, and )

Details of the compression yield/plastic potential function is as follows:

All the 2nd derivative of the above equation is zero or null:

#### 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 (Table 2) 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.

It 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 Figures 5–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.

**(a)**

**(b)**

**(a)**

**(b)**

Figures 12(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/m^{2}, 400 kN/m^{2}, and 800 kN/m^{2}, 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.

**(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

#### Data Availability

No data were used to support this study.

#### Disclosure

The research paper was not funded by the King Abdulaziz University. It is carried out as part of the employment of King Abdulaziz University.

#### Conflicts of Interest

The authors declare that they have no conflicts of interest.