Elastic-viscoplastic deformation models of salt rocks

Results of secondary creep tests for salt rock specimens were mathematically described. Elastic-viscoplastic modeling of viscous properties of salt rocks was based on elastoplastic models using a non-associated Mohr-Coulomb plastic flow and an associated volumetric yield criterion. In combination with these models the one-parameter Bingham and Duvaut-Lions as well as the two-parameter Perzyna and Perić viscoplasticity laws were analyzed. It’s been found that all laws being considered reflect an increase in longitudinal deformation of salt specimens during creep. It should be noted though that the yield criteria combined with the Perić law are characterized by relatively simple, stable parametric support and acceptable description of steady-state creep. The use of viscoplastic deformation model makes it possible to point out the fracture regions formed during creep according to a specified strength criterion.


INTRODUCTION
alt rocks are classified as geomaterials which are characterized by strongly pronounced rheological properties [1,2,3].Development of potash and salt deposits is generally carried out using room-and-pillar mining system.The overlying strata are supported by rib pillars.The rib pillars are subjected to constant, almost time invariable load, which causes the creep activation [4].In this case the geomaterial deformation depends not only on the level of applied load, but also on the time of its action.The low rate of mining advance practically eliminates the transient creep phenomenon (primary creep) in rib pillars.In this regard, only the secondary creep was analyzed.The study of the deformation process of large cubic salt specimens under uniaxial compression [5] made it possible to identify fracture (yield) criteria and associated with them plastic flow rules, which most accurately describe both their loading diagram and transverse-longitudinal deformations.In other words, the suitable yield functions and plastic potentials that adequately simulate the behavior of salt rocks under time-independent load were selected.These turned out to be the nonassociated Mohr-Coulomb criterion [6] and the associated volumetric strength criterion of rocks [7].

S
Creep is one of the phenomena of viscous properties of a material [8].Material viscosity is commonly studied by the standard laboratory relaxation and creep tests [9,10].The results are used to calibrate viscoelastic and/or viscoplastic models.As part of the research on transient effects being appeared during deformation of salt rocks, a series of creep tests under uniaxial compression were performed.Red sylvinite specimens of small size (60×30×30 mm) was used, as shown in Fig. 1.
The creep behavior of specimens was studied at various load levels (ratio of an effective stress  to an ultimate uniaxial compression strength  c -p =  / c ) [11].The results are demonstrated in Fig. 2. A steady-state creep (secondary creep) moment is assumed as the reference point.The purpose of the current research is the generalization of elastoplastic deformation model of salt rocks by including a viscous component and to describe the results of creep tests.Selected based on the results [5] the non-associated Mohr-Coulomb criterion and the associated volumetric criterion were used as yield criteria.

MATHEMATICAL MODEL
ransient effects in the deformation behavior of salt rocks were described using an isotropic viscoplastic model of media [12,13,14].Elastic straining was governed by Hooke's law.The basic constitutive relation of the general elastic-viscoplastic model describing the stress-strain state of a material is written in rate form as vp ( ) D : where D e is the fourth-order elastic stiffness tensor, the colon sign : denotes a double contraction, and the inelastic strains T vp ( , , ), ( , ) 0 0, ( , ) 0 Here,  and  are the yield function and the plastic flow potential, respectively.In contrast to the elastoplastic model, the plastic multiplier time derivative is now the function of stress tensor and material internal state parameters-A (plastic) and P (viscous).Therefore, this may violate the Kuhn-Tucker condition, i.e. the stress point could be outside the yield surface.
To simplify the model construction, the sets of internal parameters are assumed constant (A = const and P = const).Thus, hardening and temperature dependence of the parameters P were not considered.
The yield function and plastic potential for the non-associated Mohr-Coulomb criterion in the principal stress space are written as The corresponding sets of internal state parameters where c is the cohesion,  is the frictional angle, and  is the dilatancy angle.For the associated volumetric criterion of rock strength, the yield function and the plastic potential are defined as follows: , The internal state parameters vol { , } Here  c and  t are the uniaxial compression strength and uniaxial tensile strength, respectively.
The differential equations of the elastic-viscoplastic model were solved using the displacement-based finite element method.The geometry of salt specimens was represented by a finite element mesh of 8-node isoparametric hexahedral elements with 8 integration points.The solution domain (60×30×30 mm) was meshed by cubic elements with 1 mm side.The constructed model of salt specimens' deformation was calibrated according to the experimental data using multivariant modeling by varying only parameters P.
The integration schemes are required to be more accurate both at the global (within the time step) and local (at the integration point) levels to account viscoplastic properties of the geomaterial.Thus, the construction of the elasticviscoplastic model was carried out using the following integration schemes.An automatic Newton-Raphson scheme with substepping and error control [15] was used to implement the global time integration.For the local integration procedure of elastic-viscoplastic relations the implicit Euler scheme of the return-mapping algorithm was used [12,13,14].The multisurface representation of yield surface and plastic potential describing the evolution of viscoplastic strains was modeled using the Koiter's generalization [16].However, the choice of the plastic potential or linear combination of them was based on the elastoplastic solution only.Thus, during the local integration procedure the elastoplastic problem is solved first.Then, its solution is used to define which surface/edge of the overall yield surface the stress point is related to.After that the corresponding equation or the combination of equations are included in the local integration system of viscoplastic relations.
The local integration of the constitutive relations was carried out in the principal stress space using the spectral decomposition of symmetric tensor [12].The implicit Euler scheme of the return mapping algorithm reduces the integration to the system In most cases, system (7) is nonlinear.Since the sets of internal parameters are assumed constant, A n = A and P n = P.The solution of such nonlinear systems with respect to  n and  was carried out using the Newton-Raphson method.According to it, the system ( 7) is represented by residuals in matrix form (Voigt notation) as trial D N( , ) 0 ( , , ) 0 where D e is now the elastic stiffness matrix.The iterative Newton-Raphson process continues until the following condition is satisfied where j denotes the number of its iteration.Here STOL is the specified tolerance of the residuals' norm from zero.The number of Newton-Raphson iterations is limited by the MAXITL value.If there is no solution within MAXITL iterations, then the global solution is considered unreachable, the time increment is reduced and the global iterative process with the updated increment is started over.The overall integration scheme contains five tune parameters of the automatic integration algorithm for a nonlinear system: three are at the global level (DTOL, ITOL and MAXIT-described in [15]) and two (STOL and MAXITL-described above) are at the local level.The choice of their values largely depends on a particular problem and specifies the accuracy and efficiency of the solution.In this paper, the following values were used: DTOL = 10 -3 , ITOL = 10 -5 , STOL = 10 -6 , MAXIT = 10, MAXITL = 10.

