A viscoelastic, viscoplastic, and viscodamage constitutive model of salt rock for underground energy storage cavern

Computers and Geotechnics xxx (xxxx) xxxx

Contents lists available at ScienceDirect

Computers and Geotechnics journal homepage: www.elsevier.com/locate/compgeo

Research Paper

A viscoelastic, viscoplastic, and viscodamage constitutive model of salt rock for underground energy storage cavern ⁎

Jianqiang Denga,b, , Yaoru Liub, Qiang Yangb, Wei Cuia, Yinbang Zhua, Yi Liua, Bingqi Lia a b

State Key Laboratory of Simulation and Regulation of Water Cycle in River Basin, China Institute of Water Resources and Hydropower Research, Beijing 100038, China State Key Laboratory of Hydroscience and Engineering, Tsinghua University, Beijing 100084, China

A R T I C LE I N FO

A B S T R A C T

Keywords: Underground energy storage Salt rock Creep damage model Viscoelasticity Viscoplasticity Viscodamage

The underground energy storage in salt rock is of great strategic significance for ensuring the national energy safety and economic development. Due to the complex geological conditions, the construction of underground energy storage cavern in salt rock encounters great challenge and the long-term stability and safety issue of storage cavern in salt rock needs to be solved urgently. Based on the effective stress and the assumption of a small deformation, a unified constitutive model incorporating viscoelasticity, viscoplasticity, and viscodamage is developed by coupling the Kelvin-Voigt model, Duvaut-Lions model, and an improved Darabi viscodamage model to simulate the time-dependent mechanical behavior of salt rock for underground energy storage. The constitutive model is calibrated by experimental data to determine the parameters of viscoelasticity, viscoplasticity, and viscodamage. The model verification and discussion indicate that the constitutive model is able to effectively simulate the mechanical response of salt rock under different loading conditions and provide a safety guarantee for underground energy storage. The constitutive model is numerically implemented and successfully applied in long-term time-dependent deformation and failure analysis of underground energy storage cavern.

1. Introduction Salt rock has the features of low permeability, high ductility, microdamage healing, and water solubility [1–4] and has become an ideal medium for storing natural gas, petroleum, and other hydrocarbons, as well as radioactive waste [5–8]. Underground energy storage involves many key engineering and geological issues. The construction and operation of huge storage caverns creates great engineering disturbances in geological bodies of salt rock and dramatically changes the stress state and loading condition. A disturbed geological body of salt rock will undergo long-term time-dependent deformation and damage evolution that inevitably result in ground subsidence, cavern collapse, and even fire caused by energy leakage. These complex environmental problems and extreme natural disasters will threaten the safety of strategic energy reserves and the safety of people’s lives and property. Therefore, a study of the mechanical properties of salt rock under different stress states and loading conditions is crucial when undertaking a project of underground energy storage [9–12]. As described above, it was determined by the features of salt rock that the mechanical behavior of salt rock is time dependent and behaves as recoverable and unrecoverable (i.e., viscoelastic and viscoplastic) deformations. Numerous experimental investigations [13,14]

⁎

indicated that these mechanical responses generally behave as primary creep and steady-state creep under stress loading, acceleration creep under stress loading above the damage threshold, and creep recovery under stress unloading. Moreover, numerous studies have shown that the acceleration creep under stress loading and strain softening under constant strain rate loading are usually attributed to the damage evolution [15,16]. Wu et al. [17] indicated that actual stress is not constant and increases gradually due to damage accumulation over time during accelerating creep. The increased actual stress not only generates loading strain but also causes the steady creep rate to rise. This coupling possibly explains why salt rock shows nonlinear accelerating characteristics at the accelerating creep stage. In addition to physical simulation by experiments [7,18–20], numerical simulation by constitutive models is a main method for stability and failure analysis of underground energy storage cavern in salt rock. A large number of constitutive models have been proposed to simulate the viscoelastic and viscoplastic mechanical responses of salt rock. Cristescu [21] developed a general constitutive equation for salt rock that models transient creep and stationary creep. Jin and Cristescu [22] proposed an elastic and nonassociated viscoplastic model for transient creep that simulates creep, relaxation, dilatation, and contraction, as well as work hardening. Yahya et al. [23] developed an elaborate

Corresponding author. Tel.: +86 10 15001120925. E-mail address: [email protected] (J. Deng).

https://doi.org/10.1016/j.compgeo.2019.103288 Received 30 May 2019; Received in revised form 20 September 2019; Accepted 1 October 2019 0266-352X/ © 2019 Published by Elsevier Ltd.

Please cite this article as: Jianqiang Deng, et al., Computers and Geotechnics, https://doi.org/10.1016/j.compgeo.2019.103288

Computers and Geotechnics xxx (xxxx) xxxx

J. Deng, et al.

investigations of constitutive models for salt rock that couple viscodamage to both viscoelasticity and viscoplasticity in order to model the time and rate effects on damage evolution. The Kelvin-Voigt model is widely employed to simulate the viscoelastic mechanical response of soils [37]. As demonstrated by many researchers, the viscoelastic mechanical response of materials can be exactly simulated by the Kelvin-Voigt model [38,39]. The viscoplastic (rate-dependent) response of materials is widely modeled using overstress concept models such as the Perzyna model [40] and the Duvaut-Lions model [41]. The viscoplastic strain rate of the Perzyna model is directly determined by the value of the yield function. However, different from the Perzyna model, the viscoplastic strain rate of the Duvaut-Lions model is determined by the overstress, which is defined as the difference between the current stress state and the projected inviscid stress state onto the yield surface. For this reason, the Duvaut-Lions model has a remarkable advantage over the Perzyna model in that it can be easily adapted to multisurface plasticity and converges to the inviscid and elastic solution of rate-independent plasticity in limiting cases of the viscoplastic solution. However, the Perzyna model does not necessarily converge to the inviscid and elastic solution when it is applied to multisurface plasticity [42]. Pellet et al. [43] coupled viscoplasticity and viscodamage and developed a constitutive model to model the time-dependent behavior of rock, the damage evolution equation of which used the effective viscoplastic strain as the damage internal variable. Darabi et al. [44] developed a constitutive model for asphaltic material by coupling viscoelasticity, viscoplasticity, and viscodamage. Notably, the damage evolution law that Darabi et al. [44] proposed used the effective viscostrain including the viscoelastic and viscoplastic strain as the damage internal state variable. Therefore, the viscoelasticity and viscodamage were implicitly coupled. However, the Perzyna model accounting for viscoplasticity in the constitutive model developed by Darabi et al. [44] was not adapted to nonsmooth multisurface plasticity such as in the cone apex of the Drucker-Prager yield surface. Motivated by the above studies, the main purpose in our paper is to propose a unified constitutive model that incorporates viscoelasticity, viscoplasticity, and viscodamage by coupling the Kelvin-Voigt model, Duvaut-Lions model, and Darabi viscodamage model to simulate the time-dependent mechanical response and damage evolution of salt rock under a variety of loading conditions. In this unified constitutive model, the Kelvin-Voigt viscoelastic model is utilized to model the primary creep and recoverable viscoelastic strain, whereas the Duvaut-Lions overstress viscoplastic model is employed to model the steady-state creep and unrecoverable viscoplastic strain. Darabi’s viscodamage model is improved and coupled to the Kelvin-Voigt model and DuvautLions model to simulate the accelerated creep under stress loading and strain softening of the stress-strain relationship curve under a constant strain rate loading. The constitutive model is calibrated by experimental data of salt rock in order to determine the parameters of the viscoelastic, viscoplastic, and viscodamage model, and is then validated by comparing the experimental data with the model prediction. Finally, the constitutive model is numerically implemented successfully based on user material subroutine UMAT in Abaqus software and employed to analyze longterm time-dependent deformation and failure of underground energy storage cavern. The tensor notations and operations used throughout the paper are explained in detail in Appendix A.

