Microsoft Word - numero_70_art_11_5113.docx A. Baryakh et alii, Frattura ed Integrità Strutturale, 70(2024) 191-209; DOI: 10.3221/IGF-ESIS.70.11 192 Elastic-viscoplastic deformation models of salt rocks A. Baryakh, A. Tsayukov Mining Institute of the Ural Branch of the Russian Academy of Sciences, Russia bar@mi-perm.ru, http://orcid.org/0000-0003-2737-6166 andrei.tsayukov@mi-perm.ru, http://orcid.org/0000-0002-9982-4776 KEYWORDS. Salt rocks, Rib pillars, Elastic-viscoplastic model, Yield criteria, Creep, Mathematical modeling. Citation: Baryakh, A., Tsayukov, A., Elastic- viscoplastic deformation models of salt rocks, Frattura ed Integrità Strutturale, 70 (2024) 191-209. Received: 30.07.2024 Accepted: 15.08.2024 Published: 23.08.2024 Issue: 10.2024 Copyright: © 2024 This is an open access article under the terms of the CC-BY 4.0, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited. 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 non- associated Mohr-Coulomb criterion [6] and the associated volumetric strength criterion of rocks [7]. S https://youtu.be/unHOmhNinCE A. Baryakh et alii, Frattura ed Integrità Strutturale, 70(2024) 191-209; DOI: 10.3221/IGF-ESIS.70.11 193 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. Figure 1: Red sylvinite specimens. 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. Figure 2: Averaged steady-state creep diagrams for red sylvinite specimens at various load levels. 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 :e     (1) where De is the fourth-order elastic stiffness tensor, the colon sign : denotes a double contraction, and the inelastic strains T A. Baryakh et alii, Frattura ed Integrità Strutturale, 70(2024) 191-209; DOI: 10.3221/IGF-ESIS.70.11 194 vp ( , , ), ( , ) 0 0, ( , ) 0 f          A Ρ A A          (2) 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       max min max min, ( )sin 2 cos , , c           A A A          (3) The corresponding sets of internal state parameters { , }, { , }c c  A A  (4) 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:       2 2 2 vol 1 2 2 3 3 1 1 2 3 vol vol 1 , ( ) ( ) ( ) ( )( ) 2 , , c t c t                 A A A                 (5) The internal state parameters vol { , }c tA   (6) 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 elastic- viscoplastic 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 A. Baryakh et alii, Frattura ed Integrità Strutturale, 70(2024) 191-209; DOI: 10.3221/IGF-ESIS.70.11 195 trial D : N( , ) ( , , ) e n n n n n nf t         A A Ρ       (7) In most cases, system (7) is nonlinear. Since the sets of internal parameters are assumed constant, An = A and Pn = 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 e n n nf t              R A R A Ρ        (8) where De is now the elastic stiffness matrix. The iterative Newton-Raphson process continues until the following condition is satisfied   2 , STOL T j R R  (9) 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 1 ( , )  A   (10) 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  Ρ  (11) 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 T A. Baryakh et alii, Frattura ed Integrità Strutturale, 70(2024) 191-209; DOI: 10.3221/IGF-ESIS.70.11 196 trial 1 1 D N 0 ( , ) 0 e n n t               R R A        (12) for which the corresponding Jacobian is 1 1 P NT t         I J   (13) where I is the unit matrix 3×3, and P1 = De N1 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 N1 and the normal Ñ1— equations are written as in [5]       1 1 3 1 3 1 1 , ( )sin 2 cos N 1 sin , 0, 1 sin N 1 sin , 0, 1 sin T T c              A             (14) Figure 3: Multisurface representation of the Mohr-Coulomb yield surface in the deviatoric plane. 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   1 2 ,6 trial 1 1 2,6 2,6 1 1 2,6 2,6 D N N 0 ( , ) 0 ( , ) 0 e n n n t t                          R R A R A             (15) 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: A. Baryakh et alii, Frattura ed Integrità Strutturale, 70(2024) 191-209; DOI: 10.3221/IGF-ESIS.70.11 197             2 2 3 2 3 6 1 2 1 2 2 2 6 6 , ( )sin 2 cos , ( )sin 2 cos N 0, 1 sin , 1 sin , N 0, 1 sin , 1 sin N 1 sin , 1 sin , 0 , N 1 sin , 1 sin , 0 T T T T c c                               A A                         (16) The Jacobian of the system (15) is written as 1 2,6 1 2,6 P P N 0 N 0 T T t t             I J     (17) Note that (12) and (15) are linear, therefore a single iteration is required for the local iterative Newton-Raphson process to be converged. 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 Table 1. Load level Young's modulus, GPa Poisson's ratio Cohesion, MPa Frictional angle, degree. Dilatancy angle, degree Viscosity, MPahour 0.3 1.5 0.3 1.7 30 18 35 0.4 1.5 0.3 1.7 30 18 35 0.5 1.5 0.3 1.7 30 18 20 0.6 1.5 0.3 1.7 30 18 20 0.7 1.5 0.3 1.7 30 18 10 0.8 1.5 0.3 1.7 30 18 6 Table 1: “Mohr-Coulomb + Bingham” model parameters 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 trial vol vol D N( , ) 0 ( , ) 0 e n n n t               R A R A         (18) The plastic flow direction in R now depends on the stress tensor and the material internal parameters. The Jacobian is written as N D P N e T t           I J    (19) where the plastic flow vector and its derivative in the principal stress space are defined as A. Baryakh et alii, Frattura ed Integrità Strutturale, 70(2024) 191-209; DOI: 10.3221/IGF-ESIS.70.11 198 1 2 3 2 1 3 3 1 2 2 1 12 N 1 2 1N 2 , 1 1 22 c t c t c t                                                (20) 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 and Table 2, respectively. Load level Young's modulus, GPa Poisson's ratio Uniaxial tensile strength, MPa Uniaxial compressive strength, MPa Viscosity, MPahour 0.3 1.5 0.3 1 5 35 0.4 1.5 0.3 1 5 50 0.5 1.5 0.3 1 5 45 0.6 1.5 0.3 1 5 60 0.7 1.5 0.3 1 5 50 0.8 1.5 0.3 1 5 45 Table 2: “Volumetric + Bingham” model parameters Figure 4: The results of creep simulation at various load levels—Bingham's law Duvaut-Lions law Proposed 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 vp 1 ˆC : ( )e     (21) Here Ce = (De)-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  Ρ  (22) 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: A. Baryakh et alii, Frattura ed Integrità Strutturale, 70(2024) 191-209; DOI: 10.3221/IGF-ESIS.70.11 199 trial trial ˆ ˆ 0 ˆ ˆD N( , ) 0 ˆ( , ) 0 n n e n n n t t                     R R A R A               (23) 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 N ˆD P ˆ N̂ 0 e T t t                   I I 0 J 0 I 0      (24) 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 ˆ ˆˆ ˆP D N( , ), N N( , )e A A  (25) 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) trial trial ˆ 1 1 ˆ 0 ˆ D N 0 ˆ( , ) 0 n n e n n t t                       R R R A              (26) 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: trial 1 1 1 trial 1 trial ( , ) N D N ˆ D N ˆ e e n n n t t t                     A           (27) 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 2,6 2,6 ˆ( , ) 0n   R A  (28) The obtained system is solved in a closed form A. Baryakh et alii, Frattura ed Integrità Strutturale, 70(2024) 191-209; DOI: 10.3221/IGF-ESIS.70.11 200   trial trial 1 2,6 1 trial trial 2,6 1 2,6 trial 1 1 2,6 2,6 ( , ) ( , ) ( , ) ( , ) ˆ D N Ne n a b ad bc a c ad bc                      A A A A           (29) where the coefficients are defined as 1 2,61 1 2,6 1 2,6 2,6 N D NN D N , N D N , N D N ee e e ba c d       (30) 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. Load level Young's modulus, GPa Poisson's ratio Cohesion, MPa Frictional angle, degree. Dilatancy angle, degree Relaxation time, hour 0.3 1.5 0.3 1.7 30 18 9 0.4 1.5 0.3 1.7 30 18 9 0.5 1.5 0.3 1.7 30 18 4.6 0.6 1.5 0.3 1.7 30 18 4.6 0.7 1.5 0.3 1.7 30 18 2.4 0.8 1.5 0.3 1.7 30 18 1.5 Table 3: “Mohr-Coulomb + Duvaut-Lions” model parameters. 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.     trial vol vol ˆD N , 0 ˆ , 0 e n n n t t                R A R A         (31) where the backbone stress function is trial trialˆ ( , , , )n n n t t t t         P       (32) The Jacobian is written more concisely as A. Baryakh et alii, Frattura ed Integrità Strutturale, 70(2024) 191-209; DOI: 10.3221/IGF-ESIS.70.11 201 N ˆD P ˆ N̂ 0 e T t t t t                       I J     (33) 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. Load level Young's modulus, GPa Poisson's ratio Uniaxial tensile strength, MPa Uniaxial compressive strength, MPa Relaxation time, hour 0.3 1.5 0.3 1 5 13 0.4 1.5 0.3 1 5 13 0.5 1.5 0.3 1 5 8 0.6 1.5 0.3 1 5 8 0.7 1.5 0.3 1 5 5 0.8 1.5 0.3 1 5 3 Table 4: “Volumetric + Duvaut-Lions” model parameters Figure 5: The results of creep simulation at various load levels—Duvaut-Lions law 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 m e y             (34) 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. A. Baryakh et alii, Frattura ed Integrità Strutturale, 70(2024) 191-209; DOI: 10.3221/IGF-ESIS.70.11 202 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 1 ( , )1 ( ) m y         A A     (35) where ( , ) ( , ) ( )e n y  A A A    (36) 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  ,mΡ  (37) 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 m n y t        R A A    (38) 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: trial 1 1 D N 0 ( , ) ( ) 0 e n m n y t                  R R A A         (39) where the yield point  y for the Mohr-Coulomb criterion is written as ( ) 2 cosy c A  (40) The corresponding Jacobian is 1 1 P N m yT m t              I J     (41) 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 A. Baryakh et alii, Frattura ed Integrità Strutturale, 70(2024) 191-209; DOI: 10.3221/IGF-ESIS.70.11 203 trial 1 1 1 D N 0 ( , ) 0 ( ) e n m n y t                       R A R A         (42) Now the Jacobian is represented as 1 1 1 1 1 P 1 N m T ym t                 I J    (43) 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   1 2 ,6 trial 1 1 2,6 2,6 1 1 1 1 2,6 2,6 D N N 0 ( , ) 0 ( ) ( , ) 0 ( ) e n m n y m n y t t                                         R A R A A R A               (44) The corresponding Jacobian is 1 2,6 1 1 1 1 1 2,6 2,6 2,6 P P 1 N 0 1 N 0 m T y m T y m t m t                             I J       (45) 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's modulus, GPa Poisson's ratio Cohesion, MPa Frictional angle, degree. Dilatancy angle, degree Viscosity, hour Rate- sensitivity 0.3 1.5 0.3 1.7 30 18 3.5·103 0.4 0.4 1.5 0.3 1.7 30 18 104 0.4 0.5 1.5 0.3 1.7 30 18 104 0.4 0.6 1.5 0.3 1.7 30 18 1.6·104 0.4 0.7 1.5 0.3 1.7 30 18 2.2·104 0.3 0.8 1.5 0.3 1.7 30 18 2.8·104 0.3 Table 5: “Mohr-Coulomb + Perzyna” model parameters A. Baryakh et alii, Frattura ed Integrità Strutturale, 70(2024) 191-209; DOI: 10.3221/IGF-ESIS.70.11 204 Associated volumetric criterion For the associated volumetric criterion and the Perzyna viscoplastic law, the local integration scheme in matrix form is written as trial vol 1 vol vol D N( , ) 0 ( , ) 0 ( ) e n n m n y t                     R A A R A          (46) 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( )y c tA   (47) The corresponding Jacobian takes the following form: 1 N D P 1 N e m T ym t                I J     (48) 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 Uniaxial tensile strength, MPa Uniaxial compressive strength, MPa Viscosity, hour Rate- sensitivity 0.3 1.5 0.3 1 5 3.0·104 0.4 0.4 1.5 0.3 1 5 1.6·105 0.4 0.5 1.5 0.3 1 5 3.4·105 0.4 0.6 1.5 0.3 1 5 8.6·105 0.4 0.7 1.5 0.3 1 5 1.1·106 0.4 0.8 1.5 0.3 1 5 1.6·106 0.4 Table 6: “Volumetric + Perzyna” model parameters. Figure 6: The results of creep simulation at various load levels—Perzyna’s law. A. Baryakh et alii, Frattura ed Integrità Strutturale, 70(2024) 191-209; DOI: 10.3221/IGF-ESIS.70.11 205 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 m e y                (49) 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 ( , ) ( , ) ( )e y A A A    (50) In view of (50), equation (49) for a discrete time increment could be written in residual form as ( , ) ( ) 1 0 m y t t               R A A     (51) 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 trial 1 1 D N 0 ( , ) ( ) 1 0 e n m n y t t                          R R A A          (52) Also, the corresponding Jacobian is 1 1 1 P N m yT m t t t                 I J     (53) 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   1 2 ,6 trial 1 1 2,6 2,6 1 1 2,6 2,6 D N N 0 ( , ) ( ) 1 0 ( , ) ( ) 1 0 e n m n y m n y t t t t                                                  R R A A R A A                (54) A. Baryakh et alii, Frattura ed Integrità Strutturale, 70(2024) 191-209; DOI: 10.3221/IGF-ESIS.70.11 206 The Jacobian of the system is as follows 1 2,6 1 1 1 1 2,6 2,6 P P N 0 N 0 m yT m yT m t t t m t t t                              I J         (55) 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 level Young's modulus, GPa Poisson's ratio Cohesion, MPa Frictional angle, degree. Dilatancy angle, degree Viscosity, hour Rate- sensitivity 0.3 1.5 0.3 1.7 30 18 105 0.25 0.4 1.5 0.3 1.7 30 18 105 0.3 0.5 1.5 0.3 1.7 30 18 105 0.28 0.6 1.5 0.3 1.7 30 18 105 0.3 0.7 1.5 0.3 1.7 30 18 105 0.28 0.8 1.5 0.3 1.7 30 18 105 0.28 Table 7: “Mohr-Coulomb + Perić” model parameters 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 trial vol vol vol D N( , ) 0 ( , ) ( ) 1 0 e n n m n y t t                         R A R A A           (56) The corresponding Jacobian is 1 N D P N e m yT m t t t                  I J       (57) 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 A. Baryakh et alii, Frattura ed Integrità Strutturale, 70(2024) 191-209; DOI: 10.3221/IGF-ESIS.70.11 207 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. Load level Young's modulus, GPa Poisson's ratio Uniaxial tensile strength, MPa Uniaxial compressive strength, MPa Viscosity, hour Rate- sensitivity 0.3 1.5 0.3 1 5 2.2·105 0.3 0.4 1.5 0.3 1 5 1.9·106 0.3 0.5 1.5 0.3 1 5 5.9·106 0.3 0.6 1.5 0.3 1 5 2.1·107 0.3 0.7 1.5 0.3 1 5 3.9·107 0.3 0.8 1.5 0.3 1 5 6.8·107 0.3 Table 8: “Volumetric + Perić” model parameters Figure 7: The results of creep simulation at various load levels— Perić’s law Figure 8: Fracture along the vertical central section of the specimen due to shear ( ) and tear ( ) in accordance with the associated volumetric yield criterion and the Perić viscoplasticity law at load levels: a) 0.6 and b) 0.8. A. Baryakh et alii, Frattura ed Integrità Strutturale, 70(2024) 191-209; DOI: 10.3221/IGF-ESIS.70.11 208 CONCLUSION s part of the research on deformation nature of salt rocks with time the results of laboratory creep tests were mathematically described. Creep tests were carried out using red sylvinite specimens. Viscous properties of salt rocks were described using elastic-viscoplastic models as an extension of the previously constructed elastoplastic deformation models of salt specimens for the non-associated Mohr-Coulomb plastic flow and the associated volumetric criterion. In conjunction with these models the following viscoplasticity laws were considered. 1. Bingham's law. It has one parameter—viscosity. Multivariant numerical simulation of the creep of salt specimens allowed us to describe the growth rates of longitudinal deformation for all load levels considered. Analysis of the calibrated parameters showed that when the creep processes are described with the “non-associated Mohr-Coulomb + Bingham” model, the viscosity decreases with increase in the load level. Different results were obtained using the “associated Volumetric + Bingham” model—the viscosity varies in a certain range as the load level changes. 2. Duvaut-Lions law. It is notable for exploiting the solution of elastoplastic model. There is a single parameter involved— the relaxation time. Similar to Bingham's law, it is linear. Using the non-associated Mohr-Coulomb criterion, the solution of the systems of residuals is expressed in a closed form, which is an advantage for the numerical implementation. The results obtained by multivariant numerical simulation of the creep in salt specimens showed that for both yield criteria—the non- associated Mohr-Coulomb and the associated volumetric—there is a tendency for relaxation time to decrease with the load level of salt specimens to be increased. 3. Perzyna viscoplasticity law. It contains two parameters: the viscosity and the rate-sensitivity (strain rate hardening), which responsible for the material reaction to changes in strain rate. The main disadvantage of the law is the yield point being doubled when the rate-sensitivity in the limit is reduced to zero. Analysis of the calibrated parameters of the elastic- viscoplastic model by means of multivariant numerical simulation of creep has shown that for both plastic potentials considered, an increase in the load level leads to an increase in the viscosity parameter. In this case, the rate-sensitivity remains practically unchanged. 4. Perić's law. It eliminates the disadvantages of Perzyna's law for a similar set of input parameters. The analysis of the calibrated parameters for the elastic-viscoplastic model allowed us to conclude that when using the non-associated Mohr- Coulomb criterion, with increase in load level only the rate-sensitivity parameter (strain rate hardening) varies in a narrow range. The viscosity does not change. Simulation of the creep in salt specimens using the associated volumetric criterion as a plastic potential gives different results—the viscosity increases with the load level to be increased, while the rate-sensitivity does not change. Thus, all viscoplasticity laws considered in this paper reflect the growth of longitudinal deformation of salt specimens during creep. It should be noted that the yield criteria combined with Perić viscoplastic law are characterized by relatively simple, stable parametric support and acceptable description of steady-state creep. The use of viscoplastic deformation model allows us to point out the fracture regions formed during creep according to a specified strength criterion. ACKNOWLEDGEMENT he work was carried out as part of a major scientific project funded by the Ministry of Science and Higher Education of the Russian Federation (Agreement No. 075-15-2024-535 dated 23 April 2024). REFERENCES [1] Udartsev, A., Pankov, I. (2021). The study of silvinite and carnallite deformations under instaneous and sustained loading. Procedia Structural Integrity, 32(22). pp. 187-193. DOI: 10.1016/j.prostr.2021.09.027. [2] Li, Z., Kang, Y., Fan, J., Fourmeau, M., Jiang, D., Nelias, D. (2023). Creep–fatigue mechanical characteristics of salt rocks under triaxial loading: An experimental study. Engineering Geology, 322. DOI: 10.1016/j.enggeo.2023.107175. [3] Osipov, Yu.V., Voznesensky, A.S. (2022). Determination of rheological properties of bischofite from triaxial tests. Journal of Mining Science, 58(6). pp. 886-895. DOI: 10.1134/s1062739122060023. [4] Konstantinova, S.A., Chernopazov, S.A. (2007). Estimating salt rocks rheology by in-situ observations of vertical deformation of interchamber pillars. Journal of Mining Science. 43(3). pp. 17-23. DOI: 10.1007/s10913-007-0003-9. A T A. Baryakh et alii, Frattura ed Integrità Strutturale, 70(2024) 191-209; DOI: 10.3221/IGF-ESIS.70.11 209 [5] Baryakh, A., Tsayukov, A. (2022). Justification of fracture criteria for salt rocks. Frattura Ed Integrità Strutturale, 16(62), pp. 585–601. DOI: 10.3221/IGF-ESIS.62.40. [6] Labuz, J.F., Zang, A. (2012). Mohr–Coulomb Failure Criterion, Rock Mech Rock Eng, (45). pp. 975–979. DOI: 10.1007/s00603-012-0281-7. [7] Baryakh, A.A. and Samodelkina N.A. (2017). About one criteria of strength of rocks, Chebyshevskii Sbornik, 18(3), pp. 72-87. (In Russ.) DOI: 10.22405/2226-8383-2017-18-3-72-87. [8] Findley, W., Lai, J., Onaran, K. (1989). Creep and Relaxation of Nonlinear Viscoelastic Materials with an Introduction to Linear Viscoelasticity. New York, Dover Publications. ISBN: 0-486-66016-8. [9] Singh, A., Kumar, C., Kannan, L. G., Seshagiri Rao, K., Ayothiraman, R. (2018). Estimation of creep parameters of rock salt from uniaxial compression tests. Int. J. Rock Mech. Min. Sci, (107). pp. 243-248. DOI: 10.1016/j.ijrmms.2018.04.037. [10] Paraskevopoulou, C., Perras, M., Diederichs, M., Amann, F., Löw, S., Lam, T., Jensen, M. (2017). The three stages of stress relaxation - Observations for the time-dependent behaviour of brittle rocks based on laboratory testing. Engineering Geology, (216). pp. 56-75. DOI: 10.1016/j.enggeo.2016.11.010. [11] Morozov, I.A., Udarcev A.A., Pankov I.L. (2020). Laboratory deformation testing of salt rocks from the Gremyachinsk and Upper Kama deposits. Mining Inf. Anal. Bull, (10). pp. 16-28. (In Russ.) DOI: 10.25018/0236-1493-2020-10-0-16- 28. [12] de Souza Neto, E.A., Perić, D. and Owen, D.R.J. (2008). Computational Methods for Plasticity: Theory and Applications, Chichester, John Wiley & Sons Ltd. DOI: 10.1002/9780470694626. [13] Zienkiewicz, O.C., Taylor, R.L. and Fox, D. (2014). The Finite Element Method for Solid and Structural Mechanics (7th edition), Butterworth-Heinemann. DOI: 10.1016/C2009-0-26332-X. [14] de Borst, R., Crisfield, M.A., Remmers, J.J.C., Verhoosel, C.V. (2012). Non-Linear Finite Element Analysis of Solids and Structures, 2nd Edition, Chichester, John Wiley & Sons Ltd. DOI: 10.1002/9781118375938. [15] Quevedo, R., Firme, P., Roehl, D. (2019). Integration schemes with substepping algorithms for creep analysis in geomaterials. Int J Numer Anal Methods Geomech, 43, pp. 1467–1487. DOI: 10.1002/nag.2910. [16] Ottosen, N.S., Ristinmaa, M. (1996). Corners in plasticity—Koiter's theory revisited, IJSS, 33(25), pp. 3697-3721. DOI: 10.1016/0020-7683(95)00207-3. [17] ANSYS Mechanical APDL. Theory Reference. Release 2021 R2. << /ASCII85EncodePages false /AllowTransparency false /AutoPositionEPSFiles true /AutoRotatePages /None /Binding /Left /CalGrayProfile (Dot Gain 20%) /CalRGBProfile (sRGB IEC61966-2.1) /CalCMYKProfile (U.S. Web Coated \050SWOP\051 v2) /sRGBProfile (sRGB IEC61966-2.1) /CannotEmbedFontPolicy /Error /CompatibilityLevel 1.4 /CompressObjects /Tags /CompressPages true /ConvertImagesToIndexed true /PassThroughJPEGImages true /CreateJobTicket false /DefaultRenderingIntent /Default /DetectBlends true /DetectCurves 0.0000 /ColorConversionStrategy /CMYK /DoThumbnails false /EmbedAllFonts true /EmbedOpenType false /ParseICCProfilesInComments true /EmbedJobOptions true /DSCReportingLevel 0 /EmitDSCWarnings false /EndPage -1 /ImageMemory 1048576 /LockDistillerParams false /MaxSubsetPct 100 /Optimize true /OPM 1 /ParseDSCComments true /ParseDSCCommentsForDocInfo true /PreserveCopyPage true /PreserveDICMYKValues true /PreserveEPSInfo true /PreserveFlatness true /PreserveHalftoneInfo false /PreserveOPIComments true /PreserveOverprintSettings true /StartPage 1 /SubsetFonts true /TransferFunctionInfo /Apply /UCRandBGInfo /Preserve /UsePrologue false /ColorSettingsFile () /AlwaysEmbed [ true ] /NeverEmbed [ true ] /AntiAliasColorImages false /CropColorImages true /ColorImageMinResolution 300 /ColorImageMinResolutionPolicy /OK /DownsampleColorImages true /ColorImageDownsampleType /Bicubic /ColorImageResolution 300 /ColorImageDepth -1 /ColorImageMinDownsampleDepth 1 /ColorImageDownsampleThreshold 1.50000 /EncodeColorImages true /ColorImageFilter /DCTEncode /AutoFilterColorImages true /ColorImageAutoFilterStrategy /JPEG /ColorACSImageDict << /QFactor 0.15 /HSamples [1 1 1 1] /VSamples [1 1 1 1] >> /ColorImageDict << /QFactor 0.15 /HSamples [1 1 1 1] /VSamples [1 1 1 1] >> /JPEG2000ColorACSImageDict << /TileWidth 256 /TileHeight 256 /Quality 30 >> /JPEG2000ColorImageDict << /TileWidth 256 /TileHeight 256 /Quality 30 >> /AntiAliasGrayImages false /CropGrayImages true /GrayImageMinResolution 300 /GrayImageMinResolutionPolicy /OK /DownsampleGrayImages true /GrayImageDownsampleType /Bicubic /GrayImageResolution 300 /GrayImageDepth -1 /GrayImageMinDownsampleDepth 2 /GrayImageDownsampleThreshold 1.50000 /EncodeGrayImages true /GrayImageFilter /DCTEncode /AutoFilterGrayImages true /GrayImageAutoFilterStrategy /JPEG /GrayACSImageDict << /QFactor 0.15 /HSamples [1 1 1 1] /VSamples [1 1 1 1] >> /GrayImageDict << /QFactor 0.15 /HSamples [1 1 1 1] /VSamples [1 1 1 1] >> /JPEG2000GrayACSImageDict << /TileWidth 256 /TileHeight 256 /Quality 30 >> /JPEG2000GrayImageDict << /TileWidth 256 /TileHeight 256 /Quality 30 >> /AntiAliasMonoImages false /CropMonoImages true /MonoImageMinResolution 1200 /MonoImageMinResolutionPolicy /OK /DownsampleMonoImages true /MonoImageDownsampleType /Bicubic /MonoImageResolution 1200 /MonoImageDepth -1 /MonoImageDownsampleThreshold 1.50000 /EncodeMonoImages true /MonoImageFilter /CCITTFaxEncode /MonoImageDict << /K -1 >> /AllowPSXObjects false /CheckCompliance [ /None ] /PDFX1aCheck false /PDFX3Check false /PDFXCompliantPDFOnly false /PDFXNoTrimBoxError true /PDFXTrimBoxToMediaBoxOffset [ 0.00000 0.00000 0.00000 0.00000 ] /PDFXSetBleedBoxToMediaBox true /PDFXBleedBoxToTrimBoxOffset [ 0.00000 0.00000 0.00000 0.00000 ] /PDFXOutputIntentProfile () /PDFXOutputConditionIdentifier () /PDFXOutputCondition () /PDFXRegistryName () /PDFXTrapped /False /CreateJDFFile false /Description << /ARA /BGR /CHS /CHT /CZE /DAN /DEU /ESP /ETI /FRA /GRE /HEB /HRV (Za stvaranje Adobe PDF dokumenata najpogodnijih za visokokvalitetni ispis prije tiskanja koristite ove postavke. Stvoreni PDF dokumenti mogu se otvoriti Acrobat i Adobe Reader 5.0 i kasnijim verzijama.) /HUN /ITA /JPN /KOR /LTH /LVI /NLD (Gebruik deze instellingen om Adobe PDF-documenten te maken die zijn geoptimaliseerd voor prepress-afdrukken van hoge kwaliteit. De gemaakte PDF-documenten kunnen worden geopend met Acrobat en Adobe Reader 5.0 en hoger.) /NOR /POL /PTB /RUM /RUS /SKY /SLV /SUO /SVE /TUR /UKR /ENU (Use these settings to create Adobe PDF documents best suited for high-quality prepress printing. Created PDF documents can be opened with Acrobat and Adobe Reader 5.0 and later.) >> /Namespace [ (Adobe) (Common) (1.0) ] /OtherNamespaces [ << /AsReaderSpreads false /CropImagesToFrames true /ErrorControl /WarnAndContinue /FlattenerIgnoreSpreadOverrides false /IncludeGuidesGrids false /IncludeNonPrinting false /IncludeSlug false /Namespace [ (Adobe) (InDesign) (4.0) ] /OmitPlacedBitmaps false /OmitPlacedEPS false /OmitPlacedPDF false /SimulateOverprint /Legacy >> << /AddBleedMarks false /AddColorBars false /AddCropMarks false /AddPageInfo false /AddRegMarks false /ConvertColors /ConvertToCMYK /DestinationProfileName () /DestinationProfileSelector /DocumentCMYK /Downsample16BitImages true /FlattenerPreset << /PresetSelector /MediumResolution >> /FormElements false /GenerateStructure false /IncludeBookmarks false /IncludeHyperlinks false /IncludeInteractive false /IncludeLayers false /IncludeProfiles false /MultimediaHandling /UseObjectSettings /Namespace [ (Adobe) (CreativeSuite) (2.0) ] /PDFXOutputIntentProfileSelector /DocumentCMYK /PreserveEditing true /UntaggedCMYKHandling /LeaveUntagged /UntaggedRGBHandling /UseDocumentProfile /UseDocumentBleed false >> ] >> setdistillerparams << /HWResolution [2400 2400] /PageSize [612.000 792.000] >> setpagedevice