Research Article  Open Access
Mohammed Saiful Alam Siddiquee, Amin Hamdi, "A TimeDependent Double Hardening Model for Soft Rock", Advances in Civil Engineering, vol. 2019, Article ID 2927178, 18 pages, 2019. https://doi.org/10.1155/2019/2927178
A TimeDependent Double Hardening Model for Soft Rock
Abstract
The timedependent 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 threecomponent elastoviscoplastic framework is used to develop a timedependent double hardening model to predict the behavior of soft rock both in compression and shear. Due to the limitation of timedependent single yield hardening model in predicting the behavior of soft rock in compression direction (volumetric deformation) of loading, another timedependent yield surface is added to the compression direction. The intersection of twoyield surface usually gives rise to singularity phenomenon, which is avoided by using Koiter’s generalization principle. A new formulation of the stressintegration of incremental stressstrain 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 stressstrain 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 stressstrain tests on sand under drained conditions; i.e., effects of strain rate and its change on the stressstrain relation, creep deformation, and stress relaxation during otherwise monotonic loading (ML) at a constant strain rate.
Within the framework of the general nonlinear threecomponent model (Figure 1), Benedetto et al. [5] and Tatsuoka et al. [7] proposed a set of stressstrain models to simulate the effects of material viscosity on the stressstrain 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 stressstrain model called the “new isotach” was proposed to describe the loading rate effects of claylike 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 threecomponent 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 ratedependent yield surface and the consistency condition is fulfilled. In this proposal, the consistency method is used to integrate the incremental viscoelastoplastic 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 elastoviscoplastic 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 strainratedependent 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 xaxis as equivalent normal stress (as in Figure 3), the yield function and plastic potential of the elastoviscoplastic 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 stressdilatancy 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 crossanisotropic 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 crossanisotropic 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 twoyield 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 ratedependent 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 wellknown 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 twoyield 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) caseII: 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 FORTRAN90 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 8element plane strain idealization.

It is interesting to note that, to make the stressdilatancy 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 strai