VISCOPLASTICITY LAWS
Bingham's law he classical and simplest law of viscoplasticity-Bingham's law [12]-is linear with respect to the yield function The viscosity  is the only parameter of model, whose dimension is the product of stress and time, e.g., MPasec.Thus, the set of parameters P consists of a single element

Non-associated Mohr-Coulomb criterion
The residual form of local integration system for the non-associated Mohr-Coulomb criterion depends on which surface/edge the stress point is related to.
Working in the principal stress space and considering the sextant  1 ≥  2 ≥  3 -Fig.3, hereinafter the tensile stress is implied to be positive-with the stress point being related to the main surface  1 , the matrix form of residuals' system is as follows for which the corresponding Jacobian is where I is the unit matrix 3×3, and P 1 = D e N 1 is the projection vector, which is constant for the Mohr-Coulomb criterion.
It should be noted that unlike the elastoplastic model, the residual R  in (12) describes a "dynamic" yield surface the stresses are projected onto as a result of the return mapping algorithm.Such surface shape depends not only on the current stress state and internal parameters, but also on the strain rate.For the yield surface  1 , the plastic flow N 1 and the normal Ñ 1equations are written as in [5]  In other case, the residuals ( 12) are complemented by one equation and the residual R  is modified as long as the stresses are related to an edge of the yield surface where indices 2 and 6 denote the number of the adjacent surface and the corresponding plastic flow, depending on the edge the stress point is related to- 12 (extension edge) or  16 (compression edge).The expressions of  2 and  6 in the principal stress space take the form: ) s i n 2c o s N 0,1 sin , 1 sin , N 0,1 sin , 1 sin N 1 sin , 1 sin , 0 , N 1 sin , 1 sin , 0 The Jacobian of the system ( 15) is written as Note that ( 12) and ( 15) are linear, therefore a single iteration is required for the local iterative Newton-Raphson process to be converged.

Associated volumetric criterion
The volumetric criterion is continuous with respect to the stress tensor [7].This implies the system of residuals for the associated plastic flow rule and viscoplastic Bingham's law to be of two equations The plastic flow direction in R  now depends on the stress tensor and the material internal parameters.The Jacobian is written as where the plastic flow vector and its derivative in the principal stress space are defined as In the context of Bingham's law and the associated plastic potential of the volumetric criterion the results obtained by multivariant numerical simulation of the creep are presented as diagrams and calibrated model parameters in Fig. 4   Here C e = (D e ) -1 is the elastic compliance tensor,  is the relaxation time, and ̂ is the backbone stress, which corresponds to the nearest projection of the stress tensor on the yield surface.Thus, the only one element is included in the set As the Duvaut-Lions viscoplastic law is of different type, the corresponding system of residuals is different from (8).Thus, its general form in matrix notation is written as: trial trial ˆˆ0 ˆD N( , ) 0 ( , ) 0 The first equation is responsible for the Duvaut-Lions law, which is obtained by substituting t̇ vp for the plastic strain tensor N( n , A) in ( 8) 1 (here and below the subscript denotes the equation ordinal from top to bottom/left to right).The second equation describes the backbone stress and the third equation makes this stress to be on the yield surface.The corresponding Jacobian is where the elements with the circumflex sign ˆ are referred to the backbone stress on the yield surface (plastic response of a material).So, the projection vector and the normal to the yield surface are expressed as

