# A Nonlinear Creep Model of Rock Salt and Its Numerical Implement in FLAC3D

Authors: Xinrong Liu; Xin Yang; Junbao Wang
Journal: Advances in Materials Science and Engineering
 (2015)

---

## Abstract

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 combined methods of theoretical analysis and numerical simulation with a nonlinear creep model and the secondary development in FLAC3D 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.


---

## Body

## 1. 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–16]. On the basis of rock salt creep model research achievements and according to the indoor creep test results [17], this study proposes and builds a nonlinear creep model of rock salt and develops a secondary corresponding nonlinear creep model with FLAC3D 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.


## 2. Rock Salt Nonlinear Creep Equation

According to research results [18–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:(1)
η
t
=
η
0
t
+
1
t
+
1
+
m
,where η
0, m refer to the material constant and t refers to creep time.This study proposes the nonlinear creep equation for enhanced perspective of the nonlinear creep properties of rock salt (Figure1).Figure 1
Nonlinear viscoelasticity plasticity creep model.Thus, by analogical method, the axial creep equation of rock salt in three-dimensional stress status can be developed as follows.Whenσ
0
<
σ
s, the nonlinear viscoelasticity plasticity creep model equation is(2)
ε
=
σ
0
E
1
+
σ
0
η
10
t
+
m
ln
⁡
t
+
1
+
σ
0
E
2
1
-
exp
⁡
-
E
2
η
20
t
+
m
ln
⁡
t
+
1
.Whenσ
0
≥
σ
s, the nonlinear viscoelasticity plasticity creep model equation is(3)
ε
=
σ
0
E
1
+
σ
0
η
10
t
+
m
ln
⁡
t
+
1
+
σ
0
E
2
1
-
exp
⁡
-
E
2
η
20
t
+
m
ln
⁡
t
+
1
+
σ
0
-
σ
s
η
30
t
1
-
n
t
2
.


## 3. The Realization of Nonlinear Creep Model in FLAC3D Software

### 3.1. The Differential Form of Nonlinear Creep Model

As each part of the model is in series connection, there is equal stress and added strain; that is,(4)
S
i
j
=
S
i
j
M
=
S
i
j
K
=
S
i
j
B
(5)
e
i
j
=
e
i
j
M
+
e
i
j
K
+
e
i
j
B
.For the nonlinear Kelvin model, deviator stressS
i
j and deviator strain e
i
j
K have the following relation:(6)
S
i
j
=
2
η
(
m
,
t
)
K
e
˙
i
j
K
+
2
G
K
e
i
j
K
,where G
K refers to the shear modulus of the nonlinear Kelvin model and η
(
m
,
t
)
K refers to the viscosity coefficient: η
(
m
,
t
)
K
=
η
0
K
t
+
1
/
t
+
1
+
m.For the nonlinear Maxwell model, deviator stressS
i
j and deviator strain rate e
˙
i
j
M have the following relation: (7)
e
˙
i
j
M
=
S
˙
i
j
2
G
M
+
S
i
j
2
η
(
m
,
t
)
M
,where G
M refers to the shear modulus of the nonlinear Maxwell model and η
(
m
,
t
)
M refers to the viscosity coefficient: η
(
m
,
t
)
M
=
η
0
M
t
+
1
/
t
+
1
+
m.For the nonlinear Bingham model, that is,(8)
e
˙
i
j
B
=
ϕ
F
2
η
(
n
,
t
)
B
,where F refers to yield function, ϕ
F refers to switch function, F
=
σ
1
-
σ
3
-
σ
s,(9)
ϕ
F
=
0
F
<
0
F
F
≥
0
,
η
(
n
,
t
)
B
=
η
0
B
1
-
n
t
3
/
1
+
n
t
,and η
(
n
,
t
)
B refers to the viscosity coefficient of nonlinear Bingham model, and for secondary development with FLAC3D software, formula (5) should first be written as increment form; that is,(10)
Δ
e
i
j
=
Δ
e
i
j
M
+
Δ
e
i
j
K
+
Δ
e
i
j
B
.Applying central difference, formula (6) should be(11)
S
-
i
j
Δ
t
=
2
η
(
m
,
t
)
K
Δ
e
i
j
K
+
2
G
K
e
-
i
j
K
Δ
t
,where Δ
e
i
j
K
=
e
i
j
K
,
N
-
e
i
j
K
,
O, S
-
i
j
=
S
i
j
N
+
S
i
j
O
/
2, and e
-
i
j
K
=
e
i
j
K
,
N
+
e
i
j
K
,
O
/
2, among which S
-
i
j, e
-
i
j
K refer to the average deviator stress and deviator strain, respectively, of the nonlinear Kelvin model; S
i
j
N, S
i
j
O refer to the deviator stress tensors of Step i and Step i-1, respectively; e
i
j
K
,
N, e
i
j
K
,
O refer to the deviator strain tensors of Step i and Step i-1, respectively.Consolidating the formula (11), the updated deviator strain formula of Step i in the nonlinear Kelvin model can be(12)
e
i
j
K
,
N
=
1
A
S
i
j
N
+
S
i
j
O
Δ
t
4
η
(
m
,
t
)
K
-
B
e
i
j
K
,
O
,where A
=
1
+
G
K
Δ
t
/
2
η
(
m
,
t
)
K and B
=
G
K
Δ
t
/
2
η
(
m
,
t
)
K
-
1.Further, consolidating formula (12), we can acquire(13)
Δ
e
i
j
K
=
e
i
j
K
,
N
-
e
i
j
K
,
O
=
1
A
S
i
j
N
+
S
i
j
O
Δ
t
4
η
m
,
t
K
-
A
+
B
e
i
j
K
,
O
.Similarly, the increment form of formula (7) can be expressed as(14)
Δ
e
i
j
M
=
Δ
S
i
j
2
G
M
+
S
-
i
j
2
η
(
m
,
t
)
M
Δ
t
.For the nonlinear Bingham model, whenF
<
0, the increment form of formula (8) is(15a)
Δ
e
i
j
B
=
0
.When F
≥
0, its increment form is(15b)
Δ
e
i
j
B
=
S
-
i
j
B
-
2
/
3
σ
s
2
η
(
n
,
t
)
B
Δ
t
.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(16)
S
i
j
N
=
1
a
Δ
e
i
j
+
B
A
+
1
e
i
j
K
,
O
+
b
S
i
j
O
F
<
0
,
1
a
Δ
e
i
j
+
B
A
+
1
e
i
j
K
,
O
+
b
S
i
j
O
+
σ
s
Δ
t
3
η
(
n
,
t
)
B
F
≥
0
,where when F
<
0, a
=
1
/
2
G
M
+
Δ
t
/
4
η
(
m
,
t
)
M
+
Δ
t
/
4
A
η
(
m
,
t
)
K and b
=
1
/
2
G
M
-
Δ
t
/
4
η
(
m
,
t
)
M
-
Δ
t
/
4
A
η
(
m
,
t
)
K; when F
≥
0
,
a
=
1
/
2
G
M
+
Δ
t
/
4
η
(
m
,
t
)
M
+
Δ
t
/
4
A
η
(
m
,
t
)
K
+
Δ
t
/
4
η
(
m
,
t
)
B and b
=
1
/
2
G
M
-
Δ
t
/
4
η
(
m
,
t
)
M
-
Δ
t
/
4
A
η
(
m
,
t
)
K
-
Δ
t
/
4
η
(
m
,
t
)
B; A, B 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:(17)
q
=
σ
i
=
1
2
σ
1
-
σ
2
2
+
σ
2
-
σ
3
2
+
σ
3
-
σ
1
2
0.5
.Known from the classic elastic-plastic mechanics,(18)
q
=
3
2
S
i
j
S
i
j
S
i
j
=
σ
i
j
-
σ
m
δ
i
j
σ
m
=
1
3
σ
11
+
σ
22
+
σ
33
,where S
i
j refers to deviator tensor of stress; σ
i
j refers to stress tensor; σ
m refers to spherical tensor of stress.Each component of the stress tensor can be indicated by the corresponding pointer with FLAC3D software, and the stress intensity q 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 q
=
σ
1 and triaxial compression with equal confining pressure, q
=
σ
1
-
σ
3. As recursive function cannot be applied in FLAC3D 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 Δ
t) internally installed in the FLAC3D software.