constitutive model for rock salt to model the mechanical responses of unified plasticity, creep, and stress relaxation using an internal state variable (ISV) attached to isotropic and kinematic hardening. Heusermann et al. [24] described the nonlinear creep of rock salt including transient and steady-state creep behavior by the LUBBY2 constitutive model. Zhou et al. [15] applied a fractional derivative to develop a creep constitutive model by substituting an Abel dashpot for a Newtonian dashpot of the Nishihara model. Moghadam et al. [16] employed an elastoviscoplastic constitutive model to simulate transient and steady-state creep and investigate the dilatation, short-term failure, and long-term failure of rock salt. However, these models only model the initial transient and second steady-state creep, whereas the accelerated creep induced by damage evolution is not treated. The effect of damage evolution on the creep mechanical response has been investigated in several constitutive models for rock salt [1,25]. These models couple the creep and damage to model the time-dependent inelastic mechanical behavior and rupture of rock salt. In addition, Wang [26] introduced a damage accelerating limit into a creep model proposed by Carter and Hansen [27] to model the damage creep of salt rock. This model includes transient creep and steady-state creep modeled by Carter’s model, and an additional tertiary creep caused by damage evolution. Based on the definition of fractional derivative, Zhou et al. [28] developed a new damage mechanism-based creep model by substituting the variable-viscosity Abel dashpot for the Newtonian dashpot in the classical Nishihara model. Wu et al. [17] derived a new nonlinear creep damage constitutive model by assuming that the acceleration creep phase is a coupling process of loading and creeping, and then explained the mechanisms of nonlinear accelerating characteristics at accelerating creep stage. Müller et al. [29] used discrete element method to model the deformation and damage of rock salt by explicitly considering the fracture processes at grain scale. He et al. [30] carried out experimental investigation and damage modeling of salt rock subjected to fatigue loading. Li et al. [31] investigated the influence of temperature on the deformation, strength and cracking of rock salt by experiments and DEM simulations. Different approaches relevant to time-dependent deformation of rocks are also developed. Farhat et al. [32] proposed a micro-mechanics based viscoplastic model to describe time-dependent deformation for a class of clayey rocks. They extended the effective plastic yield criterion and used it as a viscoplastic loading function together with a suitable hardening law and a non-associated flow rule to develop the viscoplastic model. Based on a combined homogenization/thermodynamics framework, Zhao et al. [33] developed a unified micromechanics-based damage model for instantaneous and time-dependent inelastic behaviors of brittle rocks subjected to compressive stresses. The assumed that inelastic deformation is induced by frictional sliding along closed cracks, and strongly coupled with damage evolution by crack growth. Bikong et al. [34] considered mesoscopic and microscopic material scales and proposed a micro-macro model for the time-dependent behavior of clayey rocks. They assumed that creep deformation is induced by the time-dependent propagation of anisotropic microcracks inside the clay matrix. Zhou et al. [35] proposed a unified approach for modelling of elastic–plastic and viscoplastic behaviour coupled with induced damage in quasi-brittle rocks. They thought that rock materials may develop plastic deformation at different time scale due to evolution of material microstructure and material viscosity. Huang et al. [36] proposed an extension of Lahellec incremental variational approach for Callovo-Oxfordian argillites and considered the compressibility-dilatancy transition, influence of confining pressure and isotropic kinematic hardening effect of geomaterials. A great deal of effort has been devoted by the above researchers to developing constitutive models that couple creep and damage for salt rock or rock. However, the researchers did not propose a unified constitutive model that fully incorporates viscoelasticity, viscoplasticity, and viscodamage to model the time-dependent effects on the creep response and damage evolution. In fact, there have been few

2. Theoretical equations of constitutive model 2.1. Viscoelasticity constitutive equations The viscoelastic Kelvin-Voigt model is defined by the parallel arrangement of a spring element and a dashpot element [38]. Usually, a chain of several Kelvin-Voigt elements called the generalized Kelvin2

Computers and Geotechnics xxx (xxxx) xxxx

J. Deng, et al.

Voigt model is used to model more complicated creep behavior [39]. Moreover, the initial instantaneous elastic strain is also included in the generalized Kelvin-Voigt model. Then, the viscoelastic strain of the Kelvin-Voigt model under a three-dimensional stress state is written in integral form as

ε ve n

= D−1: σ +

∑ ⎡⎢ε0vei exp ⎛− ∫t i=1

+

∫t

t

D−1:

0

t

⎜

⎣

0

⎝

E σ 0 exp ⎜⎛− ηi ⎝

∫τ

Ei ⎞ dτ ⎟ ηi ⎠

Ei ⎞ ⎤ dτ ′⎟ dτ ηi ⎠ ⎥ ⎦

t

(1)

ε0vei

is the viscoelastic strain of the i-th Kelvin-Voigt element at where initial time t0 . τ and τ ′ are the integration variables. Ei denotes the viscoelastic modulus of the i-th spring element, and ηi denotes the viscosity coefficient of the i-th dashpot element. ε e = D−1: σ is the initial instantaneous elastic strain. σ denotes the stress tensor. D denotes the fourth-order elastic tensor, and E0 denotes the elastic modulus. For a constant stress σ , viscoelastic modulus Ei , and viscoelastic coefficient ηi , Eq. (1) is simplified to n

ε ve = D−1: σ +

E

⎧

∑ ⎨ε0vei exp ⎛− ηi Δt⎞ + ⎜

i=1

⎩

⎝

⎟

i

⎠

Fig. 1. Closest point projection method (CPPM) for nonassociated flow rule [43].

E0 −1 ⎡ E ⎫ D : σ ⎢1 − exp ⎛⎜− i Δt⎞⎟ ⎤ ⎥⎬ Ei η ⎝ i ⎠⎦⎭ ⎣ (2)

where Δt = t − t 0 . Setting φi = 1/ Ei , ri = Ei/ ηi , t0 = 0 , and ε0vei = 0 , we obtain n

ε ve = D−1: σ +

∑ φi E0 D−1: σ [1 − exp(−ri t )] i=1

(3)

Fig. 2. D-P viscoplasticity of nonassociated flow rule (CPPM) [43].

Correspondingly, the viscoelastic strain rate for multiaxial stress conditions is

function. In other words, the inviscid stress tensor σ̂ is unknown. Therefore, determining the inviscid stress tensor σ ̂ is crucial to the implementation of the Duvaut-Lions model. Generally, the inviscid stress tensor σ ̂ can be derived according to the CPPM, the generalized flow rule of plastic potential theory, and the consistency condition f (σ̂ ) = 0 . The viscoplastic strain rate expressed by the Duvaut-Lions model and the generalized flow rule is written in the following form:

n

ε ̇ve =

∑

(φi ri E0 D−1: σ − ri εive )

i=1

(4)

2.2. Viscoplasticity constitutive equations 2.2.1. Overstress model In the overstress theory of viscoplasticity, it is assumed that the stress σ only satisfies the equilibrium equation and is permitted to relax onto the yield surface with the time. The viscoplastic strain rate in the Duvaut-Lions model is defined as

ε ̇vp = Γ vpD−1: (σ − σ̂ )

ε ̇vp = Γ vpD−1: (σ − σ̂ ) = Γ vpωm

(6)

