A Nonlinear Creep Model of Rock Salt and Its Numerical Implement in FLAC 3 D

Creep characteristics are integral mechanical properties of rock salt and are related to both long-term stability and security of rock salt repository. Rock salt creep properties are studied in this paper through employing combinedmethods of theoretical analysis and numerical simulation with a nonlinear creep model and the secondary development in FLAC software. A numerical simulation of multistage loading creep was developed with the model and resulting calculations were found consequently to coincide with previously tested data.


Preface
Rock salt creep research has been the subject of a large number of studies: King [1] produced a series of experimental research on rock salt creep properties under various temperatures and different axial stress levels, identifying the relationship between deformation and time.Paul [2] researched the creep properties of four types of rock salt.Hunsche [3,4] and Cristescu [5,6] employed diversified loading methods to conduct research on rock salt properties in multiple creep stages.Chan et al. [7,8] applied continuum damage mechanics in the creep analysis of rock salt to study the unelastic flow caused by damage from the transition of creep stage to the damage accumulation in accelerated creep stage and to the final destruction process.Based on the rheological model (Lubby2), using continuum damage mechanics and considering the displacement, hardening, softening and recovery mechanism, and the damage and damage recovery mechanism of rock salt, Lux and Hou [9] proposed the Hou/Lux creep damage model.Xiande et al. [10,11] researched Changshan and Qiaohu rock salt and observed that, following analysis of the evolution and destruction process of cracks on rock salt specimens, the two rock salt types differed relatively significantly in the creep damage process, mainly due to difference in NaCl content, grain size, and cementation between grains.Chen et al. [12] conducted creep property experiments of two rock salt types, rock salt from greisens salt mine and rock salt with interlayer.Creep constitutive relation was developed under different confining pressure and axial pressure for each rock salt type.
Rock salt creep modeling research has experienced enormous progress marked by numerous research achievements.Experimental conditions and technical limitations, however, have slowed the study of rock salt nonlinear creep modeling as related to the obvious nonlinear properties in the attenuation creep and steady-state creep rock salt stages, especially soft rock [13][14][15][16].On the basis of rock salt creep model research achievements and according to the indoor creep test results [17]

Rock Salt Nonlinear Creep Equation
According to research results [18][19][20], the viscosity coefficient in the creep model varies over time according to rules, which have the following three features: (1) the initial value is not 0; (2) there is monotone escalation with the increasing of time; (3) there exists upper limit value.
Assuming the viscosity coefficient of the Newton model varies with time according to the following formula: where  0 ,  refer to the material constant and  refers to creep time.
This study proposes the nonlinear creep equation for enhanced perspective of the nonlinear creep properties of rock salt (Figure 1).Thus, by analogical method, the axial creep equation of rock salt in three-dimensional stress status can be developed as follows.
When  0 <   , the nonlinear viscoelasticity plasticity creep model equation is When  0 ≥   , the nonlinear viscoelasticity plasticity creep model equation is (3)