### 3.2. Verification of the Model Program

#### 3.2.1. 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], critical load is unknown σ
s, and the model parameter η
30 and n 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 σ
s
=
35 MPa, η
30
=
30000 MPa·h, and n
=
0.2.Table 1
Calculating parameter of creep model.
K/MPa
G
1/MPa
G
2/MPa
η
10/MPa·h
η
2
0/MPa·h
m
2409
1112
3553
356650
50678
40
η
30/MPa·h
n
c/MPa
φ/°
σ
t/MPa
σ
s/MPa
300000
0.2
5.8
40
1
35Figure2 is the Y-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.Figure 2
TheY-direction (vertical) displacement nephogram for the test sample under different stress levels.
(a)
Axial pressure 20 MPa
(b)
Axial pressure 25 MPa
(c)
Axial pressure 30 MPa
(d)
Axial pressure 37 MPaFigure3 is 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 time and the creep rate of the test sample, with the increase of triaxial compression, increases gradually. Simulated creep experiment in the program increases gradually aligning with the indoor experiment result. Under the fourth stage of loading, as the axial stress exceeds the critical load of accelerated creep, the creep curve experiences three stages, attenuation creep, steady state creep, and accelerated creep. The test sample is vulnerable to damage upon entering the accelerated creep stage resulting in the nonlinear accelerated increase of deformation over time. This also agrees with the actual situation, proving the rationality and validity of the secondary development program for the nonlinear viscoelasticity plasticity creep model in this paper.Figure 3
Creep curve for test sample under different stress levels.