where the second-order direction tensor m = ∂g/ ∂σ governs the direction of the plastic flow. g = g (σ ) is a plasticity potential function. The positive scalar ω governs the magnitude of the plastic flow depending on the minimum distance between the stress state and the yield surface in a sense of norm. It is termed the associated flow rule or normality flow rule if the plasticity potential function is consistent with the yield function, i.e., g = f (σ ) , and the nonassociated flow rule if g ≠ f (σ ) . For the smooth and convex plasticity potential function surface g = g (σ ) , the direction tensor m = ∂g/ ∂σ always exists uniquely. When CPPM [45] is used, the direction tensor m can be calculated by σ̂ (implicit Euler method for CPPM) or by σ (explicit Euler method for CPPM). The implicit method is usually more robust than the explicit method. However, the implicit method usually requires more complicated calculation and is more difficult to implement.

(5)

Γ vp

= 1/ τr denotes a positive viscoplastic viscosity coefficient, τr where represents the viscoplasticity relaxation time, D denotes the fourthorder elastic tensor, and σ̂ = P(σ ) denotes the second-order inviscid stress tensor, which can be viewed as the closest point projection of σ onto the yield surface f (σ ) = 0 in the sense of the energy norm ∥σ − σ̂ ∥D−1 = (σ − σ̂ ): D−1: (σ − σ̂ ) . The closest point projection method (CPPM) for the nonassociated flow rule is shown in Fig. 1. It follows that σ̂ is unique for any σ provided that the yield and plastic potential function surfaces are convex. The execution of CPPM is illustrated by the vertical symbol in Fig. 2. However, it is worth noting that the stress σ̂ is not the geometrical closest point of σ in stand principle stress space but the extreme point stress to minimize the energy norm ∥σ − σ̂ ∥D−1. It is noteworthy that the viscoplastic strain rate of the Duvaut-Lions model is not directly calculated by the yield function value but is calculated by the overstress, which is in contrast to the Perzyna model. Therefore, the Duvaut-Lions model has a remarkable advantage in that it is able to easily account for multisurface plasticity. Unlike the Perzyna model, the Duvaut-Lions model does not directly give the explicit expression of the viscoplastic strain rate and the yield

2.2.2. Duvaut-Lions overstress model for Drucker–Prager viscoplasticity The Drucker–Prager (D-P) yield function and viscoplastic potential function similar to the D-P yield function are written as follows. (It is postulated that tensile stress is positive and compressive stress is negative.)

3

f = α I1 +

J2 − κ (p) = 0

(7)

g = βI1 +

J2

(8)

Computers and Geotechnics xxx (xxxx) xxxx

J. Deng, et al.

where I1 = tr(σ ) is the first stress tensor invariant, and J2 = tr(s2)/2 is the second deviatoric stress tensor invariant. s = σ − (I1/3) I denotes the deviatoric stress tensor, and I is the second-order unit tensor. α denotes a pressure-sensitive parameter associated with the internal friction angle, β denotes a parameter accounting for the dilation (β < α ) or contraction (β > α ), and κ (p) denotes the function of isotropic hardening associated with both the cohesive strength and the internal friction angle. This function depends on the effective viscoplastic strain p and can be expressed as [46]

κ (p) = κ 0 + κ1 [1 − exp (−κ2 p)]

∂g ∂g = ∂σ̂ ∂σ

This is an important property for the CPPM of the D-P yield function and nonassociated plasticity potential function. Therefore, the CPPM of the Duvaut-Lions model for the D-P yield function and nonassociated plasticity potential function is essentially the explicit integration algorithm that has the first-order accuracy of implicit integration. The CPPM is unconditionally convergent and stable for the convex yield surface and therefore possesses a high computational efficiency for numerical analyses of large load increments.

(9)

where κ 0 denotes the initial yield stress, κ1 denotes the upper-limit stress for the fully hardened material, and κ2 describes the transition rate between κ 0 and κ 0 + κ1. The material parameters α and κ 0 in this paper are given by

sinθ + 3 (3 + sinθ)

α=

3ccosθ + 3 (3 + sinθ)

κ0 =

sinθ 3 (3 − sinθ) 3ccosθ 3 (3 − sinθ)

2.3. Viscodamage constitutive equations A number of studies [47,48] indicated that the damage evolution under different loading conditions depends on the stress level, accumulated strain, and other factors such as the environmental temperature and humidity. Therefore, the general damage evolution equation is proposed as a function of the stress and viscous strain, in which the effect of stress is embodied in a power function of the so-called D-P type damage driving force, and the effect of the viscous strain is included in an exponential function of equivalent viscoplastic strain. When the temperature effect is neglected, the damage evolution equation is written as [44]

(10)

(11)

where θ and c refer to the internal friction angle and the cohesive strength, respectively. The D-P viscoplasticity of nonassociated flow rules (CPPM) is plotted in Fig. 2. At the cone apex of the D-P yield surface, I1∗ = k / α . The corresponding stress state is σ ∗ = [κ /(3α )] I . For the Drucker-Prager viscoplasticity, the following equations can be derived (for details, see Appendix B):

I1 ̂ = I1 − 9Kβω

J2 ̂ =

q

(14)

Eventually, the overstress of the D-P yield function is

σ − σ̂ =

〈f (σ ) 〉 ∂g 9αβK + G ∂σ

= σ=σ

〈f (σ ) 〉 ⎛ 1 ⎞ s⎟ ⎜3KβI + G 9αβK + G ⎝ J2 ⎠

(15)

If β = α , Eq. (15) represents the overstress for the associated flow rule. According to Eq. (13), the applicable condition of Eq. (15) is J2 ̂ ⩾ 0 , i.e.,

αI1 − γ J2 − κ⩽0

(16)

κ 1 κ I = ⎛I1 − ⎞ I + s 3α 3⎝ α⎠

q

⎜

〈f(σ ) 〉 ⎛ 1 ⎞ s⎟, (α I1 − γ J2 − κ⩽0) ⎜β I + 9αβ K + G ⎝ 2 J2 ⎠

1 κ 1 ⎤ ⎛I1 − ⎞ I + s , (α I1 − γ J2 − κ> 0) ε ̇vp = Γ vp ⎡ α⎠ 2G ⎦ ⎣ 9K ⎝

⎟

(23)

The first term in Eq. (23) is the convex and increasing function of the accumulative equivalent viscoplastic strain, whereas the second term in Eq. (23) is the convex and decreasing function of the accumulative equivalent viscoplastic strain.

(17)

Substituting Eqs. (15) and (17) into Eq. (5), respectively, yields

ε ̇vp = Γ vp

(22)

J2

Y ϕ ̇ = ⎛ ⎞ [Γ vd1 exp(k1 p) + Γ vd2 exp(−k2 p)] ⎝ Y0 ⎠

where γ = 9αβK / G . Notably, when αI1 − γ J2 − κ> 0 (the area denoted by the shadow in Fig. 2), the inviscid stress of CPPM is at the apex of the cone, and the overstress of the D-P yield function is

σ − σ̂ = σ −

Y = αI1 +

where denotes the viscosity coefficient of the viscodamage, q denotes a stress dependency parameter, and k is a material parameter. Y denotes the damage driving force calculated by the nominal stress, and Y0 denotes the reference damage driving force determined under a reference stress level. p is the accumulative equivalent viscoplastic strain t and can be calculated by p = ∫0 p ̇ dt . α is a constant of the D-P yield function, and I1 and J2 are invariants of σ . ϕ denotes the scalar damage variable ranging between ϕ = 0 (intact) and ϕ = 1 (eventual failure). When damage self-healing is not considered, the damage always increases with time or the accumulative equivalent viscoplastic strain. However, the damage rate may not always increase with time or the accumulative equivalent viscoplastic strain under different stress levels or at different damage evolution stages. It is noteworthy that the damage rate expressed by Eq. (21) is nonnegative and always increases with the accumulative equivalent viscoplastic strain. However, in some cases, this equation poorly fits the experimental data, as will be shown in Section 3. To characterize the damage evolution more accurately, the damage evolution equation is improved as