The Realization of Nonlinear Creep Model
in FLAC 3D Software For the nonlinear Kelvin model, deviator stress   and deviator strain    have the following relation: where   refers to the shear modulus of the nonlinear Kelvin model and   (,) refers to the viscosity coefficient:   (,) =   0 (( + 1)/( + 1 + )).
For the nonlinear Maxwell model, deviator stress   and deviator strain rate ė   have the following relation: where   refers to the shear modulus of the nonlinear Maxwell model and   (,) refers to the viscosity coefficient:   (,) =   0 (( + 1)/( + 1 + )).For the nonlinear Bingham model, that is, where  refers to yield function, () refers to switch function, and   (,) refers to the viscosity coefficient of nonlinear Bingham model, and for secondary development with FLAC 3D software, formula (5) should first be written as increment form; that is, Applying central difference, formula (6) should be where Consolidating the formula (11), the updated deviator strain formula of Step i in the nonlinear Kelvin model can be where  = 1 +   Δ/2  (,) and  =   Δ/2  (,) − 1.Further, consolidating formula (12), we can acquire Similarly, the increment form of formula ( 7) can be expressed as For the nonlinear Bingham model, when  < 0, the increment form of formula ( 8) is When  ≥ 0, its increment form is Substitute formula (10) with formulas ( 13), (14), and (15a) and (15b) and consolidate.The updated stress formula of Step i in nonlinear viscoelasticity plasticity creep model is where when  < 0,  = 1/2 ; ,  mean the same as those appearing in formula (12).
Thus, stress-strain relations of the nonlinear viscoelasticity plasticity creep model can be expressed by formula (16) in the program.
Yield function is derived when the software automatically reads and recognizes the loaded stress level, introduced in this study as the concept of stress intensity: Known from the classic elastic-plastic mechanics, where   refers to deviator tensor of stress;   refers to stress tensor;   refers to spherical tensor of stress.Each component of the stress tensor can be indicated by the corresponding pointer with FLAC 3D software, and the stress intensity  can be calculated according to formula (18).After reading stress intensity with relevant yield function, software can automatically identify yield condition according to the definition of stress intensity; if there is uniaxial compression  =  1 and triaxial compression with equal confining pressure,  =  1 −  3 .As recursive function cannot be applied in FLAC 3D software, the nonstationary normalizing of model viscosity coefficient in this study is derived mainly by the accumulation of ps->Creep pointer (indicating the creep time increment Δ) internally installed in the FLAC 3D software.

Experiment Simulation. A triaxial compression creep
experiment is simulated with the program in this section to verify validity of the nonlinear viscoelasticity plasticity creep model computing program.The simulated test sample is a cylindrical standard sample with diameter of 50 mm and height of 100 mm and divided into 2560 units and 2827 points.Apply normal constraint to the bottom of the sample, vertical and uniformly distributed load to the top, and confining pressure to the circumferential surface area, respectively.Nonlinear viscoelasticity plasticity model is then applied as the creep model with experiment results [17] utilized as creep parameters and the loading scheme applied using fixed confining pressure and four-stage axial pressure (20 MPa, 25 MPa, 30 MPa, and 37 MPa), respectively.As the triaxial compression creep of rock salt did not enter into the accelerated creep stage in the experiment [17]      Accelerated creep calculation with the Cvisc model is not available; thus, a very large value is assumed as the critical load of accelerated creep for calculation in the study model to prevent the surrounding rock of the storage cavern from entering the accelerated creep stage.Deducting to a certain extent with the method provided in literature [21], we convert the elastic-plastic mechanics parameters obtained from the indoor compression experiment of rock salt samples into the parameter of the rock mass with the same creep parameters shown in Table 1.Refer to Table 2 for the specific parameter values.Calculations with the nonlinear model of this study utilize all parameters in Table 2 while calculations with the Cvisc model utilize all parameters except  in Table 2. Comparisons then reveal the influence of nonlinear coefficient  to the calculation results under identical parameters.

Analysis of Examples
(2) Comparison of the Calculation Results. Figure 5 displays the displacement contour line calculated by two models for the surrounding rock following two years of creep.Displacement distribution rules then, as calculated with the two models, are essentially consistent.Calculations by the study model are slightly larger, with maximum displacement of 1.27 m while maximum displacement calculated by the Cvisc model was 1.21 m.Displacement distribution range of the cavity walls varies as the red zone (parts with displacement above 1.1 m) of the study model nearly completely covers all cavity walls, while the red zone of the Cvisc model mainly concentrates at the top and shoulder of the cavity.
Figure 6 displays curves, calculated by the two models, for vertical displacement changes over time at the vault point of the cavity.Instantaneous deformation at the cavern excavation moment is nearly consistent as calculated by the two models.The two curves transform over time similarly through the attenuation deformation stage to the steady deformation stage.Compared with the Cvisc model, the attenuation deformation stage in the nonlinear model persists for a shorter time with larger deformation and higher deformation rate in the same period.The above indicates that when all the other parameters are the same, the nonlinear coefficient  in this model works to larger calculation results.
Figure 7 illustrates the distribution of plastic zones around the cavern as calculated by the two models.The distribution range of plastic zones calculated by the study model is slightly larger than the range calculated by the Cvisc model.Utilizing the Fish language programming, volumes of the plastic zones calculated by the study model were 13649 m 3 and 12542 m 3 as calculated by the Cvisc model, indicating that the nonlinear coefficient  in this model affects the calculation.