#### 3.2.2. Analysis of Examples

(1) Calculation Model and Parameter. Viability of this model calculation program for practical engineering is demonstrated with a simple rock salt gas storage model and by comparing the calculation results with the Cvisc model installed in FLAC3D. A portion of the model (1/4) was symmetrically set up to conserve calculation time for calculating the creep during a period of 2 years. Gas storage occurred spherically with a radius of 20 m and the origin of coordinates placed in the center of the sphere. Horizontal, upper, and lower boundaries of the model were at 200 m resulting in a rectangular calculation area of 200 m × 200 m × 400 m as displayed in Figure 4. Normal constraint was applied to the bottom and four vertical surfaces of the model and 20 MPa was applied to equally distribute load to the top surface of the model.Figure 4
Rock salt gas storage model.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 m in Table 2. Comparisons then reveal the influence of nonlinear coefficient m to the calculation results under identical parameters.Table 2
Calculation parameters of the rock salt gas storage example.
K/MPa
G
1/MPa
G
2/MPa
η
10/MPa·h
η
20/MPa·h
c/MPa
φ/°
σ
t/MPa
m
1927
890
3553
356650
50678
1.74
34
1
40(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 5
Contour lines for the displacement around the cavern calculated by the two models.
(a)
The nonlinear model of this paper
(b)
Cvisc modelFigure6 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 m in this model works to larger calculation results.Figure 6
Curves for displacement at the vault point calculated by the two models.Figure7 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 m in this model affects the calculation.Figure 7
Distribution of plastic zones around the cavern calculated by the two models.
(a)
The nonlinear model of this paper
(b)
Cvisc model


## 3.1. The Differential Form of Nonlinear Creep Model

As each part of the model is in series connection, there is equal stress and added strain; that is,(4)
S
i
j
=
S
i
j
M
=
S
i
j
K
=
S
i
j
B
(5)
e
i
j
=
e
i
j
M
+
e
i
j
K
+
e
i
j
B
.For the nonlinear Kelvin model, deviator stressS
i
j and deviator strain e
i
j
K have the following relation:(6)
S
i
j
=
2
η
(
m
,
t
)
K
e
˙
i
j
K
+
2
G
K
e
i
j
K
,where G
K refers to the shear modulus of the nonlinear Kelvin model and η
(
m
,
t
)
K refers to the viscosity coefficient: η
(
m
,
t
)
K
=
η
0
K
t
+
1
/
t
+
1
+
m.For the nonlinear Maxwell model, deviator stressS
i
j and deviator strain rate e
˙
i
j
M have the following relation: (7)
e
˙
i
j
M
=
S
˙
i
j
2
G
M
+
S
i
j
2
η
(
m
,
t
)
M
,where G
M refers to the shear modulus of the nonlinear Maxwell model and η
(
m
,
t
)
M refers to the viscosity coefficient: η
(
m
,
t
)
M
=
η
0
M
t
+
1
/
t
+
1
+
m.For the nonlinear Bingham model, that is,(8)
e
˙
i
j
B
=
ϕ
F
2
η
(
n
,
t
)
B
,where F refers to yield function, ϕ
F refers to switch function, F
=
σ
1
-
σ
3
-
σ
s,(9)
ϕ
F
=
0
F
<
0
F
F
≥
0
,
η
(
n
,
t
)
B
=
η
0
B
1
-
n
t
3
/
1
+
n
t
,and η
(
n
,
t
)
B refers to the viscosity coefficient of nonlinear Bingham model, and for secondary development with FLAC3D software, formula (5) should first be written as increment form; that is,(10)
Δ
e
i
j
=
Δ
e
i
j
M
+
Δ
e
i
j
K
+
Δ
e
i
j
B
.Applying central difference, formula (6) should be(11)
S
-
i
j
Δ
t
=
2
η
(
m
,
t
)
K
Δ
e
i
j
K
+
2
G
K
e
-
i
j
K
Δ
t
,where Δ
e
i
j
K
=
e
i
j
K
,
N
-
e
i
j
K
,
O, S
-
i
j
=
S
i
j
N
+
S
i
j
O
/
2, and e
-
i
j
K
=
e
i
j
K
,
N
+
e
i
j
K
,
O
/
2, among which S
-
i
j, e
-
i
j
K refer to the average deviator stress and deviator strain, respectively, of the nonlinear Kelvin model; S
i
j
N, S
i
j
O refer to the deviator stress tensors of Step i and Step i-1, respectively; e
i
j
K
,
N, e
i
j
K
,
O refer to the deviator strain tensors of Step i and Step i-1, respectively.Consolidating the formula (11), the updated deviator strain formula of Step i in the nonlinear Kelvin model can be(12)
e
i
j
K
,
N
=
1
A
S
i
j
N
+
S
i
j
O
Δ
t
4
η
(
m
,
t
)
K
-
B
e
i
j
K
,
O
,where A
=
1
+
G
K
Δ
t
/
2
η
(
m
,
t
)
K and B
=
G
K
Δ
t
/
2
η
(
m
,
t
)
K
-
1.Further, consolidating formula (12), we can acquire(13)
Δ
e
i
j
K
=
e
i
j
K
,
N
-
e
i
j
K
,
O
=
1
A
S
i
j
N
+
S
i
j
O
Δ
t
4
η
m
,
t
K
-
A
+
B
e
i
j
K
,
O
.Similarly, the increment form of formula (7) can be expressed as(14)
Δ
e
i
j
M
=
Δ
S
i
j
2
G
M
+
S
-
i
j
2
η
(
m
,
t
)
M
Δ
t
.For the nonlinear Bingham model, whenF
<
0, the increment form of formula (8) is(15a)
Δ
e
i
j
B
=
0
.When F
≥
0, its increment form is(15b)
Δ
e
i
j
B
=
S
-
i
j
B
-
2
/
3
σ
s
2
η
(
n
,
t
)
B
Δ
t
.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(16)
S
i
j
N
=
1
a
Δ
e
i
j
+
B
A
+
1
e
i
j
K
,
O
+
b
S
i
j
O
F
<
0
,
1
a
Δ
e
i
j
+
B
A
+
1
e
i
j
K
,
O
+
b
S
i
j
O
+
σ
s
Δ
t
3
η
(
n
,
t
)
B
F
≥
0
,where when F
<
0, a
=
1
/
2
G
M
+
Δ
t
/
4
η
(
m
,
t
)
M
+
Δ
t
/
4
A
η
(
m
,
t
)
K and b
=
1
/
2
G
M
-
Δ
t
/
4
η
(
m
,
t
)
M
-
Δ
t
/
4
A
η
(
m
,
t
)
K; when F
≥
0
,
a
=
1
/
2
G
M
+
Δ
t
/
4
η
(
m
,
t
)
M
+
Δ
t
/
4
A
η
(
m
,
t
)
K
+
Δ
t
/
4
η
(
m
,
t
)
B and b
=
1
/
2
G
M
-
Δ
t
/
4
η
(
m
,
t
)
M
-
Δ
t
/
4
A
η
(
m
,
t
)
K
-
Δ
t
/
4
η
(
m
,
t
)
B; A, B 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:(17)
q
=
σ
i
=
1
2
σ
1
-
σ
2
2
+
σ
2
-
σ
3
2
+
σ
3
-
σ
1
2
0.5
.Known from the classic elastic-plastic mechanics,(18)
q
=
3
2
S
i
j
S
i
j
S
i
j
=
σ
i
j
-
σ
m
δ
i
j
σ
m
=
1
3
σ
11
+
σ
22
+
σ
33
,where S
i
j refers to deviator tensor of stress; σ
i
j refers to stress tensor; σ
m refers to spherical tensor of stress.Each component of the stress tensor can be indicated by the corresponding pointer with FLAC3D software, and the stress intensity q 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 q
=
σ
1 and triaxial compression with equal confining pressure, q
=
σ
1
-
σ
3. As recursive function cannot be applied in FLAC3D 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 Δ
t) internally installed in the FLAC3D software.