Non-associated Mohr-Coulomb criterion
As before, for the Mohr-Coulomb criterion in the principal stress space we consider the following cases-the stress point is related to the surface and to one of the edges.In the first case, the system of residuals generally corresponds to (23 The only difference is the plastic flow.Here it is a constant.Thus, the system is linear and its solution could be obtained in a closed form.The local integration scheme is reduced to sequential calculations: In the case of the stress point being related to an edge, the R  in the system (23) is replaced similarly to (15) 1 and one of the equations of adjacent surface is added The obtained system is solved in a closed form where the coefficients are defined as The final stress tensor is calculated via (27) 3 .The results of multivariant numerical simulations of the creep in salt specimens obtained using the Duvaut-Lions law and the non-associated Mohr-Coulomb criterion are presented in Fig. 5.The calibrated parameters of the elastic-viscoplastic model for each numerical experiment are given in Table 3.

Associated volumetric criterion
As the volumetric criterion is nonlinear with respect to the stress tensor, the system of residuals is nonlinear as well.The system of residuals and the corresponding Jacobian are represented completely as the general form (23)-(24).The system of residuals could be reduced by removing one equation in order to get the solution more efficient.Thus, an expression for n could be obtained from R  and used as a function.So then it is introduced to other equations of the system.As a result, the system of residuals for the Duvaut-Lions viscoplastic model is reduced to two equations.
where the backbone stress function is trial trial ˆ( , , , ) The Jacobian is written more concisely as The following results were obtained via the multivariant numerical simulation of the creep in salt specimens using the Duvaut-Lions viscoplastic law and the associated volumetric criterion.The creep diagrams at various load levels are presented in Fig. 5.The calibrated parameters of the elastic-viscoplastic model for each numerical experiment are given in Table 4.

Perzyna's law
One of the most widely used laws in computational viscoplasticity was proposed by Perzyna [12,13,14].It could be expressed as 1 1 1 where  is the viscosity-related parameter of time dimension,  e is an equivalent stress,  y is the corresponding yield point, and m is the rate-sensitivity/strain rate hardening parameter.Typically, the von Mises effective stress is used as  e .In this case the yield point  y is the ultimate von Mises effective stress.The associated von Mises criterion-also is referred to as the Prandtl-Reiss law [12]-is commonly used as the plastic potential.
However, any yield criterion could be reduced to (34) by assuming as  e the yield function with its free term being added, and taking the free term of the yield function as  y .Given this, the Perzyna viscoplastic model could be written in terms of the yield function where ( , ) ( , ) ( ) and  y (A) is the free term of the yield function representing a certain limit of equivalent stress.
The set of viscoplasticity parameters now consists of two elements The parameter m specifies the material reaction to a strain rate change.The closer it is to zero, the more the material behavior is consistent with plastic straining.The limit m → 0 is equivalent to a perfect plastic behavior, i.e. the transient effects are completely absent.
The constraint equation of the plastic multiplier in the system of residuals according to the Perzyna viscoplastic law is written as follows: ( , ) ( ) 0 The main drawback of the Perzyna model also follows from equation (38).In the limit m = 0 the model doubles the yield point  e (, A) = 2 y (A).

Non-associated Mohr-Coulomb criterion
As before, we consider the cases when the stress tensor is related to the surface and when to one of the edges.In the first case the system of residuals in matrix notation takes the form: where the yield point  y for the Mohr-Coulomb criterion is written as The corresponding Jacobian is Note that the diagonal element in the second row of the Jacobian cannot be determined from the initial approximation  = 0.In order to eliminate this computational uncertainty, we rearrange the equation R  as follows Now the Jacobian is represented as Such form of the tangent operator of the system of residuals is possible to be determined from the initial approximation  = 0. Nevertheless, it should be noted that the Jacobian is undetermined on the yield surface.Thus, the solution of the system might be unstable as the stress tensor approaches the yield surface as well as at small values of the parameter m.In the second case when the stress tensor is related to an edge of the Mohr-Coulomb yield surface, the system of residuals is written in a similar way ) The corresponding Jacobian is The results of multivariant numerical simulation of the creep of salt specimens are shown in Fig. 6.The corresponding calibrated parameters of the elastic-viscoplastic model are given in Table 5.