Conclusions
Secondary development of the nonlinear viscoelasticity plasticity creep model of rock salt is produced in this study and verifies the obtained model calculation program, drawing the following conclusions.
(1) Based on finite difference theory, this study deduces the finite difference expression form in detail for the nonlinear viscoelasticity plasticity creep model of rock salt in combination with the favorable secondary

Figure 2 :
Figure 2: The -direction (vertical) displacement nephogram for the test sample under different stress levels.

Figure 2
Figure2is the -direction (vertical) displacement contour line nephogram for the test sample under different stress levels.The calculated time for creep in the first three stages of loading is 15 h and is 4.5 h for the fourth stage of loading.Observations of the figure indicate that vertical displacement is the largest at the upper end and reduces with progression toward the bottom.Under the first three stages of loading, maximum vertical displacement of the test sample gradually increases with the increase of axial stress, at 0.696 mm, 0.975 mm, and 1.254 mm, respectively.Under the fourth stage of loading, however, as axial stress exceeds the critical load setting of accelerated creep, vertical displacement of the test sample increases rapidly.Figure3is the curve indicating axial pressure values for the relationship between time and the vertical displacement at the upper surface center of the test sample under different stress levels.Creep curve under the first three stages of loading alters with similar rules undergoing the two stages of attenuation creep and steady state creep.All instantaneous deformation, including creep deformation within the same

Figure 3
Figure2is the -direction (vertical) displacement contour line nephogram for the test sample under different stress levels.The calculated time for creep in the first three stages of loading is 15 h and is 4.5 h for the fourth stage of loading.Observations of the figure indicate that vertical displacement is the largest at the upper end and reduces with progression toward the bottom.Under the first three stages of loading, maximum vertical displacement of the test sample gradually increases with the increase of axial stress, at 0.696 mm, 0.975 mm, and 1.254 mm, respectively.Under the fourth stage of loading, however, as axial stress exceeds the critical load setting of accelerated creep, vertical displacement of the test sample increases rapidly.Figure3is the curve indicating axial pressure values for the relationship between time and the vertical displacement at the upper surface center of the test sample under different stress levels.Creep curve under the first three stages of loading alters with similar rules undergoing the two stages of attenuation creep and steady state creep.All instantaneous deformation, including creep deformation within the same

Figure 4 :
Figure 4: Rock salt gas storage model.

Figure 5 :
Figure 5: Contour lines for the displacement around the cavern calculated by the two models.

Figure 6 :( 2 )( 3 )
Figure 6: Curves for displacement at the vault point calculated by the two models.
, this study proposes and builds a nonlinear creep model of rock salt and develops a secondary corresponding nonlinear creep model with FLAC 3D software.Following numerical simulation of the creep test with the developed model, a comparison of simulated results and experimental results indicates that the developed creep model is workable and may even provide reference for creep calculation in engineering practice.
refer to the average deviator stress and deviator strain, respectively, of the nonlinear Kelvin model;    ,    refer to the deviator stress tensors of Step i and Step i-1, respectively;  ,  ,  ,  refer to the deviator strain tensors of Step i and Step i-1, respectively.

Table 1 :
Calculating parameter of creep model.
, critical load is unknown   , and the model parameter  30 and  in its accelerated creeping stage.Validity of the nonlinear viscoelasticity plasticity creep model computing program, with specific parameter values provided in Table 1, is established as   = 35 MPa,  30 = 30000 MPa⋅h, and  = 0.2.

Table 2 :
Calculation parameters of the rock salt gas storage example.