## 3.2. Verification of the Model Program

### 3.2.1. 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], critical load is unknown σ
s, and the model parameter η
30 and n 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 σ
s
=
35 MPa, η
30
=
30000 MPa·h, and n
=
0.2.Table 1
Calculating parameter of creep model.
K/MPa
G
1/MPa
G
2/MPa
η
10/MPa·h
η
2
0/MPa·h
m
2409
1112
3553
356650
50678
40
η
30/MPa·h
n
c/MPa
φ/°
σ
t/MPa
σ
s/MPa
300000
0.2
5.8
40
1
35Figure2 is the Y-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.Figure 2
TheY-direction (vertical) displacement nephogram for the test sample under different stress levels.
(a)
Axial pressure 20 MPa
(b)
Axial pressure 25 MPa
(c)
Axial pressure 30 MPa
(d)
Axial pressure 37 MPaFigure3 is 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 time and the creep rate of the test sample, with the increase of triaxial compression, increases gradually. Simulated creep experiment in the program increases gradually aligning with the indoor experiment result. Under the fourth stage of loading, as the axial stress exceeds the critical load of accelerated creep, the creep curve experiences three stages, attenuation creep, steady state creep, and accelerated creep. The test sample is vulnerable to damage upon entering the accelerated creep stage resulting in the nonlinear accelerated increase of deformation over time. This also agrees with the actual situation, proving the rationality and validity of the secondary development program for the nonlinear viscoelasticity plasticity creep model in this paper.Figure 3
Creep curve for test sample under different stress levels.