(13)

〈f (σ ) 〉 9αβK + G

(21)

⎟

Γ vd

Substituting Eq. (12) and Eq. (13) into the consistency condition f (I1,̂ J2)̂ = 0 , the nonnegative scalar ω is obtained:

ω=

Y ϕ ̇ = Γ vd ⎛ ⎞ exp(kp) ⎝ Y0 ⎠ ⎜

(12)

1 ⎞ J2 ⎜⎛1 − G ω⎟ J2 ⎠ ⎝

(20)

2.4. Coupling of viscodamage, viscoelastic, and viscoplastic models 2.4.1. Total strain decomposition For small deformation assumptions, the total strain under a constant nominal stress level can be split into elastic (instantaneous) and viscous (time-dependent) strains. Then, the viscous strain is split into viscoelastic and viscoplastic strains, i.e.,

(18)

(19)

It is indicated by Eq. (15) that if the flow rule of the smooth and convex plastic potential function is satisfied, the Duvaut-Lions model degenerates into the Perzyna model. Moreover, it can be proven that (for details, see Appendix C):

(24)

ε = εe + εv εv

=

ε ve

+

where ε , 4

(25)

ε vp

εe,

ε v,

ε ve

, and

ε vp

denote the total, elastic, viscous, viscoelastic,

Computers and Geotechnics xxx (xxxx) xxxx

J. Deng, et al.

and viscoplastic strain tensors, respectively.

vp ε¯ ̇ = Γ vp

2.4.2. Effective stress According to the effective stress concept in continuous damage mechanics (CDM), the relationship between the effective and nominal stress is written as [49]

σ σ¯ = 1−ϕ

σ (1 − ϕ)2

(26)

q

ε → ε¯ = ε

i=1

(35)

J¯2

3.1. Viscoelasticity and viscoplasticity model parameters Generally, the damage has not yet evolved or is insignificant and can be neglected under low stress levels and short-time loading. Therefore, the salt rock merely exhibits primary and steady-state creep. The parameters for the viscoelastic and viscoplastic model are calibrated by fitting this primary and steady-state creep curve. [For the fitting equations, see Eqs. (D.1)–(D.3) in Appendix D.] When the nonassociated viscoplasticity and isotropic hardening are not considered, β = α and κ (p) = κ 0 . The experimental data of salt rock and the fitted values of creep strain without damage under short-time uniaxial stress loading are illustrated in Fig. 3(a). In this paper, two Kelvin-Voigt elements are utilized to model the viscoelastic strain. The determined parameters for the viscoelastic and viscoplastic model are presented in Table 2.

(28)

3.2. Viscodamage model parameters When time is sufficient for damage evolution under a high stress level, accelerated creep in addition to the primary and steady-state creep occurs. Because the accelerated creep is caused by damage evolution, the parameters for the viscodamage model are determined by this primary, steady, and accelerated creep curve. The determined parameters for the viscoelastic and viscoplastic model without damage presented in Table 2 are utilized to predict the total creep including the primary, steady, and accelerated creep. However, the result predicted by the viscoelastic and viscoplastic model without coupling a viscodamage model greatly underestimates the total creep, including the accelerated creep induced by damage evolution. However, coupling the viscodamage model to the viscoelastic and viscoplastic model will compensate for this underestimated prediction.

(29) (30)

Table 1 Static mechanical parameters of salt rock.

n

∑ (φi ri E0 D−1: σ¯ − ri ε¯ive)

Y¯ = αI¯1 +

In this part, the procedure to calibrate the viscoelastic, viscoplastic, and viscodamage constitutive model is presented according to the experimental data of salt rock. The experimental data of salt rock originates from Özşen [13]. The basic mechanical parameters are presented in Table 1. For the convenience of plotting, the compression stress and strain are defined as positive, and the tensile stress and strain are defined as negative in the following figures.

The viscoelastic strain rate in the effective stress configuration of the material is ve ε¯ ̇ =

(34)

3. Calibration method of constitutive model

2.4.3. Viscoelastic, viscoplastic, and viscodamage models in effective stress space On the basis of the effective stress concept and strain equivalence hypothesis, in order to couple viscoelasticity, viscoplasticity, and viscodamage, all of the constitutive equations of the viscoelastic and viscoplastic models should be written using the effective stress σ¯ , which generates larger viscoelastic and viscoplastic deformation, rather than using the nominal stress tensor. Then, the coupling of the viscoelastic, viscoplastic, and viscodamage models is realized by substituting the effective stress tensor σ¯ for the nominal stress tensor σ :

σ (1 − ϕ)2

Y¯ (1 − ϕ)2 ⎤ vd1 ϕ̇ = ⎡ [Γ exp(k1 p¯) + Γ vd2 exp(−k2 p¯)] ⎢ ⎥ Y0 ⎣ ⎦

All of the variables with the superimposed notation “–” are calculated in the effective configuration.

(27)

There are two types of methods to combine viscoelasticity, viscoplasticity, and viscodamage. The first is to formulate constitutive equations in the effective stress space, in which the effective stress is the microscopic stress exerted on the intact configuration between microcracks and/or microvoids. The second is to express the constitutive equations in the nominal stress space, in which the nominal stress is the macroscopic stress exerted on the entire configuration including damaged and intact parts. Nevertheless, coupling the viscoelasticity, viscoplasticity, and viscodamage models by the effective stress appears to be more effective and stable in its numerical implementation [50,52].

σ → σ¯ =

(33)

As noted by Darabi, Al-Rub [44], to reflect the historical damage effects, the damage driving force Y is written using the nominal stress instead of the effective stress when the viscodamage model is coupled to the viscoelastic and viscoplastic models. Eq. (23) is then rewritten as

where ε and ε¯ denote the nominal and effective total strain tensors, respectively. ε e and ε¯e denote the nominal and effective elastic strain tensors, respectively. ε ve and ε¯ ve denote the nominal and effective viscoelastic strain tensors, respectively. ε vp and ε¯ vp denote the nominal and effective viscoplastic strain tensors, respectively. However, the strain equivalence hypothesis usually brings about a linear relationship between the stiffness modulus degradation and the damage evolution, which is not experimentally motivated. To remedy this issue and model the stiffness degradation with damage more accurately, Cicekli, Voyiadjis [50] and Konartakhteh [51] suggested that Eq. (26) be modified as

σ¯ =

(32)

1 ⎤ 1 ¯ κ vp ⎛I1 − ⎞ I + s¯ , (α ¯I1 − γ J¯2 − κ> 0) ε¯ ̇ = Γ vp ⎡ 2G 9K α ⎝ ⎠ ⎦ ⎣

where σ¯ and σ denote the effective and nominal stress tensors, respectively. It was indicated by Darabi, Al-Rub [44] that the difference between the nominal strain and effective strain is small and can be neglected when isotropic damage and small deformations are assumed. Therefore, the strain equivalence hypothesis is used here:

ε¯ = ε, ε¯e = ε e, ε¯ ve = ε ve, ε¯ vp = ε vp

〈f(σ¯ ) 〉 ⎛ 1 ⎞ βI + s¯⎟, (α ¯I1 − γ J¯2 − κ⩽0) 9αβ K + G ⎜ 2 J¯2 ⎠ ⎝