Load level
Young

Associated volumetric criterion
For the associated volumetric criterion and the Perzyna viscoplastic law, the local integration scheme in matrix form is written as Here the plastic flow is no longer constant unlike the Mohr-Coulomb plastic potential.The yield point for the volumetric criterion is given as vol ( ) The corresponding Jacobian takes the following form: The results of numerical creep experiments for salt specimens are shown in Fig. 6.The corresponding calibrated parameters of the elastic-viscoplastic model for each experiment are given in Table 6.

Load level
Young's modulus, GPa Poisson's ratio

Perić's law
The main drawback of the Perzyna viscoplasticity law is the double yield point with transient effects being completely absent in the limit m → 0. In addition, small values of the rate-sensitivity (strain rate hardening) parameter could lead to an unstable solution.The Perić model [8] eliminates these drawbacks.In the ANSYS engineering software package, the model is referred to as the Pierce model [17].Its analytical form could be written as 1 1 1 The model parameters are the same as (37).In contrast to the Perzyna's law, here the normalized effective stress is raised to a power.The expression for  e could be represented in terms of the yield function

  
In view of (50), equation (49) for a discrete time increment could be written in residual form as Clearly, equation ( 51) is recovered the yield surface equation in the limit m = 0.

Non-associated Mohr-Coulomb criterion
Considering the case when the stress tensor is related to the main surface of the yield surface, the system of residuals is written as Also, the corresponding Jacobian is It is clear from this that the Jacobian is determined for all values of the viscoplastic parameters, as well as for the initial approximation of the local Newton-Raphson iterative process.
In the case when the stress tensor is related to one of the edges of the overall Mohr-Coulomb yield surface, the system of residuals in the principal stress space for the sextant  1 ≥  2 ≥  3 is written as The Jacobian of the system is as follows The results of multivariant creep simulation of salt specimens are presented as diagrams in Fig. 7.The calibrated parameters of the elastic-viscoplastic "non-associated Mohr-Coulomb + Perić" model for each numerical experiment are given in Table 7. Load

Associated volumetric criterion
Considering the volumetric criterion as a yield surface and plastic potential, the local system of the implicit Euler integration scheme of the return-mapping algorithm is given as The corresponding Jacobian is The results of multivariant numerical experiments of the creep at various load levels along with the results of laboratory tests are shown in Fig. 7.The corresponding parameters obtained due to model calibration are presented in Table 8.
In contrast to the viscoelastic deformation model, the viscoplastic one makes it possible to point out the fracture regions formed in the material according to a specified strength criterion.The fracture due to shear and tear of a salt specimen along its vertical section at different load levels are illustrated in Fig. 8.The action of all negative (compressive) principal stresses ( 1 < 0,  2 < 0,  3 < 0) reached the yield point was considered as a fracture due to shear while the fracture due to tear was taken place at  1 > 0. In the numerical implementation, if at more than half of integration points of a finite element the yield condition was met-with the corresponding stresses, it was considered to be fractured.It is seen from the figure that in the central section of the specimen, fracture due to tear is dominant.Shear fracture is localized in the upper and lower halves of the specimen.The volume of shear fracture at load level of 0.6 in a given section is lower than at level of 0.8.With increasing load level, the shape of shear fracture regions is changed symmetrical to the center of the specimen along its section.So, for the load level of 0.6 the shear fracture regions are almost plane shaped in the cross-section while for the load level of 0.8 they are shaped as semicircles.The regions of fracture due to tear formed in here are almost identical for both load levels.

Figure 2 :
Figure 2: Averaged steady-state creep diagrams for red sylvinite specimens at various load levels.

Figure 3 :
Figure 3: Multisurface representation of the Mohr-Coulomb yield surface in the deviatoric plane.

Figure 4 :
Figure 4: The results of creep simulation at various load levels-Bingham's lawDuvaut-Lions lawProposed by Duvaut and Lions[14] an alternative viscoplasticity law is interesting in that it is an extension of the elastoplastic model and the only one parameter is involved-the relaxation time in its pure form.The Duvaut-Lions law-just as Bingham's law-is linear, yet relative to the backbone stress, and is written as follows

Figure 5 :
Figure 5: The results of creep simulation at various load levels-Duvaut-Lions law

Figure 6 :
Figure 6: The results of creep simulation at various load levels-Perzyna's law.

Figure 7 :Figure 8 :
Figure 7: The results of creep simulation at various load levels-Perić's law

Table 1 :
The results of creep simulation using Bingham's law and the non-associated Mohr-Coulomb plastic potential are presented in Fig.4.The model parameters obtained by multivariant simulation are given in Table1.