### 3.2.2. Analysis of Examples

(1) Calculation Model and Parameter. Viability of this model calculation program for practical engineering is demonstrated with a simple rock salt gas storage model and by comparing the calculation results with the Cvisc model installed in FLAC3D. A portion of the model (1/4) was symmetrically set up to conserve calculation time for calculating the creep during a period of 2 years. Gas storage occurred spherically with a radius of 20 m and the origin of coordinates placed in the center of the sphere. Horizontal, upper, and lower boundaries of the model were at 200 m resulting in a rectangular calculation area of 200 m × 200 m × 400 m as displayed in Figure 4. Normal constraint was applied to the bottom and four vertical surfaces of the model and 20 MPa was applied to equally distribute load to the top surface of the model.Figure 4
Rock salt gas storage model.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 m in Table 2. Comparisons then reveal the influence of nonlinear coefficient m to the calculation results under identical parameters.Table 2
Calculation parameters of the rock salt gas storage example.
K/MPa
G
1/MPa
G
2/MPa
η
10/MPa·h
η
20/MPa·h
c/MPa
φ/°
σ
t/MPa
m
1927
890
3553
356650
50678
1.74
34
1
40(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 5
Contour lines for the displacement around the cavern calculated by the two models.
(a)
The nonlinear model of this paper
(b)
Cvisc modelFigure6 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 m in this model works to larger calculation results.Figure 6
Curves for displacement at the vault point calculated by the two models.Figure7 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 m in this model affects the calculation.Figure 7
Distribution of plastic zones around the cavern calculated by the two models.
(a)
The nonlinear model of this paper
(b)
Cvisc model


## 3.2.1. 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], critical load is unknown σ
s, and the model parameter η
30 and n 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 σ
s
=
35 MPa, η
30
=
30000 MPa·h, and n
=
0.2.Table 1
Calculating parameter of creep model.
K/MPa
G
1/MPa
G
2/MPa
η
10/MPa·h
η
2
0/MPa·h
m
2409
1112
3553
356650
50678
40
η
30/MPa·h
n
c/MPa
φ/°
σ
t/MPa
σ
s/MPa
300000
0.2
5.8
40
1
35Figure2 is the Y-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.Figure 2
TheY-direction (vertical) displacement nephogram for the test sample under different stress levels.
(a)
Axial pressure 20 MPa
(b)
Axial pressure 25 MPa
(c)
Axial pressure 30 MPa
(d)
Axial pressure 37 MPaFigure3 is 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 time and the creep rate of the test sample, with the increase of triaxial compression, increases gradually. Simulated creep experiment in the program increases gradually aligning with the indoor experiment result. Under the fourth stage of loading, as the axial stress exceeds the critical load of accelerated creep, the creep curve experiences three stages, attenuation creep, steady state creep, and accelerated creep. The test sample is vulnerable to damage upon entering the accelerated creep stage resulting in the nonlinear accelerated increase of deformation over time. This also agrees with the actual situation, proving the rationality and validity of the secondary development program for the nonlinear viscoelasticity plasticity creep model in this paper.Figure 3
Creep curve for test sample under different stress levels.


## 3.2.2. Analysis of Examples

(1) Calculation Model and Parameter. Viability of this model calculation program for practical engineering is demonstrated with a simple rock salt gas storage model and by comparing the calculation results with the Cvisc model installed in FLAC3D. A portion of the model (1/4) was symmetrically set up to conserve calculation time for calculating the creep during a period of 2 years. Gas storage occurred spherically with a radius of 20 m and the origin of coordinates placed in the center of the sphere. Horizontal, upper, and lower boundaries of the model were at 200 m resulting in a rectangular calculation area of 200 m × 200 m × 400 m as displayed in Figure 4. Normal constraint was applied to the bottom and four vertical surfaces of the model and 20 MPa was applied to equally distribute load to the top surface of the model.Figure 4
Rock salt gas storage model.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 m in Table 2. Comparisons then reveal the influence of nonlinear coefficient m to the calculation results under identical parameters.Table 2
Calculation parameters of the rock salt gas storage example.
K/MPa
G
1/MPa
G
2/MPa
η
10/MPa·h
η
20/MPa·h
c/MPa
φ/°
σ
t/MPa
m
1927
890
3553
356650
50678
1.74
34
1
40(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 5
Contour lines for the displacement around the cavern calculated by the two models.
(a)
The nonlinear model of this paper
(b)
Cvisc modelFigure6 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 m in this model works to larger calculation results.Figure 6
Curves for displacement at the vault point calculated by the two models.Figure7 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 m in this model affects the calculation.Figure 7
Distribution of plastic zones around the cavern calculated by the two models.
(a)
The nonlinear model of this paper
(b)
Cvisc model


## 4. 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 development platform of FLAC3D software. The study realizes the secondary development of creep model with vc++ program, obtaining the dynamic linking calculation program of this model.
(2)
Results of the triaxial compression creep experiment simulation indicate that when axial stress is lower than critical load of accelerated creep, the test sample’s creep curve undergoes two stages, attenuation creep and steady state creep; instantaneous strain, creep strain, and creep rate of the samples gradually increase with the increase of stress level; when axial stress attains critical load of accelerated creep, the creep curve undergoes three stages: attenuation creep, steady state creep, and accelerated creep. Simulation results align with the actual situation, proving the validity and rationality of the secondary development calculation program for the study model.
(3)
Comparing engineering examples simulated by the study model in this paper to samples simulated by the Cvisc model installed in FLAC3D software, both displacement and plastic zones calculated by the study model were larger than those calculated by the Cvisc model.



---

*Source: 285158-2015-03-25.xml*