(31)

The viscoelastic strain rate in the effective stress configuration of the material is 5

E0 (MPa)

v

θ (˚)

c (MPa)

18,000

0.3

30

1

Computers and Geotechnics xxx (xxxx) xxxx

J. Deng, et al.

Fig. 4. Experimental data of salt rock by Özşen [13] and model predictions of creep strain without damage for low stress level under uniaxial stress loading. (The nonassociated viscoplasticity and isotropic hardening are not considered).

rock by Özşen, Özkan [13] and the fitted values of the creep strain including damage under long-time uniaxial stress loading are illustrated in Fig. 3(b). Finally, the determined parameters of the viscodamage model are presented in Table 3. The computational equations for the viscoelastic and viscoplastic strain and the total creep strain caused purely by damage evolution variable ϕ are shown as Eqs. (D.8)–(D.10) in Appendix D. These equations clearly indicate the interaction between damage evolution variable ϕ , the viscoelastic and viscoplastic strain components, and the acceleration process of the creep strain by damage evolution variable ϕ .

4. Verification and discussion of constitutive model In this part, the newly proposed viscoelastic, viscoplastic, and viscodamage constitutive model is verified by a series of comparisons between model predictions and experimental data. The experimental data of salt rock by Özşen, Özkan [13] and model predictions of the creep strain without damage for a low stress level under uniaxial stress loading are shown in Fig. 4. It is observed that salt rock only exhibits primary and steady creep under a low stress level. Only the primary and steady-state creep can be modeled without the viscodamage model. The model predictions for the primary and steady creep show strong agreement with the data of the experiment overall. The experimental data of salt rock by Özşen, Özkan [13] and the model predictions of creep strain including damage at different highstress levels under uniaxial stress loading are shown in Fig. 5. It is observed that salt rock exhibits primary and steady creep as well as accelerated creep under a high stress level. When the viscodamage model is coupled to the viscoelastic and viscoplastic model, the accelerated creep can be accurately modeled. The model predictions for the primary, steady, and accelerated creep and final creep rupture time show strong agreement with the data of the experiment overall. It is observed that an obvious accelerated creep starts when the time is earlier for a higher stress level. The corresponding damage is insignificant at the initial time under low stress levels, and increases with time. The higher the stress level, the faster the damage evolves. The experimental data of salt rock by Wang [26] and model predictions of strain without damage under multistage stress loading in which the confining pressure is zero are shown in Fig. 6. The experimental data of salt rock by Yang [53] and the model predictions of strain or creep strain without damage under multistage stress loading in which the confining pressure is nonzero are shown in Fig. 7. It is observed that the model predictions of the strain or creep strain show

Fig. 3. Experimental data of salt rock by Özşen [13] and fitted values of creep strain (a) without damage under short time uniaxial stress loading and (b) including damage under long time uniaxial stress loading. (The nonassociated viscoplasticity and isotropic hardening are not considered).

Table 2 Visco-elastic and visco-plastic model parameters of salt rock without damage. φ1 (MPa−1)

r 1 (h−1)

φ2 (MPa−1)

r2 (h−1)

Γvp (h−1)

2.2 × 10−4

4.4 × 10−3

3.61 × 10−4

4.08 × 10−1

3.92 × 10−2

Table 3 Visco-damage model parameters of salt rock. Γvd1 (h−1)

Y0 (MPa)

q

k1

Γvd2 (h−1)

k2

2.26 × 10−11

27

1

378

6.04 × 10−5

3

Hence, the stresses σ1, σ2 , and σ3 that are used to calculate the viscoelastic and viscoplastic strain are replaced by the effective stresses σ¯1, σ¯2 , and σ¯3 , respectively. Then, the damage variable ϕ at each time during the creep process can be calculated. By analyzing the damage curve against the accumulative equivalent viscoplastic strain, the parameters of the viscodamage model can be determined. The fitting equations are shown as Eqs. (D.4)–(D.7) in Appendix D. The experimental data of salt

6

Computers and Geotechnics xxx (xxxx) xxxx

J. Deng, et al.

Fig. 6. Experimental data of salt rock by Wang [22] and model predictions of strain without damage under multi-stage stress loading. (The nonassociated viscoplasticity and isotropic hardening are not considered).

Fig. 5. Experimental data of salt rock by Özşen [13] and model predictions of creep strain including damage at different high stress levels under uniaxial stress loading. (The nonassociated viscoplasticity and isotropic hardening are not considered).

good agreement with the experimental data overall under different loading steps. When the confining pressure remains constant, the creep strain rate increases as the axial stress increases under different loading steps. The viscoplastic strain rate is higher under a higher stress level or deviatoric stress. Therefore, the slope of the creep strain curve is steeper (Fig. 6). In fact, the increase in the creep strain has a great relationship with the deviatoric stress. The creep strain rates of two loading steps are nearly the same when the deviatoric stresses are nearly the same in Fig. 7(a). The creep strain rate tends to be constant when the deviatoric stresses are higher in Fig. 7(b). The experimental data of salt rock by Roberts, Buchholz [54] and model predictions of strain without damage under stress loading and unloading are shown in Fig. 8. It is observed that the model predictions of the strain growth during loading and the strain recovery during unloading show strong agreement with the experimental data overall under loading and unloading steps. The axial stress of the first loading step is 3 MPa, and the axial stress of the second loading step is 37 MPa. However, the lateral stresses of both loading steps are 20 MPa, which is between 3 MPa and 37 MPa. Therefore, the deviatoric stress of the first loading step is −17 MPa, and the deviatoric stress of the second loading

Fig. 7. Experimental data of salt rock by Yang [41] and model predictions of strain or creep strain without damage under multi-stage stress loading. (The nonassociated viscoplasticity and isotropic hardening are not considered). 7

Computers and Geotechnics xxx (xxxx) xxxx

J. Deng, et al.

Fig. 8. Experimental data of salt rock by Roberts [42] and model predictions of strain without damage under stress loading and unloading. (The nonassociated viscoplasticity and isotropic hardening are not considered).

step is 17 MPa. The loading directions of the loading steps are opposite to each other. Hence, the strain increases in the first loading step and decreases in the second loading step. The viscoplastic strain rate decreases during partial unloading or becomes zero during complete unloading, whereas the viscoelastic strain recovers partially during partial unloading or recovers to zero during complete unloading. The experimental data of salt rock by Wang [26] and model predictions of the creep strain including damage under triaxial stress loading are shown in Fig. 9. It is shown that the acceleration creep develops earlier when the lateral stress is relatively low and the deviatoric stress is relatively high. The experimental data of salt rock by Yang [53] and model predictions of the creep strain without damage under triaxial stress loading are shown in Fig. 10. There are three groups of experiments at different confining pressures. There are at least two test samples of different axial stress levels for each group of experiments. Therefore, this is a comprehensive set of experimental data. A more powerful and robust constitutive model of salt rock is needed to model this series of experiments. It is observed that the model predictions of the creep strain show good agreement with the experimental data overall under different axial loading stress levels and

Fig. 10. Experimental data of salt rock by Yang [41] and model predictions of creep strain without damage under triaxial stress loading. (The nonassociated viscoplasticity and isotropic hardening are not considered).

confining pressures. The experimental data of salt rock by Yang [53] and model predictions of the axial and lateral creep strain without damage under triaxial stress loading are shown in Fig. 11. There are some differences between the model predictions and experimental data of the lateral creep strain when the nonassociated viscoplasticity and isotropic hardening are not considered. However, the model predictions of both axial and lateral creep strains show good agreement with the experimental data overall as shown in Fig. 12. The cause may be the anisotropic properties of salt rock, which need further investigation. 5. Long-term time-dependent deformation and failure analysis of underground energy storage cavern The developed viscoelastic, viscoplastic, and viscodamage constitutive model is numerically implemented based on user material subroutine UMAT in Abaqus software and used to conduct long-term time-dependent deformation and failure analysis of underground energy storage cavern The height and diameter of ellipsoidal cavern are 300 m and 180 m, respectively. The burial depth of cavern center is 750 m. The domain of

Fig. 9. Experimental data of salt rock by Wang [22] and model predictions of creep strain including damage under triaxial stress loading. (The nonassociated viscoplasticity and isotropic hardening are not considered). 8

Computers and Geotechnics xxx (xxxx) xxxx

J. Deng, et al.

Fig. 11. Experimental data of salt rock by Yang [41] and model predictions of axial and lateral creep strain without damage under triaxial stress loading. (The nonassociated viscoplasticity and isotropic hardening are not considered).

Fig. 12. Experimental data of salt rock by Yang [41] and model predictions of axial and lateral creep strain without damage under triaxial stress loading. (The nonassociated viscoplasticity and isotropic hardening are considered).

finite element model is 450 × 450 × 1500 m. The parameters of salt rock are shown in Fig. 9. The long-term stability and failure of cavern at internal pressure 0 MPa and 2 MPa are analyzed. The distributions of maximum principle strain and damage of cavern at internal pressure 0 MPa and 2 MPa are shown in Figs. 13 and 14, respectively. When internal pressure is 2 MPa, the maximum values of maximum principle strain and damage at the time of 10 years is 2.5% and 0.031, respectively. It can be seen that the damage is very small. The damage zone is mainly located in the middle position of cavern. The damage at the edge of the cavern is greater than that in the surrounding rock of the cavern. However, when internal pressure is 0 MPa, the maximum values of maximum principle strain and damage at the time of about 1 year is 7.5% and nearly 1, respectively. It can be seen that the damage has almost reached its maximum value at the edge of the middle position of cavern. The large damage zone locally concentrates in the central position of cavern. The damage at the edge of the cavern is significantly greater than that in the surrounding rock of the cavern. This indicates that failure is beginning to occur in this area. The curves of maximum principle strain and damage versus time at four positions (as shown in Fig. 15) of cavern at internal pressure 0 MPa and 2MPaa are shown in Figs. 16 and 17, respectively.

It can be seen that the time-dependent deformation and damage at the edge of the middle position of cavern (point A) develop fastest. The next fastest are that at the edge of the lower position of cavern (point B). The third fastest are that at the edge of the upper position of cavern (point C). The least fast are that at the inner surrounding rock of the middle position of cavern (point D). The greater the internal pressure, the greater the difference between the four positions. The time-dependent deformation and damage tend to be constant and the cavern is stable when internal pressure is 2 MPa. However, the time-dependent deformation and damage develop rapidly and the cavern is instable when internal pressure is 0 MPa. This indicates that internal pressure has a great influence on the stability of cavern. The internal pressure should not be too low. It is necessary to maintain an appropriate internal pressure at the operation period of underground energy storage cavern. 6. Conclusions In this paper, a unified constitutive model of salt rock incorporating viscoelasticity, viscoplasticity, and viscodamage was first proposed based on the Kelvin-Voigt model, Duvaut-Lions model, and Darabi viscodamage model. Darabi’s viscodamage model was improved and 9

Computers and Geotechnics xxx (xxxx) xxxx

J. Deng, et al.

Fig. 13. Distribution of (a) maximum principle strain and (b) damage (time = 10 years, internal pressure = 2 MPa).

Fig. 14. Distribution of (a) maximum principle strain and (b) damage (time = 1.065 years, internal pressure = 0 MPa).

coupled to the Kelvin-Voigt model and the Duvaut-Lions model according to the effective stress concept of CDM. The overstress of the Duvaut-Lions model accounting for Drucker–Prager type viscoplasticity was derived. It is noteworthy that the Duvaut-Lions model used in this paper, which generally defines the viscoplastic strain rate by the overstress, has a remarkable advantage over the Perzyna model in that it is able to easily account for nonsmooth yield surfaces. The constitutive model was calibrated by experimental data on salt rock to determine the viscoelastic, viscoplastic, and viscodamage model parameters. The verification and discussion indicated that this newly developed constitutive model can reasonably model the primary creep, steady-state creep, accelerated creep in constant stress loading, and stress step-by-step loading and creep recovery in stress loading-unloading. Based on user material subroutine UMAT in Abaqus software, a three-dimensional finite element implementation of this constitutive model is realized. The constitutive model can effectively simulate mechanical responses including the long-term time-dependent deformation and damage evolution of geological bodies of salt rock and successfully conduct stability evaluation and failure analysis of cavern of large time span. This will provide a safety guarantee for underground energy storage. Fig. 15. Locations of elements A, B, C and D.

Declaration of Competing Interest No conflict of interest exits in the submission of this manuscript, and 10

Computers and Geotechnics xxx (xxxx) xxxx

J. Deng, et al.

Fig. 16. Curves of (a) maximum principle strain and (b) damage versus time of elements A, B, C and D (internal pressure = 2 MPa).

Fig. 17. Curves of (a) maximum principle strain and (b) damage versus time of elements A, B, C and D (internal pressure = 0 MPa).

manuscript is approved by all authors for publication. I would like to declare on behalf of my co-authors that the work described was original research that has not been published previously, and not under consideration for publication elsewhere, in whole or in part. All the authors listed have approved the manuscript that is enclosed.

National Natural Science Foundation for Young Scientists of China, China (Grant No. 51709281), Open Research Fund Program of the State Key Laboratory of Hydroscience and Engineering of Tsinghua University, China (Grant No. SS0203A062018), the National Natural Science Foundation of China, China (Grant No. 11572174 and 51779277), and Special Scientific Research Project of the China Institute of Water Resources and Hydropower Research, China (Grant No. SS0145B112018, SS0145B612017 and SS0145B392016).

Acknowledgements The authors acknowledge the financial support provided by the Appendix A

Throughout the paper, boldface italic letters are used to indicate tensors. The components of the tensors are written in a Cartesian coordinate system and have indices ranging from 1 to 3. Unless otherwise stated, the repeated Latin indices denote the Einstein summation convention. For the second-order tensors A and B , as well as the fourth-order tensor D , the operation of the tensors is defined as follows. The dot product of two second-order tensors is (A ·B )ij = Aik Bkj. The power of a second-order symmetrical tensor is (A2 )ij = Aik Akj. The double dot product of two second-order tensors is A : B = AijBji . The trace of a second-order tensor is tr(A) = Aii . The tensor product of two second-order tensors is (A ⊗ B )ijkl = Aij Bkl . The square product of two second-order tensors is (A⊠B)ijkl = Aik Bjl . The double dot product of a fourth-order tensor and a second-order tensor is (D : A)ij = Dijkl Alk . Appendix B Generally, the yield function f (σ ) and plasticity potential function g (σ ) can be written as f = f (I1, J2) and g = g (I1, J2) , respectively. The partial 11

Computers and Geotechnics xxx (xxxx) xxxx

J. Deng, et al.

derivatives of two stress invariants with respect to σ are ∂I1/ ∂σ = I and ∂J2 / ∂σ = s . Then, the partial derivative of the plasticity potential function g (σ ) with respect to σ can be derived:

∂g ∂J2 ∂g ∂I1 ∂g 1 = βI + + = s ∂J2 ∂σ ∂I1 ∂σ ∂σ 2 J2

(B.1)

The elastic tensor D can be written as (B.2)

D = λI ⊗ I + 2GI ⊠I

where λ = νE0/[(1 + ν )(1−2ν )] and G = E0/[2(1 + ν )] denote the Lame constants. E0 and ν denote the elastic modulus and Poisson’s ratio, respectively. The following equation can be obtained:

σ̂ = σ − ωD :

∂g ∂σ

(B.3)

Substituting Eq. (B.1) and Eq. (B.2) into Eq. (B.3) yields

1 σ̂ = σ − ω ⎡ (3λ + 2G ) βI + G ⎢ J2 ⎣

s⎤ ⎥ ⎦

(B.4)

For tr(I) = 3 and tr(s) = 0 , the first invariant of stress tensor σ̂ can be determined:

I1 ̂ = tr(σ̂ ) = tr(σ ) − ω [(3λ + 2G) β tr(I)] = I1 − 9Kβω

(B.5)

where K = E0/[3(1−2ν )] denotes the bulk modulus. The deviatoric component of stress tensor σ̂ can also be determined:

1 1 ⎞ s = ⎜⎛1 − ωG ⎟s J2 J2 ⎠ ⎝

s ̂ = s − ωG

(B.6)

Then, the following equation can be obtained:

1 ⎞ J2 ⎜⎛1 − G ω⎟ J2 ⎠ ⎝

J2 ̂ =

(B.7)

Appendix C For the D-P yield function f = αI1 + obtained:

J2 − κ (εevp) and nonassociated plasticity potential function g = βI1 +

1 ⎞ s ̂ = ⎜⎛1 − ωG ⎟s J2 ⎠ ⎝

J2 ̂ =

J2 , the following equations can be

(C.1)

J2 − ωG

(C.2)

∂g s = βI + ∂σ 2 J2

(C.3)

∂g ŝ = βI + ∂σ̂ 2 J2 ̂

(C.4)

According to Eq. (C.1) and Eq. (C.2), the following equation can be derived:

ŝ J2 ̂

=

(1 − ωG / J2 ) s J2 − ωG

=

s J2

(C.5)

Then, it follows that

∂g ∂g = ∂σ ∂σ̂

(C.6)

Appendix D The viscoelastic and viscoplastic strain and total creep strain without damage for a true triaxial stress state of σ1 ≠ 0 , σ2 ≠ 0 , and σ3 ≠ 0 can be calculated as follows: n

ε1ve =

∑ (σ1 − νσ2 − νσ3 ) φi [1 − exp(−ri t )]

(D.1)

i=1

12

Computers and Geotechnics xxx (xxxx) xxxx

J. Deng, et al.

ε1vp =

Γ vp 〈αI1 + 9αβK + G

J2 − κ 0 〉 ⎛⎜β + ⎝

2σ1 − σ2 − σ3 ⎞ ⎟t 6 J2 ⎠

(D.2)

ε1v = ε1ve + ε1vp

(D.3)

)2

)2

)2].

where I1 = σ1 + σ2 + σ3 and J2 = (1/6)[(σ1 − σ2 + (σ2 − σ3 + (σ3 − σ1 The viscoelastic and viscoplastic strain and total creep strain including damage for a true triaxial stress state of σ1 ≠ 0 , σ2 ≠ 0 , and σ3 ≠ 0 can be calculated as follows:

ε¯1ve =

ε¯1vp =

1 (1 − ϕ)2

n

∑ (σ1 − νσ2 − νσ3 ) φi [1 − exp(−ri t )]

(D.4)

i=1

Γ vp 9αβK + G

αI1 + J2 2σ − σ2 − σ3 ⎞ − κ 0 ⎜⎛β + 1 ⎟t (1 − ϕ)2 6 J2 ⎝ ⎠

(D.5)

q

αI1 + J2 ⎞ vd1 ϕ ̇ = ⎜⎛ exp(k1 p) + Γ vd2 exp(−k2 p)] ⎟ [Γ Y0 ⎝ ⎠

(D.6)

ε¯1v = ε¯1ve + ε¯1vp

(D.7)

where p =

t 0

∫

ε ̇vp : ε ̇vp 3β2 + 1 / 2

dt

The viscoelastic and viscoplastic strain and total creep strain caused purely by damage for a true triaxial stress state of σ1 ≠ 0 , σ2 ≠ 0 , and σ3 ≠ 0 can be calculated as follows:

1 ε1ved = ε¯1ve − ε1ve = ⎡ − 1⎤ 2 ⎢ ⎥ ⎣ (1 − ϕ) ⎦

n

∑ (σ1 − νσ2 − νσ3 ) φi [1 − exp(−ri t )]

(D.8)

i=1

Γ vp (αI1 + J2 ) ⎛ 1 2σ1 − σ2 − σ3 ⎞ − 1⎤ ε1vpd = ε¯1vp − ε1vp = ⎡ ⎜β + ⎟t 2 ⎢ (1 − ϕ) ⎥ 6 J2 ⎣ ⎦ 9αβK + G ⎝ ⎠

(D.9)

ε1vd = ε¯1v − ε1v = ε1ved + ε1vpd

(D.10)

References [19] [1] Chan K, Bodner S, Munson D. Recovery and healing of damage in WIPP salt. Int J Damage Mech. 1998;7(2):143–66. [2] Hunsche U, Hampel A. Rock salt—the mechanical properties of the host rock material for a radioactive waste repository. Eng Geol. 1999;52(3):271–91. [3] Stormont J. Conduct and interpretation of gas permeability measurements in rock salt. Int J Rock Mech Min. 1997;34(3):303. e1-e11. [4] Urai JL, Spiers CJ, Zwart HJ, Lister GS. Weakening of rock salt by water during longterm creep. 1986. [5] Ozarslan A. Large-scale hydrogen energy storage in salt caverns. Int J Hydrogen Energ. 2012;37(19):14265–77. [6] Deng J, Yang Q, Liu Y. Time-dependent behaviour and stability evaluation of gas storage caverns in salt rock based on deformation reinforcement theory. Tunn Undergr Space Technol 2014;42:277–92. [7] Deng J, Yang Q, Liu Y, Pan Y. Stability evaluation and failure analysis of rock salt gas storage caverns based on deformation reinforcement theory. Comput Geotech 2015;68:147–60. [8] Bracke G, Fischer-Appelt K. Methodological approach to a safety analysis of radioactive waste disposal in rock salt: An example. Prog Nucl Energ. 2015;84:79–88. [9] Wang T, Yang C, Chen J, Daemen JJK. Geomechanical investigation of roof failure of China's first gas storage salt cavern. Eng Geol 2018;243:59–69. [10] Wei L, Li Y, Yang C, Daemen JJK, Yun Y, Zhang G. Permeability characteristics of mudstone cap rock and interlayers in bedded salt formations and tightness assessment for underground gas storage cavern. Eng Geol 2015;193:212–23. [11] Singh A, Kumar C, Kannan LG, Rao KS, Ayothiraman R. Engineering properties of rock salt and simplified closed-form deformation solution for circular opening in rock salt under the true triaxial stress state. Engineering Geology. 2018;243. [12] Schulze O, Popp T, Kern H, Langer M, Talbot CJ. Development of damage and permeability in deforming rock salt. Eng Geol 2001;61(2–3):163–80. [13] Özşen H, Özkan İ, Şensöğüt C. Measurement and mathematical modelling of the creep behaviour of Tuzköy rock salt. Int J Rock Mech Min. 201466):128-35. [14] Yang C, Daemen J, Yin J-H. Experimental investigation of creep behavior of salt rock. Int J Rock Mech Min. 1999;36(2):233–42. [15] Zhou H, Wang C, Han B, Duan Z. A creep constitutive model for salt rock based on fractional derivatives. Int J Rock Mech Min. 2011;48(1):116–21. [16] Moghadam SN, Mirzabozorg H, Noorzad A. Modeling time-dependent behavior of gas caverns in rock salt considering creep, dilatancy and failure. Tunn Undergr Sp Technol. 2013;33:171–85. [17] Wu F, Chen J, Zou Q. A nonlinear creep damage model for salt rock. Int J Damage Mech 2018;105678951879264. [18] Liu W, Zhang Z, Chen J, Fan J, Jiang D, Jjk D, et al. Physical simulation of

[20] [21]

[22] [23] [24]

[25]

[26] [27] [28]

[29]

[30] [31]

[32] [33]

[34]

[35]

13

construction and control of two butted-well horizontal cavern energy storage using large molded rock salt specimens. Energy. 2019;185. Fan J, Jiang D, Liu W, Wu F, Chen J, Daemen J. Discontinuous fatigue of salt rock with low-stress intervals. Int J Rock Mech Min Sci 2019;115. Mansouri H, Ajalloeian R. Mechanical behavior of salt rock under uniaxial compression and creep tests. Int J Rock Mech Min Sci 2018;110. Cristescu N. A general constitutive equation for transient and stationary creep of rock salt. International journal of rock mechanics and mining sciences & geomechanics abstracts. Elsevier 1993:125–40. Jin J, Cristescu N. An elastic/viscoplastic model for transient creep of rock salt. Int J Plasticity. 1998;14(1):85–107. Yahya O, Aubertin M, Julien M. A unified representation of the plasticity, creep and relaxation behavior of rocksalt. Int J Rock Mech Min. 2000;37(5):787–800. Heusermann S, Rolfs O, Schmidt U. Nonlinear finite-element analysis of solution mined storage caverns in rock salt using the LUBBY2 constitutive model. Comput Struct. 2003;81(8):629–38. Cristescu N, Hunsche U. Determination of Nonassociated Constitutive Equation for Rock Salt from Experiments. Finite Inelastic Deformations—Theory and Applications: Springer, 1992. p. 511-23. Wang G. A new constitutive creep-damage model for salt rock and its characteristics. Int J Rock Mech Min. 2004;41(61-7). Carter NL, Hansen FD. Creep of rocksalt. Tectonophysics 1983;92(4):275–333. Zhou HW, Wang CP, Mishnaevsky L, Duan ZQ, Ding JY. A fractional derivative approach to full creep regions in salt rock. Mechanics of Time-Dependent Materials. 2013;17(3):413–25. Müller C, Frühwirt T, Haase D, Schlegel R, Konietzky H. Modeling deformation and damage of rock salt using the discrete element method. Int J Rock Mech Min Sci 2018;103:230–41. He M, Li N, Zhu C, Chen Y, Wu H. Experimental investigation and damage modeling of salt rock subjected to fatigue loading. Int J Rock Mech Min Sci 2019;114. Li W, Cheng Z, Yang C, Kang D, Hu W. Experimental and DEM investigations of temperature effect on pure and interbedded rock salt. J Nat Gas Sci Eng 2018;56. (S1875510018302208). Farhat F, Shen WQ, Shao JF. A micro-mechanics based viscoplastic model for clayey rocks. Comput Geotech 2017;89:92–102. Zhao LY, Zhu QZ, Xu WY, Dai F, Shao JF. A unified micromechanics-based damage model for instantaneous and time-dependent behaviors of brittle rocks. Int J Rock Mech Min Sci 2016;84:187–96. Bikong C, Hoxha D, Shao JF. A micro-macro model for time-dependent behavior of clayey rocks due to anisotropic propagation of microcracks. Int J Plast 2015;69:73–88. Zhou H, Jia Y, Shao JF. A unified elastic–plastic and viscoplastic damage model for quasi-brittle rocks. Int J Rock Mech Min Sci 2008;45(8).

Computers and Geotechnics xxx (xxxx) xxxx

J. Deng, et al.

[46] Lemaitre J, Chaboche J-L. Mechanics of solid materials. Cambridge University Press; 1994. [47] Shakiba M, Darabi MK, Al-Rub RKA, Masad EA, Little DN. Microstructural modeling of asphalt concrete using a coupled moisture–mechanical constitutive relationship. Int J Solids Struct. 2014;51(25):4260–79. [48] Al-Rub RKA, Darabi MK. A thermodynamic framework for constitutive modeling of time-and rate-dependent materials. Part I: Theory. Int J Plasticity. 2012;34(61-92. [49] Rabotnov IUrN. Creep problems in structural members: North-Holland Pub. Co., 1969. [50] Cicekli U, Voyiadjis GZ, Al-Rub RKA. A plasticity and anisotropic damage model for plain concrete. Int J Plasticity. 2007;23(10):1874–900. [51] Konartakhteh MD. Thermo-viscoelastic-viscoplastic-viscodamage-healing modeling of bituminous materials: Theory and computation [PhD dissertation]. USA: Texas A &M University; 2011. [52] Voyiadjis GZ, Al-Rub RKA, Palazotto AN. Thermodynamic framework for coupling of non-local viscoplasticity and non-local anisotropic viscodamage for dynamic localization problems using gradient theory. Int J Plasticity. 2004;20(6):981–1038. [53] Yang C. Time -dependent behavior of rock salt Experimental investigation and theoretical analysis [Doctoral dissertation]: University of Nevada, Reno, 2000. [54] Roberts LA, Buchholz SA, Mellegard KD, Düsterloh U. Cyclic Loading Effects on the Creep and Dilation of Salt Rock. Rock Mech Rock Eng 2015;48(6):2581–90.

[36] Huang Y, Guéry AC, Shao JF. Incremental variational approach for time dependent deformation in clayey rock. Int J Plast 2015;64:88–103. [37] Bratosin D, Sireteanu T. A nonlinear Kelvin-Voigt model for soils. Proceedings of the Romanian Academy. 2002;3(3):99-104. [38] Zienkiewicz O, Watson M, King I. A numerical method of visco-elastic stress analysis. Int J Mech Sci. 1968;10(10):807–27. [39] Argyris J, Doltsinis IS. Constitutive modelling and computation of non-linear viscoelastic solids. Part I: Rheological models and numerical integration techniques. Comput Method Appl M. 1991;88(2):135–63. [40] Perzyna P. Flmdamental Problems in Viseoplasticity. Adv Appl Mech. 1966;9(243. [41] Duvaut G, Lions JL. Inequalities in mechanics and physics. Springer; 1976. [42] Simo J, Kennedy J, Govindjee S. Non-smooth multisurface plasticity and viscoplasticity. Loading/unloading conditions and numerical algorithms. Int J Numer Meth Eng 1988;26(10):2161–85. [43] Pellet F, Hajdu A, Deleruyelle F, Besnus F. A viscoplastic model including anisotropic damage for the time dependent behaviour of rock. Int J Numer Anal Met. 2005;29(9):941–70. [44] Darabi MK, Al-Rub RKA, Masad EA, Huang C-W, Little DN. A thermo-viscoelastic–viscoplastic–viscodamage constitutive model for asphaltic materials. Int J Solids Struct. 2011;48(1):191–207. [45] Simo JC, Hughes TJ. Computational inelasticity. Springer Science & Business Media; 2006.

14