Parametric assessment of salt cavern performance using a creep model describing dilatancy and failure

International Journal of Rock Mechanics & Mining Sciences 79 (2015) 250–267

Contents lists available at ScienceDirect

International Journal of Rock Mechanics & Mining Sciences journal homepage: www.elsevier.com/locate/ijrmms

Parametric assessment of salt cavern performance using a creep model describing dilatancy and failure Saeed Nazary Moghadam a,b,n, Kimia Nazokkar c, Richard J. Chalaturnyk a,b,d, Hasan Mirzabozorg e a Department of Civil and Environmental Engineering, Geotechnical Engineering Center, University of Alberta, GeoREF/GeoCERF, 6-380 Donadeo Innovation Centre for Engineering, Edmonton, Alberta, Canada T6H 1G9 b Reservoir Geomechanics Research Group, University of Alberta, Edmonton, Alberta, Canada c Department of Civil, Geological and Mining Engineering, École Polytechnique de Montréal, Montreal, Quebec, Canada d Foundation CMG Research Chair in Reservoir Geomechanics for Unconventional Resources, University of Alberta, Rm 6-222 Donadeo Innovation Centre for Engineering, Edmonton, Alberta, Canada T6G 1H9 e Faculty of Civil Engineering, Department of Structural Engineering, K.N. Toosi University of Technology, Tehran, Iran

art ic l e i nf o

a b s t r a c t

Article history: Received 31 January 2013 Received in revised form 19 June 2015 Accepted 30 June 2015

The present paper discusses a systematic approach to study the time-dependent performance of underground storage caverns in rock salt. An elasto-viscoplastic constitutive model capable of describing dilatancy, short-term failure and long-term failure during creep of rock salt was implemented in a Lagrangian finite element formulation. Using the finite element formulation, a series of numerical analyses was performed to investigate the influence of the cavern geometry and non-salt interbeds on the performance of underground storage caverns excavated in rock salt in the framework of large displacements. The attained results demonstrated the capability of the proposed approach in performance assessment of the underground caverns. & 2015 Elsevier Ltd. All rights reserved.

Keywords: Nonlinear finite element Creep Rock mechanics Rock salt Underground cavern Oil and gas storage

1. Introduction As a consequence of the world's growing energy demands, underground storage caverns are being used for economical and safe storage of petroleum products. In this regard, rock salt is considered to be an ideal material for underground storage due to its low permeability, healing capacity and availability. Another benefit of the utilization of salt formations as media for underground storages is that caverns can be created in these formations by solution mining techniques, which are more economical than other conventional excavation techniques. As a result, underground storage caverns excavated in rock salt are playing a crucial role in providing storage capacity for petroleum products. n Corresponding author at: Department of Civil and Environmental Engineering, Geotechnical Engineering Center, University of Alberta, GeoREF/GeoCERF, 6-380 Donadeo Innovation Centre for Engineering, Edmonton, Alberta, Canada. T6H 1G9. Tel.: þ 1 780 492 3953; fax: þ1 780 492 0249. E-mail addresses: [email protected], [email protected] (S. Nazary Moghadam), [email protected] (K. Nazokkar), [email protected] (R.J. Chalaturnyk), [email protected] (H. Mirzabozorg).

http://dx.doi.org/10.1016/j.ijrmms.2015.06.012 1365-1609/& 2015 Elsevier Ltd. All rights reserved.

The stability of underground caverns excavated in rock salt is of vital importance for safe storage. The specific factors influencing the cavern stability are cavern internal pressure, cavern geometry, cavern depth and the mechanical behavior of the rock formation around the underground cavern. In addition, in the case of heterogeneous salt formations that are interspersed with non-salt sedimentary interbeds, the mechanical behavior of the interbeds has a dramatic effect on the stability of the caverns. Experimental evidence reveals that the crucial mechanical property of rock salt is the inelastic time-dependent or creep behavior that involves inelastic volumetric changes. The consequence of this behavior is that inelastic ground movement around an underground opening can progress for a long period of time, which brings about significant closure of the opening. Also, due to the viscoplastic behavior, stress relaxation can occur around the opening, which may lead to the instability of the cavern. Hence, in order to assure a safe storage, the influence of the above-mentioned factors should be taken into account in the stability analyses of the underground caverns. So far, numerous viscoplastic constitutive models have been developed to describe the time-dependent behavior of rock salt. In this regard, Herrmann et al.1,2 proposed an empirical creep model

S. Nazary Moghadam et al. / International Journal of Rock Mechanics & Mining Sciences 79 (2015) 250–267

of rock salt for the Waste Isolation Pilot Plant known as WIPP creep model. The WIPP creep model is based on the existing creep models for metals, which neglect the inelastic volumetric changes. Aubertin et al.3 and Yahya et al.4 formulated a viscoplastic constitutive model in which a single set of equations and material parameters were used to represent the short-term and long-term ductile behavior of rock salt. In that research, it was assumed that the stress state remains far enough from the failure condition, and therefore the inelastic volumetric straining was not considered. Furthermore, in this regard, works of Shao et al.,5 Zhou et al.6 and others7,8 are worth mentioning. Moreover, the development of efficient computational models for the time-dependent analysis of underground storage caverns has been one of the most important research activities.9–20 In particular, Dawson and Munson9 and Munson et al.10 carried out research on the time-dependent deformations around an opening excavated in rock salt. In that research, a constitutive model of rock salt neglecting inelastic volumetric changes was developed and implemented into a finite element procedure neglecting large deformation effects to numerically simulate the creep deformation of rock salt around an opening. Heusermann et al.11 presented a numerical procedure using finite element method and a creep model of rock salt neglecting inelastic volumetric changes to simulate the time-dependent behavior of underground caverns in rock salt and established a framework in which a method for the assessment of the stability and usability of the caverns is provided. Mostly, creep models of rock salt have been developed based on existing viscoplastic models of metals. In these models, first, uniaxial constitutive equations are defined, and after wards they are generalized for triaxial stress state based on classical theory of plasticity. One of the drawbacks of these models is that they do not take into account the effects of confining stress, represented by the first invariant of stress tensor, on the creep behavior. Therefore, these models do not take into account the relaxation of confining stress during viscoplastic flow. According to these models, since only deviatoric stress relaxes during viscoplastic flow, the failure of rock salt as a function time has not been taken into account properly. In contrast, triaxial experimental measurements reveal that not only deviatoric stress, but also confining stress relaxes by time.21–24 This means that unlike the prediction of conventional creep models, the stress state in rock salt can approach the failure criteria after a period of time. Another issue is that the conventional creep models cannot describe the inelastic volumetric changes including inelastic compressibility and dilatancy due to the opening and closing of microcracks during creep. Nevertheless, this behavior plays a vitally important role in long-term safety assessment of repositories in rock salt in terms of the seepage of stored fluid and cavern stability, particularly in case of the presence of non-salt interbeds. Therefore, these models cannot realistically reflect semi-brittle behavior of rock salt, microfracturing mechanism as well as failure mechanism in short-term and longterm. It should also be mentioned that large deformation effects are neglected in the majority of numerical models developed to analyze the behavior of underground caverns in rock salt. However, in long-term, the dimensions of deep underground openings in rock salt dramatically change which leads to large deformation of material elements near the opening surface. So far, several research works have been done to tackle the mentioned issues.21–28 In some of these investigations, the aim is to incorporate damage mechanics into the conventional viscoplastic models to account for semi-brittle creep behavior. In this regard, investigations conducted by Chan et al.27,28 are worth mentioning. Another remedy for the mentioned issues is to directly develop triaxial viscoplastic models based on triaxial experimental data. These models are capable of describing inelastic volumetric changes, ductile and semi-brittle creep behavior and failure in short-term and long-

251

term.21–26 In this regard, Cristescu21 and Cristescu and Hunsche22– proposed a viscoplastic constitutive model for rock salt to describe transient and steady state creep, inelastic volumetric changes, short-term failure and damage in the framework of small deformations. Jin and Cristescu25 derived a transient creep model for rock salt which was more appropriate for numerical analyses. Paraschiv-Munteanu and Cristescu26 formulated a semi-analytical procedure using finite element discretization to analyze the variation with time of the stress state in rock salt around a plane strain borehole during transient and steady state creep considering inelastic volumetric changes. This kind of constitutive model comprehensively describes rock salt behavior up to failure and provides an effective tool for the performance assessment of repositories in rock salt for long-term. However, so far, this model has not been implemented into a general numerical formulation, e.g. based on finite element or finite difference method, with the capability of modeling problems with different boundary and initial conditions to systematically conduct performance assessment of repositories in rock salt for different conditions. Another issue is that large deformation effects are neglected in this model. In this regard, Wang et al.16 utilized the WIPP creep model neglecting semi-brittle behavior, inelastic dilatancy and large deformation effects to analyze time-dependent behavior of repositories in rock salt with the presence of non-salt interbeds using finite difference method. Then, the criteria developed by Cristescu21 for inelastic volumetric behavior of rock salt during viscoplastic flow were employed by the authors to assess the obtained numerical results to discuss about the caverns performance in terms of seepage and stability. However, the creep model utilized in that research could not evaluate inelastic volumetric changes and relaxation of confining stress during viscoplastic flow, which forms the basis of the utilized criteria. Also, in that research, since a comprehensive failure criterion was not defined, very large displacements are predicted without indicating any failure in rock salt. In the present paper, efforts have been made to formulate a numerical model capable of resolving the mentioned issues. Also, systematic numerical analyses were performed to show the capabilities of the model in design and performance assessment of repositories in rock salt in long-term. In the present paper, the performance of underground caverns in rock salt as a function of different factors is investigated. To this end, an elasto-viscoplastic constitutive model is employed to represent creep, dilatancy, short-term failure as well as long-term failure behavior of rock salt. The constitutive model is employed by a Lagrangian finite element algorithm to simulate the timedependent response of underground openings excavated in rock salt in the framework of large displacements. A series of finite element analyses is carried out to investigate the influence of different parameters, namely side wall slope, cavern floor and roof shape, height-to-diameter ratio as well as non-salt interbed number and location, on the performance of the underground caverns. In these analyses, the mentioned parameters are varied to determine their effects on the cavern performance.

24

2. Elasto-viscoplastic constitutive model Based on experimental evidence, it can be postulated that the creep deformation of rock salt is comprised of three phases, namely transient (initial), steady state (secondary or stationary) and accelerated (tertiary) creep.29 In the transient phase, creep deformation continuously decreases in time up to a stable state, while, in the steady state phase, creep deformation progresses with a constant rate under constant or nearly constant stress state. Also, if rock salt reaches the creep failure condition, tertiary or accelerated creep phase is encountered where creep deformation

252

S. Nazary Moghadam et al. / International Journal of Rock Mechanics & Mining Sciences 79 (2015) 250–267

with an increasing rate is occurred. The elasto-viscoplastic model of geomaterials proposed by Cristescu21 and Cristescu and Hunsche22–24 was formulated based on a great number of experimental data obtained by performing true triaxial tests on cubic samples. In these experiments, first, the samples are loaded hydrostatically, i.e. mean stress increases while deviatoric stress is equal to zero. In the next stage, i.e. deviatoric phase, the mean stress is kept constant while the deviatoric stress is increased. Based on these experimental data, it was observed that the stress space can be divided into two separate domains, namely compressibility and dilatancy domains where geomaterials exhibit two distinct time-dependent behaviors during creep, i.e. compressible ductile and dilatant semi-brittle behavior. Therefore, the concept of compressibility–dilatancy (C–D) boundary was introduced in the stress space inside which geomaterials reflect compressible behavior, while, above this boundary, geomaterials display dilatant behavior during creep. The definition of the C–D boundary plays a crucial role in this constitutive model, particularly in describing damage during creep of rock salt. In microcracked media such as geomatrials, the mechanical behavior significantly relies on the opening and closure status of the preexisting defects such as microcracks. At high confining pressures, i.e. in compressibility region where mean stress is high relative to deviatoric stress, it is assumed that the microcracks are closed. Therefore, due to high normal stresses existing on crack surfaces, the frictional resistance is so high that the shear-induced propagation of microcracks is very unlikely. Therefore, as long as the stress state is located in the compressibility region, geomaterials can inelastically contract during creep due to closing of pre-existing microcracks which results in the decrease of damage, i.e. healing, in geomaterials. In contrast, when the stress state is located in dilatancy domain, low confining stress is encountered, i.e. mean stress is low relative to deviatoric stress. As a result, microcracks can slide by shear during creep, which leads to the development of wing-tips in some of the sliding microcracks. The wing-tip cracks can be opened during creep leading to inelastic volumetric expansion. Therefore, it can be inferred that the development of damage in geomaterials involves the generation and growth of microcracks leading to inelastic dilatancy21–28. From the above discussion, it can be concluded that the time-dependent damage in geomatrials commences only when the deviatoric stress is so high relative to mean stress that the stress state is located in the dilatancy region. Therefore, the concept of C–D boundary plays a vitally important role in defining the long-term failure criteria during creep, i.e. creep failure. In this regard, the ductile creep behavior occurs if the stress state is located in the compressibility domain. In contrast, as the stress state goes beyond the dilatancy threshold, defined by C–D boundary, so that it is located in dilatancy region, semi-brittle creep behavior is encountered. In these circumstances, the process of generation and growth of microcracks causes the progressive degradation of material properties, which can lead to long-term failure after a period of time. The criteria defined by the C–D boundary by itself cannot represent damage evolution to determine when long-term failure occurs. Therefore, a scalar internal damage variable can be defined representing damage evolution. Such a damage variable has been defined in many research works27,28 directly based on inelastic shear strain, inelastic volumetic strain and confining pressure. This kind of damage variable increases by increasing inelastic shear strain or inelastic dilatant volumetric strain and it decreases by increasing confining pressure. On the basis of experimental investigations and microfracturing mechanism described above, since semi-brittle behavior and damage in geomaterials only occurs when stress state is located in dilatancy region, a parameter representing the accumulated amount of viscoplastic dilatancy

would be an appropriate measure for the amount of damage during creep. In this regard, Cristescu21 and Cristescu and Hunsche22–24 defined the internal damage variable based on inelastic volumetric strain work which is related to volume change due to opening and closing of microcracks. This parameter represents the energy of microcracking during creep of geomaterials which is accumulated during compressibility and released during dilatancy. For the definition of the damage parameter, this energybased parameter is more capable than those defined based on pure deformation in modeling real material behavior.30 For instance, it is capable of representing not only damage, when it decreases in dilatancy region, but also healing, when it increases in compressibility region. Also, experimental measurements demonstrate that long-term failure in geomaterials is encountered as soon as the defined scalar damage parameter reaches a critical damage value denoting the net energy release due to microfracturing until longterm failure is occurred. Therefore, in order to complete the definition of the long-term failure criteria, the critical damage value should be determined. Experimental evidence shows that this parameter increases by increasing confining pressure until a certain value and after that it is almost constant under different stress conditions.21 This behavior again implies that the definition of the damage variable based on the viscoplastic volumetric strain work, i.e. the accumulated value of the viscoplastic volumetric strain multiplied by mean stress, instead of only viscoplastic volumetric strain, is more practical and physically more justifiable. Because, using this method, the dimension of the damage variable is uniform to stress dimension which properly allows the definition of a critical damage value which is dependent on the confining stress. Another concept which is defined in this model is the shortterm failure criterion which is determined by the locus of critical deviatoric stresses at different confining stresses where failure of geomaterials occurs in relatively short time intervals. To sum up, in the present model, ductile and semi-brittle creep behavior of rock salt is simulated using an elasto-viscoplastic constitutive model, C–D boundary, short-term failure criterion and long-term or creep failure criterion. In this model, although a damage parameter is defined to determine long-term failure criteria, it does not have any effects on deformation and stress calculation as long as the long-term failure has not reached. For both ductile behavior in compressibility region and semi-brittle behavior in dilatancy region, the stress and deformation calculation is performed by the utilization of the elastio-viscoplastic constitutive model. Indeed, this model does not incorporate continuum damage in the sense that material properties of geomaterials, such as elastic and viscoplastic material parameters, are affected by the damage parameter. However, using the long-term failure criterion defined based on the scalar damage parameter, the transition between the secondary and tertiary stages of creep response can be identified. In this regard, it should be mentioned that the evolution of damage during viscoplastic flow leads to the evolution of localized damage if the accumulated damage approaches the long-term failure criterion.24 In these circumstances, creep damage accelerates due to the evolution of localized damage leading to dynamic crack growth and tertiary or accelerated creep. The present model is valid as long as viscoplastic flow dominates the rock salt behavior and continuum mechanics is valid. However, as short-term failure or long-term failure is encountered, this model is no longer valid and fracture mechanics approaches should be employed to analyze the post-failure behavior. It is also worth mentioning that, although the behavior of rock salt after the occurrence of longterm failure is not considered in the present paper, the post-failure behavior can be simulated using the defined damage parameter. To this end, based on the continuum damage method, the degradation of elastic and viscoplastic material properties after long-term failure can be considered as a function of the damage parameter.

S. Nazary Moghadam et al. / International Journal of Rock Mechanics & Mining Sciences 79 (2015) 250–267

One of the major capabilities of the present model is that it can be employed for the long-term safety assessment of underground storage caverns excavated in rock salt. In these underground structures, it is crucial to assess long-term performance in terms of both stability and adequate tightness against liquids and gases determined by permeability as a function of time. In this regard, the criteria defined for short-term and long-term failure can be utilized for stability analyses. Furthermore, in terms of tightness, the capabilities of the model to determine both dilatant area around the cavern and the amount of damage caused by dilatancy can be effectively employed. In this regard, on the one hand, research works31–33 have been carried out to develop models defining the evolution of the permeability as a function of the evolution of dilatancy. These models can be applied to analyze the seepage of fluids around the underground caverns based on the coupling of stress filed and seepage field. On the other hand, we know that when the stress sate is located in dilatancy domain, semi-brittle creep behavior of rock salt is encountered where increasing permeability and progressive damage leading to creep failure are inevitable. In contrast, as long as the stress state is located in the compressibility region, rock salt exhibits ductile creep behavior with decreasing damage and permeability. Therefore, C– D boundary can be considered as an appropriate long-term safety boundary to ensure the desired performance of repositories in rock salt for long periods of time. Generally, the time interval considered in an analysis determines the dominant components of the creep deformation. In this paper, since the time intervals considered in the analyses are in a range that steady state creep is dominant, transient creep of rock salt is not considered. In general, rock salt exhibits very complex behavior which is influenced by stress invariants, time, temperature, humidity, anisotropy, loading path and so on. However, considering all these effects in a single constitutive model is difficult. The main objective of the present paper is to develop a numerical framework for the time-dependent analysis of caverns excavated in rock salt using a constitutive model describing dilatancy and failure behavior and show the effect of this behavior in the cavern design. In this regard, by assuming that rock salt is isotropic, all constitutive functions depend on stress invariants while the direction of loading, i.e. principal axes directions of stress tensor, does not have any effects on the constitutive behavior. Also, the effect of humidity and temperature change is not considered in this model. Therefore, the focus of the present paper is on the effects of stress invariants and time. In this regard, in the research carried out by Hunsche30 on rock salt, it is shown that the compressibility–dilatancy is not affected by the third stress invariant, i.e. Lode angle representing triaxial compression-extension. However, since short-term failure is affected by third stress invariant,30 this parameter should be considered to model the behavior of rock salt near failure. But, in this paper, for the sake of simplicity, the effect of third stress invariant was not considered. This drawback can be improved by introducing additional terms in constitutive model representing the effects of third stress invariant based on available experimental data. Consequently, in this paper, it is assumed that the time-dependent behavior of rock salt depends on first stress invariant or mean stress, representing confining stress resulting in frictional behavior, and second stress invariant which is octahedral shear stress here, representing deviatoric stress causing distortion. In order to represent steady state creep together with inelastic volumetric changes and failure of rock salt, the constitutive equations are determined based on the creep models proposed by Cristescu,21 Jin and Cristescu25 and Paraschiv-Munteanu and Cristescu26 and most of the constitutive parameters are determined based on the experimental data obtained for New Mexico rock salt1,2 as well as Avery Island rock salt.34,35 In this

253

regard, Elastic material parameter and material parameters of the viscoplastic model are obtained based on the experimental data for New Mexico and Avery Island rock salt. However, due to the lack of data, material parameters of C–D boundary and failure criteria are calibrated based on the experimental data for Gorleben rock salt.21–26 In this regard, C–D boundary practically only depends on the stress state, i.e. mean stress and octahedral shear stress, and it is not dependent on the rock salt type.21 Also, experimental data show that the mentioned kinds of rock salt exhibit comparable mechanical properties, particularly failure responses are practically the same.21 Therefore, the proposed material properties are appropriate for preliminary studies although it would be better to determine all material properties based on a single rock salt. In this paper, the constitutive model is derived such a way that it can be used for large displacements analyses. In this regard, as the second Piola Kirchhof stress and Green's strain do not change under rigid body motions, any material description developed in terms of engineering stress and strain can be utilized in large displacements and large rotations but small strains analysis, provided that second Piola Kirchhof stress and Green's strain are used.36 Therefore, in the following, constitutive equations are written in terms of second Piola Kirchhof stress tensor :ij and Green's strain tensor , ij . Also, in this section, it is assumed that compressive stresses and strains are positive while tensile stresses and strains are negative. In the present creep model of rock salt, the constitutive functions are written in terms of the stress invariants, namely mean stress sm and octahedral shear stress τ which are defined in Appendix A. It is assumed that the total strain , ij can be decomposed into elastic and viscoplastic parts as

,ij = ,ije + ,ijvp , eij

(1) , ijvp

where and are the elastic and viscoplastic strain tensors, respectively. Also, the above equation can be rewritten in the rate form as follows: e vp ,̇ ij = ,̇ ij + ,̇ ij

(2)

where over-dot indicates the derivatives with respect to time. In this paper, since the transient creep is not considered, the viscoplastic strain rate is due to steady state creep and is defined by the following equation:

∂. vp ,̇ ij = k ∂:ij

(3)

where k is a viscosity coefficient and . is the viscoplastic potential function defined in Appendix A. The C–D boundary, which represented by the equation ? = 0, separates the area of a rock compressibility, i.e. ? > 0, and the area of the rock dilatancy, i.e. ? < 0, in the sm − τ plane. The C–D boundary has been defined as

?= −

⎛ s ⎞2 ⎛s ⎞ τ + f1 ⎜ m ⎟ + f2 ⎜ m ⎟ = 0 s⁎ ⎝ s⁎ ⎠ ⎝ s⁎ ⎠

(4)

where f1 ¼  0.017 and f2 ¼ 0.9. The area of the present constitutive model validity is bounded by the short term failure surface which is the locus of maximum values of octahedral shear stresses determined with standard triaxial tests in relative short time intervals. When the stress state is located out of this area, the shortterm failure of the rock salt is encountered. The short-term failure surface is represented by the following equation:

⎛ s ⎞ τf = γ1 − γ2 exp ⎜ −γ3 m ⎟ s⁎ ⎠ ⎝

(5)

where τf is the octahedral shear stress at failure and γ1, γ2 and γ3 are material parameters which are respectively equal to 38 MPa,

254

S. Nazary Moghadam et al. / International Journal of Rock Mechanics & Mining Sciences 79 (2015) 250–267 = D (t ) = − W vp (t )

(9)

where D (t ) is the damage parameter at time t . Also, the maximum allowable energy release due to microfracturing during dilatancy of the rock salt is defined as 2

(

= Df = − W vp

) failure

⎛ s ⎞3 = d0 ⎜ m ⎟ ⎝ s⁎ ⎠

(10)

where Df denotes the damage parameter at failure and d0 is equal to 0.0727 MPa. In addition, the parameter long-term failure ratio, i.e. LFR , is defined as the ratio of the current damage parameter to the failure damage parameter as follows:

LFR =

Fig. 1. Constitutive functions including viscoplastic potential, short-term failure surface, C–D boundary.

34.9 MPa and 0.04. The constitutive functions defined in the material model of rock salt, namely C–D boundary, short-term failure surface and viscoplastic potential function, are illustrated in Fig. 1 in mean stress Sm-octahedral shear stress τ space. As it can be seen from the figure, in the compressibility region, viscoplastic potential function . represents horizontal straight lines in sm − τ plane indicating that in this region no steady state viscoplastic volumetric change and therefore no viscoplastic compressibility is possible. In contrast, in dilatancy region, viscoplastic potential function . depends on Sm. Therefore, in dilatancy region, steady state viscoplastic volumetric change and therefore viscoplastic dilatancy occurs. Consequently, the viscoplastic potential function . describes either viscoplastic dilatancy or viscoplastic incompressibility. It is because of the fact that the steady state creep model must not cause the compressibility of the rock, otherwise the volume could shrink to zero due to steady state creep deformation. In addition, the parameter short-term failure ratio, i.e. SFR , is defined as the ratio of the current octahedral shear stress to the failure octahedral shear stress as follows:

SFR =

τ (t ) τf

(6)

Also, in order to define a failure criterion which considers the time-dependent failure of rock salt, an energetic failure criterion representing the progression in time of the damage has been used as follows21: = W vp (t )

=

vp vp ,̇ = = ,̇ ii

∫0

t

= sm ,̇ vp dt

D (t ) Df

(11)

In viscoplastic problems, as the stress–strain relationship is defined in a rate form, a numerical integration method should be employed to evaluate the stresses. With this aim, the increment of viscoplastic strain (δ , ijvp)n during the time interval (δt )n can be computed using the fully explicit or forward-Euler time integration scheme as n

vp n ij

( δ , ) = ( ,̇ ) (δt ) vp ij

n

(12)

where n represents the time step counter. As the total strain was decomposed into elastic and viscoplastic components, the stresses can be evaluated at the end of the time increment (δt )n as n + 1⎤ ( :ij)n + 1 = *ijkl ⎡⎣ ( ,kl )n + 1 − ( ,vp kl ) ⎦

(13)

where *ijkl is the elastic constitutive tensor. Using this procedure, stresses are always updated from the end of the last converged equilibrium state. Also, the incremental form of the above equation is given as

⎤ ( δ :ij)n = *ijkl ⎡⎣ ( δ ,kl )n − ( δ ,vp kl ) ⎦ n

(14)

It should also be mentioned that, since the employed explicit time integration scheme is conditionally stable, an allowable limit should be imposed on the time step length to assure stable and accurate solution. To this end, two criteria are defined to determine the proper time step size in each time step. The first criterion represented by the Eq. (15) limits the time step size so that effective viscoplastic strain increment in each time step is smaller than a predefined fraction of the total accumulated effective strain until that time step. Therefore, the parameter ξ is used to limit the vp n maximum value of effective viscoplastic strain increment (δ ,kl ) as n

( δ ,¯ vp)n = ( ,¯ ̇ vp)

n (δt )n ≤ ξ ( ,¯ )

(15)

(7) ,¯ = (8)

= where W vp stands for the volumetric viscoplastic strain work per = unit volume and ,̇ vp denotes the rate of volumetric viscoplastic = does not change strain. Based on the present creep model, the W vp in compressibility domain, while this parameter decreases in dilatancy domain representing the released mechanical energy due to opening and forming microcracks. Therefore, as long as the stress is located in the compressibility domain, damage does not occur. In contrast, as the stress state is located in dilatancy domain, damage begins to increase. On the basis of this argument, the total = is used as a measure of the rock salt damage as decrease of W vp follows:

3 ,ij ,ij 2

(16)

where ,¯ is the total effective strain. The minimum of the calculated values for (δt )n at each Gaussian integration points is selected to be checked by the second criterion. By the utilization of the second criterion which is represented by Eq. (17), an effort is made to prevent probable oscillations in the solution. In this criterion, the parameter ς is utilized to limit the change of time step length between successive time steps as

(δt )n + 1 ≤ ς (δt )n

(17)

Consequently, the minimum time step size calculated based on these two criteria is selected for the current time step size in the viscoplastic strain calculation. Using these restrictions, the suitable

S. Nazary Moghadam et al. / International Journal of Rock Mechanics & Mining Sciences 79 (2015) 250–267

length of time step can be chosen such a way that accurate results can be obtained. In the research works carried out by Zienkiewicz and Cormeau37,38 and Owen and Hinton,39 it is demonstrated that choosing a value of ξ in the range 0.01 ≤ ξ ≤ 0.02 and taking ς = 1.5 is suitable in order to achieve an accurate and stable solution. This explicit time integration scheme is simple and easy to be implemented. In particular, this scheme is very practical in case of complicated viscoplastic constitutive laws such as the one used in this paper. Also, the accuracy of this scheme has been proven in.37–39 Another choice to tackle the instability issues is to employ unconditionally stable implicit time integration schemes. These implicit methods provide robust numerical time integration schemes. However, in this paper, these implicit methods are not utilized and they can be considered in future studies.

contains displacements derivatives. The matrixes 3 (D) and 9 are defined in Appendix A. Also, D denotes the displacement gradient vector which is defined as follows:

⎧ ∂u ∂v ∂w ∂u ∂v ∂w ∂u ∂v ∂w ⎫ ⎬ DT = ⎨ , , , , , , , , ⎩ ∂x ∂x ∂x ∂y ∂y ∂y ∂z ∂z ∂z ⎭

In the analysis of continuum mechanics problems, nonlinear behavior may be encountered due to not only nonlinear constitutive models, but also large deformations. The nonlinear material model of rock salt along with the numerical integration of the constitutive equations has been discussed in the previous section. In this section, the large displacements viscoplastic behavior of rock salt is formulated using a large displacements finite element method. In the present numerical model, since the analyses are performed up to failure, the displacements are not so large that the distortion of finite element meshes is occurred. Hence, updating of finite element meshes is not necessary during the analyses. However, if post-failure behavior of rock salt is incorporated in the model, the displacements can be so large that an adaptive remeshing scheme would be necessary to update the finite element mesh of deformed body during analyses. Therefore, in the present paper, finite element mesh updating was not performed during analyses which can be considered in the case of post-failure analyses in future studies. In order to formulate large displacements continuum mechanics problems, a reference configuration should be used. In this paper, a total Lagrangian description, which uses the original undeformed continuum body as the reference configuration, has been utilized. Using the finite element discretization method, a domain is discretized into elements and the displacements in each element are related to nodal values via shape functions as

u = Np

(18)

where u is the displacement vector in an arbitrary place within the element, the matrix N contains interpolation functions 5i and p is the nodal displacement vector. For large displacements analyses, a strain measure which vanishes for any rigid body motion should be used. Green's strain meets this requirement and is suitable for Lagrangian description. Using Green's strain, the nonlinear relationship between strain and displacement can be decomposed into the small and large displacements components as40

, = ,l + , nl

,l = 9D

, nl =

1 3 (D) D 2

(19)

(20)

(21)

where , is the Green's strain vector, ,l is the linear component of Green's strain vector, , nl is the nonlinear component of Green's strain vector, and 3 (D) is a suitably defined matrix operator which

(22)

By taking variation from the equation derived for Green's strain and taking into account the relation δ 3 (D) D = 3 (D) δ D , the changes in Green's strain vector can be expressed as

δ , = 9δD +

1 1 δ 3 (D) D + 3 (D) δD = [9 + 3 (D)] δD 2 2

(23)

Also, the incremental form of the displacement gradient vector δ D can be determined in terms of nodal displacements as follows:

δ D = / δp 3. Finite element formulation

255

(24)

where / contains the derivatives of the interpolation functions 5i . From the above discussion it is obvious that the incremental form of the Green's strain can be expressed as

δ , = [9 + 3 (D)] / δp = ) (p) δp

(25)

where ) (p ) stands for the strain–displacement transformation matrix. As it can be seen, the strain–displacement transformation matrix relates the increments of Green's strain vector in an arbitrary point within an element to the increments of element nodal displacement vector. The matrix ) (p ) can be decomposed into the linear and nonlinear components as

) (p) = )l + )nl (p)

(26)

)l = 9/

(27)

)nl (p) = 3 (D) /

(28)

where )l is the linear strain–displacement transformation matrix which is the same as that of small displacements analysis, and )nl is the nonlinear strain–displacement transformation matrix which is the only term in the matrix ) (p ) that depends on the displacements. It is obvious that, regardless of small or large displacements assumption, the equilibrium condition should be satisfied. In order to satisfy the equilibrium condition, the principle of virtual work under virtual displacements can be applied. On the basis of this principle, the equilibrium condition for a continuum body is satisfied, provided that the virtual work of all forces acting on the body under a virtual displacement field which is compatible with boundary conditions is equal to zero. In the Lagrangian description, the conjugate vectors used to obtain the internal virtual work are the second Piola Kirchhoff stress vector S and the increment of the virtual Green's strain vector δ , v . Thus, using the total Lagrangian description, the equilibrium equation of a deforming continuum body at time t + Δt can be written as

W = Wint − Wext =

∫=

0

ST δ , v d =0 − - Text δpv = RT δpv = 0

(29)

where Wint is the internal virtual work, Wext is the external virtual work, δpv is an arbitrary small virtual displacement field which is compatible with the boundary conditions, -ext is the external force vector, and R represents the out-of-balance or residual force vector. In the above equation, the second Piola Kirchhoff stress vector is measured at time t + Δt using the initial undeformed configuration with volume =0 . Substituting the increment of the virtual Green's strain vector from Eq. (25) into the above equation gives

256

S. Nazary Moghadam et al. / International Journal of Rock Mechanics & Mining Sciences 79 (2015) 250–267

⎡ W = Wint − Wext = ⎢ ⎣

∫=

⎤ ST ) (p) d =0⎥ δpv − - Text δpv = RT δpv ⎦

0

(30)

=0

Rearranging the above equation, the equilibrium equation is rewritten as

(31)

R = -int − -ext = 0

-int =

∫=

) (p)T S d =0

(32)

0

where -int is the internal force vector. The above nonlinear equilibrium equations are solved using a time marching method which linearizes the equilibrium equations with respect to time based on the truncated Taylor series expansion. In these circumstances, during the time step (δ :)n , the incremental form of equilibrium equations is as follows:

(R)n + 1 = (R)n + (δ R)n = 0 n

(33)

n

(δ R)n = ( δ -int ) − ( δ -ext )

( δ -int )n = ∫=

0

written to implement the computational formulation achieved based on the proposed large displacement viscoplastic material model. In this section, using the computer code, a series of finite element analyses is performed to investigate the performance of underground storage caverns in rock salt as a function of different parameters representing cavern geometry and non-salt interbed effects. In these simulations, the underground caverns are modeled using axisymmetric eight-noded isoparametric elements in which 2  2 Gaussian integration rule is utilized for the numerical integration process of stiffness matrix computation.

(34)

⎡ δ ) (p)T ⎤n (S)n d = + 0 ⎣ ⎦

∫=

0

⎡ ) (p)T ⎤n (δ S)n d = 0 ⎣ ⎦

(35)

n

4.1. Verification of the proposed finite element formulation In this section, in order to verify the proposed finite element formulation, the stress distribution around a plain strain thickwalled cylinder was modeled using the proposed formulation and the results were compared with those obtained by an analytical solution. In this regard, if Norton creep law is employed, there is an analytical solution for creep response of thick-walled cylinders under internal or external pressure.41,42 Norton creep law is usually utilized to model steady sate creep in metals and it can be defined based on the Green-Lagrange strain and second PiolaKirchhoff stress as follows:

In the above equation, (δ :) can be substituted from the matrix form of Eq. (14) as

N ,¯ ̇ vp = A :¯

n (δ S)n = C ⎡⎣ () (p))n (δp)n − ,̇ vp (δt )n⎤⎦

¯ is effective where ,¯ ̇ vp is effective viscoplastic strain rate and : stress which are defined as

( )

(36)

where C is the elastic modular matrix. Also, based on the definition of the matrix / , the following expression can be achieved:

⎡ δ ) (p)T ⎤n (S)n = / T S˜ n / (δp)n ⎣ ⎦

()

(37)

where S˜ is a 9 × 9 matrix which contains stress components and is defined in Appendix A. Based on Eqs. (33–36), the displacement increment at the end of time increment (δt )n is achieved as follows: n

( Kt )

(δp)n = (δ -)n

(38)

In this equation, the time step counter n starts from a converged solution at the previous load increment and δ - and Kt are the incremental pseudo-load vector and the tangential stiffness matrix, respectively, which are defined as follows:

δ- =

∫=

) (p)T C,̇ vp δt d =0 + δ -ext + R 0

K t = K) + K S

K) =

∫=

) (p)T C) (p) d =0 0

,¯ ̇ vp =

(39)

(40)

(41)

:¯ =

(43)

3 ,ijvp ,ijvp 2

(44)

3 :ij :ij 2

(45)

Norton creep law is also widely used to simulate steady state creep of rock salt. In this regard, WIPP steady state creep model of rock salt proposed by Herrmann et al.1,2 is developed based on Norton creep law. WIPP steady state creep model of rock salt can be defined based on the Green-Lagrange strain and second PiolaKirchhoff stress tensors by the following equation: Q ⎛ 3 ⎞ N vp ,̇ ij = DS¯ e− RT ⎜ Sij ⎟ ⎝ 2S¯ ⎠

(46)

where D and N are material constants, R is the universal gas constant, Q is the activation energy and T is the temperature in Kelvin (K). The material parameters for WIPP steady state creep model determined for New Mexico rock salt1,2 are presented in Table 1. Using the WIPP steady state creep model, the stress distribution around a plane strain thick-walled cylinder with the inner radius of a and outer radius of b can be analytically obtained for long periods of time by the following equation41,42: 2

KS =

∫=

/ T S˜ / d =0 0

Sr = − Pi

4. Parametric study In this research, a finite element C þ þ computer code was

2

(b/a) N − 1

(42)

in which K) represents the linear stiffness matrix along with the initial displacement stiffness matrix, while KS denotes the initial stress stiffness matrix.

(b/r ) N − 1 (47)

Table 1 Material parameters for WIPP steady state creep model of the rock salt. Parameter

D (MPa  M/day)

M

Q (Cal/Mol)

R (Cal/Mol K)

T (K)

Value

0.1257

4.9

12,000

1.987

304

S. Nazary Moghadam et al. / International Journal of Rock Mechanics & Mining Sciences 79 (2015) 250–267

257

Fig. 2. Verification of the numerical formulation with analytical solution; (a) finite element mesh of the plain strain thick-walled cylinder, (b) variation of radial stress along line AB, and (c) variation of tangential stress along line AB.

Fig. 3. Verification of the proposed viscoplastic model; (a) geometry, finite element mesh and boundary condition of underground cavern excavated in rock salt, (b) variation of radial stress along line y¼ 0 after 6 years, (c) variation of tangential stress along line y¼ 0 after 6 years, and (c) variation with time of relative closure, i.e. cavern volume change divided by initial cavern volume.

258

S. Nazary Moghadam et al. / International Journal of Rock Mechanics & Mining Sciences 79 (2015) 250–267

2 2 (b/r ) N N 2 (b/a) N − 1

(

1− 1− Sθ = − Pi

)

(48)

where Pi is internal pressure, Sr is radial stress, Sθ is tangential stress and r is radial coordinate which is a ≤ r ≤ b . With the aim of verification, WIPP steady state creep model was implemented into the proposed large displacement finite element formulation to simulate the time-dependent stress distribution around a plain strain thick-walled cylinder subjected to internal pressure and the numerical analysis results were compared to those obtained by the analytical solution. To this end, due to the axial symmetry of the cylinder, an axisymmetric finite element model was utilized. The finite element mesh of the model is shown in Fig. 2(a). In this Table 2 Elastic material properties for rock salt and non-salt formations. Material

Elastic modulus (MPa)

Poisson's ratio Density (kg/m3)

Rock salt overburden and interbeds underlying rock formation

31 22

0.38 0.25

2200 2500

50

0.20

2500

model, in order to take into account the plain strain condition of the cylinder, the boundaries of the model on lines EF and CD are restrained against vertical displacement. The inner and outer radii of the cylinder is equal to 2.5 m and 7.5 m, respectively, while the height of the cylinder is equal to 1 m. Also, a pressure of Pi ¼  2 MPa was applied to the internal boundary of the cylinder. It should also be mentioned that in order to prevent instability in numerical results, the time stepping control parameters are chosen as ξ = 0.01 and ς = 1.5. Using the finite element model, stress distribution for different periods of time was numerically calculated and compared with the analytical solution in Fig. 2(b, c). As it can be seen from the figures, the numerical results obtained using the finite element formulation is stable and agrees well with the analytical solution. In the next stage, time dependent behavior of an underground cavern excavated in rock salt was numerically simulated by the finite element formulation using both the proposed viscoplastic model and the verified WIPP model. The geometry, finite element mesh and boundary conditions of the model are illustrated in Fig. 3(a). In these simulations, the bottom of the model is restrained against vertical displacements and the gravity loading and internal pressure is applied to the model at time t¼ 0. Also, it is assumed that the cavern is instantaneously excavated at time t¼0 and the finite element analyses are carried out over a period of 30 years. In order to apply initial condition, it is assumed that

Fig. 4. Schematic cavern initial geometries and boundary conditions for studying the effect of cavern geometry.

Table 3 Simulation details for studying the influence of cavern geometry (cavern internal pressure ¼ 5 MPa). Simulation number

Cavern diameter D (m)

Cavern height H (m)

H/D

D1/D2

Relative closure

1 2 3 4 5 6

120 118.554 118.554 127.519 173.070 109.027

240 237.108 237.108 255.038 173.070 327.082

2 2 2 2 1 3

1 0.5 2 1 1 1

0.068 0.077 0.063 0.067 0.060 0.075

(at (at (at (at (at (at

time time time time time time

27.40 years) 27.22 years) 30.14 years) 30.17 years) 30.17 years) 30.14 years)

S. Nazary Moghadam et al. / International Journal of Rock Mechanics & Mining Sciences 79 (2015) 250–267

259

Fig. 5. LFR contours in the case of studying the effect of cavern side wall slope.

the vertical in situ stress is obtained based on overburden pressure R as sz ¼ γ dz where γ is the density of geomaterials at different depths and z represents depth. Also, maximum and minimum horizontal in situ stresses respectively denoted by sH and sh can be generally obtained as sH ¼ n1sv and sh ¼n2sv, where n1 and n2 are coefficients of earth pressure in maximum and minimum horizontal in situ stress directions, respectively. In the present paper, for simplicity, it is assumed that in situ stress state before excavation of the underground opening was hydrostatic, i.e. the primary stresses were equal in all directions and n1 ¼n2 ¼1. This assumption is logical for thick salt deposits and great depths. This

assumption implies that, over a long period of time, rock salt behavior is the same as highly viscous fluid with no shear strength. However, in the case of caverns in thin salt layers, the coefficients n1 and n2 should be chosen more carefully since in this case the assumption of hydrostatic in situ stress is usually invalid. In this regard, it should be mentioned that when in situ stress state is not hydrostatic, due to the existence of shear stress in far field, the stress state around the opening is more critical in terms of cavern stability compared with the case of hydrostatic in situ stress sate. As in situ stress state before excavation of the cavity is assumed to be hydrostatic, the boundary conditions at the outer surface of the

Fig. 6. SFR contours in the case of studying the effect of cavern side wall slope.

260

S. Nazary Moghadam et al. / International Journal of Rock Mechanics & Mining Sciences 79 (2015) 250–267

Fig. 7. SFR contours in the case of studying the effect of cavern floor and roof shape.

model are imposed by applying distributed lateral edge loading which is equal to the hydrostatic stress in the rock formation.24 In addition, three different rock formations, namely rock salt, nonsalt overburden and non-salt underlying rock are considered in the analyses. According to the experimental data presented for Avery Island rock salt4 and for non-salt formations,16,43 elastic modulus, Poisson's ratio and the density of the rock salt along with non-salt formations are given in Table 2. It should also be mentioned that, the non-salt formations

exhibit significant creep behavior even though these formations creep at lower rates in comparison with rock salt. Therefore, for simplicity, the creep model used for rock salt with second order of magnitude lower steady state creep rate compared with that of rock salt is utilized to represent the time-dependent behavior of the non-salt formations.16 The analyses results for the underground cavern using the proposed viscoplastic material model is presented and compared with the verified WIPP model in Fig. 3(b–d). In these analyses, the time stepping control

Fig. 8. LFR contours in the case of studying the effect of cavern floor and roof shape.

S. Nazary Moghadam et al. / International Journal of Rock Mechanics & Mining Sciences 79 (2015) 250–267

261

Fig. 9. SFR contours in the case of studying the effect of cavern height-to-diameter ratio.

Fig. 10. LFR contours in the case of studying the effect of cavern height-to-diameter ratio.

parameters of ξ = 0.01 and ς = 1.5 are utilized. As it can be seen from the figures, the time integration scheme employed in this paper is stable for the simulation of the time-dependent behavior of underground caverns using the mentioned values for the time stepping control parameters. Also, it can be seen that both the proposed viscoplastic model and the verified WIPP model follows similar trend which implies the validity of the formulation of the proposed viscoplastic material model. According to the figures, stress relaxation reflected by the proposed viscoplastic model is more considerable than that represented by WIPP model. In particular, since the proposed viscoplastic material model takes into account the inelastic volume change during viscoplastic flow, it shows significant stress relaxation in terms of confining stress

which causes the occurrence of short-term failure, i.e. SFR greater than or equal to 1, in one of the Gaussian integration points around the cavern after 6 years. In contrast, since no inelastic volume change is considered in WIPP model, and since no failure criteria have been defined in WIPP model to capture the failure of the rock salt, numerical analyses using WIPP creep model can be performed for long periods of time producing extreme time-dependent deformations without any indication of failure response.16 It is due to the fact that WIPP creep model is mainly developed based on the creep models of metals and it cannot simulate semi-brittle creep behavior which is an inherent time-dependent behavior of geomaterials, e.g. rock salt. Also, using both the LFR parameter defined based on the time-dependent inelastic volume change and the

262

S. Nazary Moghadam et al. / International Journal of Rock Mechanics & Mining Sciences 79 (2015) 250–267

change in the present numerical formulation, this model can be better verified by comparing its results with experimental data. 4.2. Influence of cavern geometry

Fig. 11. Schematic cavern initial geometry and boundary conditions for studying the effect of non-salt interbeds.

Table 4 Simulation details for studying the influence of non-salt interbeds (cavern internal pressure¼5 MPa). Simulation number

Interbed location

Relative closure

1 2 3 4

Line Line Line Line

0.058 0.065 0.065 0.058

AB CD EF AB, CD and EF

(at (at (at (at

time time time time

30.21 30.21 30.21 30.21

years) years) years) years)

concept of C–D boundary, it is possible to predict the area around the cavern which is susceptible to the seepage of the stored material in the cavern. Also, as it can be observed from the figures, due to the inelastic volume change considered in the proposed viscoplastic model, this model exhibits more significant relative closure, i.e. cavern internal volume change divided by the initial cavern internal volume, compared with that predicted by WIPP model. It should also be mentioned that finite element analyses with different mesh refinements were performed where the mesh size near the cavern was small and it increased gradually by increasing the distance from the cavern. In this regard, finite element analysis using the mesh illustrated in Fig. 3(a) was accurate enough for the purpose of this paper. According to the finite element results presented in Fig. 3, the stress distribution around the cavern is smooth which means that the finite element mesh near the cavern surface is fine enough. Also, in terms of SFR and LFR contours around the cavern, no significant change was observed by refining finite element mesh more than that utilized in the present analyses. Finally, it is worth mentioning that by implementing a transient creep model considering inelastic volume

In this section, finite element analyses have been performed to investigate the influence of geometric parameters on the stability and deformation of underground caverns excavated in rock salt. To this end, the underground caverns with different geometries which are schematically illustrated in Fig. 4 as well as Table 3 have been analyzed numerically. In these analyses, six different geometries are assumed for the cavern layout in which cavern side wall slope, cavern floor and roof shape and cavern height-to-diameter ratio are varied in a constant cavern internal volume condition. Then, the individual influence of the mentioned parameters in the performance of the underground cavern is investigated. In the first series of analyses, the potential influence of side wall slope on the cavern performance is investigated by comparing the finite element results of Simulations 1, 2 and 3. In Simulation 1, the cavern side wall is assumed to be vertical while in Simulation 2 the cavern side wall slope is assumed such a way that the center of cavern internal volume moves towards the cavern floor in comparison with Simulation 1. Also, in Simulation 3, the side wall slope is chosen such a way that the center of cavern internal volume moves towards the cavern roof compared with that of Simulation 1. In these simulations, in order to diminish the effect of stress concentration, the cavern corners are modeled by fillets. The simulation results comprising of the SFR and LFR contours are plotted in Figs. 5 and 6. Also, the relative closure at the last time step of the analyses is reflected in Table 3. As it can be observed, due to the stress concentration effects, the stress state in cavern corners is more critical compared to the neighboring areas. Also, from the figures, it can be concluded that the closer the center of cavern internal volume to the cavern floor, the more critical the stress state around the cavern, the faster the cavern relative closure, and vice versa. For the time interval considered in these analyses, the stress state in the case of Simulation 3 do not endanger the stability of the cavern, while in the case of Simulations 1 and 2 short-term failure is encountered at the corner of the cavern floor at time around t ¼27 years leading to cavern instability. As a consequence, a center of cavern internal volume closer to the cavern roof results in a more stable cavern. In the next analyses, the influence of cavern floor and roof shape on the performance of the cavern is investigated. To this end, in simulation 4 the initial geometry of the cavern is assumed to be comprised of a vertical cylinder and two hemispheres at top and bottom of the cavern, and then the results of Simulations 1 and 4 are compared. The SFR and LFR contours are plotted in Figs. 7 and 8, and also the relative closure at the last time step of the analyses is reflected in Table 3. From the SFR and LFR contours, it can be seen that the spherical shape of the cavern roof and floor not only brings about less critical stress state, but also considerably diminishes the viscoplastic dilatancy around the cavern. In addition, spherical roof and floor reduces the rate of cavern relative closure. It is also worth mentioning that, unlike the Simulation 1, the stress state in the case of Simulation 4 does not endanger the stability of the cavern in the time interval considered in these analyses. The effect of the cavern height-to-diameter ratio on the performance of the cavern is studied in the last series of analyses. In Simulations 4, 5 and 6, the radius and height of a cavern with spherical floor and roof are varied while the cavern volume is kept constant. The SFR and LFR contours are plotted in Figs. 9 and 10 and the relative closure at the last time step of the analyses is reflected in Table 3. According to the SFR and LFR contours, the stress state becomes

S. Nazary Moghadam et al. / International Journal of Rock Mechanics & Mining Sciences 79 (2015) 250–267

263

Fig. 12. Finite element configuration of the cavern for studying the effect of non-salt interbeds.

closer to the failure state and the viscoplastic dilatancy increases as the height-to-diameter ratio deviates from 1, although the stress states do not bring about the instability of the caverns in the case of Simulations 4, 5 and 6 in the time period considered in these analyses. Also, it can be seen that caverns with height-todiameter ratios closer to 1 show slower relative closure. Based on the analyses, as the LFR contour is equal to zero in the case of simulation 5, i.e. spherical cavern, it can be concluded that the

stress state in the rock salt is in the compressibility domain for a spherical cavern. As a result, caverns with geometries closer to spherical shape, i.e. height-to-diameter ratios closer to 1, are more stable, and also the spherical shape is the optimum geometry for the underground caverns. Finally, it should also be mentioned that as it can be observed in the SFR and LFR contours, for every simulation the most critical stress state occurs in cavern floor where the greatest difference between internal pressure and geostatic

264

S. Nazary Moghadam et al. / International Journal of Rock Mechanics & Mining Sciences 79 (2015) 250–267

Fig. 13. SFR contours in the case of studying the effect of non-salt interbeds.

pressure is encountered. 4.3. Influence of non-salt interbeds In this section, the influence of non-salt interbeds on the stability and deformation of salt caverns is studied. In order to assess the effects of non-salt interbeds number and location on the performance of the caverns, four different cases of interbed locations are considered. The thicknesses of the interbeds are assumed to be equal to 10 m in all cases, and the number and locations of the

interbeds are illustrated in Fig. 11 as well as Table 4. The boundary and initial conditions and the constitutive parameters are the same as those represented in previous section and analyses are performed in a t¼30 years period. The finite element mesh of the cavern for the simulations presented in this section is depicted in Fig. 12, and also the SFR and LFR contours are plotted in Figs. 13 and 14. The relative closure at the last time step of the analyses is reflected in Table 4. As can be observed, non-salt interbeds which are not located near the cavern surface do not affect the cavern relative closure, as

S. Nazary Moghadam et al. / International Journal of Rock Mechanics & Mining Sciences 79 (2015) 250–267

265

Fig. 14. LFR contours in the case of studying the effect of non-salt interbeds.

in the case of Simulations 2 and 3, while the presence of non-salt interbeds near the cavern surface, as in the case of Simulations 1 and 4, brings about lower relative closure. From the SFR contours, in every simulation it can be observed that beside the most critical stress state located at the bottom of the cavern, a relative critical stress state spreads across the non-salt interbeds for long distances which indicates the possible danger of the stored material leakage for long distances from the cavern for a more severe

stress state compared to the one assumed herein. As it can be seen, the mentioned issue is more considerable when the non-salt interbed is in contact with the cavern surface. However, for the time interval considered in these analyses, the stress states do not endanger the stability of the cavern. Moreover, LFR contours show the highest values at the cavern floor in the case of Simulations 2 and 3, while in the case of Simulations 1 and 4 the highest values of LFR contour occur not only at the cavern floor, but also around

266

S. Nazary Moghadam et al. / International Journal of Rock Mechanics & Mining Sciences 79 (2015) 250–267

the intersection point of the cavern surface and non-salt interbed. As a consequence, the presence of non-salt interbeds which creep at a lower rate compared to the rock salt has a considerable adverse effect on the stability of the caverns, particularly when the interbeds are in contact with the cavern surface.

5. Conclusions In the present paper, attempts have been made to propose a systematic approach to assess the time-dependent performance of underground caverns excavated in rock salt. With this aim, a computational model based on an elasto-viscoplastic creep model reflecting dilatancy and failure of rock salt has been formulated in a Lagrangian large deformation finite element framework. Using the computational model, a series of finite element analyses have been performed to determine the influence of cavern geometric parameters, namely cavern side wall slope, cavern floor and roof shape and cavern diameter-to-thickness ratio on the cavern performance. According to the analyses, a cavern side wall slope representing a center of cavern internal volume closer to the cavern roof brings about a more stable cavern. Also, it was demonstrated that the spherical shape of the cavern roof and floor not only brings about less critical stress state, but also considerably diminishes the viscoplastic dilatancy around the cavern. Finite element analyses results revealed that the stress state becomes closer to the critical state and the viscoplastic dilatancy increases as the height-to-diameter ratio deviates from 1. As a result, caverns with geometries closer to spherical shape, i.e. height-to-diameter ratios closer to 1, are more stable, and also the spherical shape is the optimum geometry for the underground caverns. Furthermore, the influence of non-salt interbeds which creep at a lower rate compared to the rock salt on the performance of the caverns was investigated. The finite element analyses results revealed that a relative critical stress state spreads across the non-salt interbeds for long distances indicating the possible danger of the stored material leakage for long distances in a severe stress state. The mentioned issue was more considerable when the non-salt interbed was in contact with the cavern surface.

The authors of this paper would like to appreciate Professor N. D. Cristescu from the University of Florida and Professor M. Aubertin from École Polytechnique de Montréal for their useful discussions and suggestions during the preparation of this paper.

The stress invariants used to define constitutive equations are as follows:

(A1)

:′ij = :ij − sm δij

τ=

1 3

:′ij :′ij

(A4)

where b, m and n are material parameters which are equal to  0.864  10  9 1/day, 5 and  0.1, respectively, and the function 8 is defined as l p⎛ τ⎞ ∂8 = ⎜ ⎟ s⁎ ⎝ s⁎ ⎠ ∂τ

(A5)

where p ¼0.3677184  10  7 MPa/day, l ¼5 and sn ¼ 1 MPa. Furthermore, τm is the maximum value of octahedral shear stress of the compressibility–dilatancy (C–D) boundary, and also parameters sa and sb are obtained by intersecting the C–D boundary with a horizontal straight line in sm–τ plane for a fixed value of τ o τm. Therefore, the parameters τm, sa and sb can be achieved as

f 2 τm = − 2 s⁎ 4f1

sb = s⁎

−f2 +

(A6)

(f

2

2

+

1 4f1 τ 2 s⁎

)

2f1

−f2 −

(f

2

2

(A7)

+

1 4f1 τ 2 s⁎

)

2f1

(A8)

Also, the matrixes 3 (D), 9 , and S˜ are defined by the following equations:

Appendix A

:ii 3

⎧ ⎛ ⎞m ⎛ ⎞n + 1 ⎪ bs⁎ ⎜ τ ⎟ ⎜ sm ⎟ IF: ( τ ≥ τm ) ⎪ n + 1 ⎝ s⁎ ⎠ ⎝ s⁎ ⎠ ⎪ or ( τ < τm and sm ⎛τ⎞ ⎪ + 8⎜ ⎟ ⎪ s ≤ a) ⎝ s⁎ ⎠ ⎪ ⎪ m + n 1 bs⁎ ⎛ τ ⎞ ⎛ sa ⎞ ⎪ IF: ( τ < τm and sa < sm < sb ) ⎪ n + 1 ⎜⎝ s⁎ ⎟⎠ ⎜⎝ s⁎ ⎟⎠ ⎪ ⎛τ⎞ ⎪ ⎪ + 8⎜ ⎟ k. = ⎨ ⎝ s⁎ ⎠ ⎪ ⎪ ⎛ ⎞m ⎪ bs⁎ ⎜ τ ⎟ IF: ( τ < τm and sm ≥ sb ) ⎪ n + 1 ⎝ s⁎ ⎠ ⎪ ⎪ ⎡ ⎛ s ⎞n + 1 ⎛ s ⎞n + 1 ⎪ ⎢ ⎜ a⎟ − ⎜ b⎟ ⎝ s⁎ ⎠ ⎪ ⎢⎣ ⎝ s⁎ ⎠ ⎪ ⎪ ⎛τ⎞ ⎛ s ⎞n + 1⎤ ⎪ + ⎜ m ⎟ ⎥ + 8⎜ ⎟ ⎥⎦ s ⎪ ⎝ s⁎ ⎠ ⎝ ⎠ ⁎ ⎩

sa = s⁎

Acknowledgements

sm =

used in this paper was defined as

(A2)

(A3)

where S′ij is the deviatoric stress tensor and δ ij denotes the Kronecker delta. Also, the viscoplastic potential for the creep model

⎡ ∂u ⎢ ⎢ ∂x ⎢ 0 ⎢ ⎢ ⎢ 0 ⎢ 3 (D) = ⎢ ∂u ⎢ ⎢ ∂y ⎢ ⎢ 0 ⎢ ⎢ ∂u ⎢⎣ ∂z

∂v ∂w ∂x ∂x

0

0

0

0

0

0

0

∂u ∂v ∂w ∂y ∂y ∂y

0

0

0

∂u ∂v ∂z ∂z

∂v ∂w ∂u ∂y ∂y ∂x ∂u 0 0 ∂z ∂v ∂w 0 ∂z ∂z

0

0

0

0

∂v ∂w 0 ∂x ∂x ∂v ∂w ∂u ∂z ∂z ∂y ∂u 0 0 ∂x

0

∂v ∂y ∂v ∂x

⎤ 0 ⎥ ⎥ 0 ⎥ ⎥ ⎥ ∂w ⎥ ∂z ⎥ ⎥ 0 ⎥ ⎥ ∂w ⎥ ⎥ ∂y ⎥ ∂w ⎥ ⎥ ∂x ⎦

(A9)

S. Nazary Moghadam et al. / International Journal of Rock Mechanics & Mining Sciences 79 (2015) 250–267

⎡1 ⎢ ⎢0 ⎢0 9=⎢ 0 ⎢ ⎢0 ⎢⎣ 0

0 1 0 1 0 0

0 0 1 0 0 1

0 0 0 1 0 0

0 0 0 0 0 0

0 0 0 0 1 0

0 0 0 0 0 0

0 0 0 0 1 0

0⎤ ⎥ 0⎥ 1⎥ 0⎥ ⎥ 0⎥ 1⎥⎦

⎡ :xI :xyI :xz I ⎤ ⎥ ⎢ S˜ = ⎢ :xyI :yI :yz I ⎥ ⎥ ⎢ ⎣ :xz I :yz I :z I ⎦

(A10)

(A11)

where I is the 3 × 3 identity matrix.

References 1 Herrmann W, Wawersik WR, Lauson HS. Analysis of Steady State Creep of Southeastern New Mexico Bedded Salt. Sandia National Laboratories, SAND800558; 1980. 2 Herrmann W, Wawersik WR, Lauson HS. Creep Curves and Fitting Parameters for Southeastern New Mexico Bedded Salt. Sandia National Laboratories, SAND800558; 1980. 3 Aubertin M, Gill DE, Ladanyi B. A unified viscoplastic model for the inelastic flow of alkali halides. Mech Mater. 1991;11:63–82. 4 Yahya OML, Aubertin M, Julien MR. A unified representation of the plasticity, creep and relaxation behavior of rocksalt. Int J Rock Mech Min Sci. 2000;37:787– 800. 5 Shao JF, Zhu QZ, Su K. Modeling of creep in rock materials in terms of material degradation. Comput Geotech. 2003;30:549–555. 6 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:1237–1251. 7 Yang C, Daemen JJK, Yin JH. Experimental investigation of creep behavior of salt rock. Int J Rock Mech Min Sci. 1999;36:233–242. 8 Chang C, Zoback MD. Viscous creep in room-dried unconsolidated Gulf of Mexico shale (II): Development of a viscoplasticity model. J Petrol Sci Eng. 2010;72:50–55. 9 Dawson PR, Munson DE. Numerical simulation of creep deformations around a room in a deep potash mine. Int J Rock Mech Min Sci. 1983;20:33–42. 10 Munson DE, Fossum AF, Senseny PE. Approach to first principles model prediction of measured WIPP (waste isolation pilot plant) in-situ room closure in salt. Tunnel Undergr Space Technol. 1990;5:135–139. 11 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:629–638. 12 Fuenkajorn K, Phueakphum D. Effects of cyclic loading on mechanical properties of Maha Sarakham salt. Eng Geol. 2010;112:43–52. 13 Zhu WS, Li XJ, Zhang QB, Zheng WH, Xin XL, Sun AH, Li SC. A study on sidewall displacement prediction and stability evaluations for large underground power station caverns. Int J Rock Mech Min Sci. 2010;47:1055–1062. 14 Hatzor YH, Wainshtein I, Mazor DB. Stability of shallow karstic caverns in blocky rock masses. Int J Rock Mech Min Sci. 2010;47:1289–1303. 15 Zhu WS, Li Y, Li S, Wang S, Zhang Q. Quasi-three-dimensional physical model tests on a cavern complex under high in-situ stresses. Int J Rock Mech Min Sci. 2011;48:199–209. 16 Wang G, Guo K, Christianson M, Konietzky H. Deformation characteristics of rock salt with mudstone interbeds surrounding gas and oil storage cavern. Int J Rock Mech Min Sci. 2011;48:871–877.

267

17 Nazary Moghadam S, Mirzabozorg H, Noorzad A. Modeling time-dependent behavior of gas caverns in rock salt considering creep, dilatancy and failure. Tunnel Undergr Space Technol. 2013;33:171–185. 18 Tongtao Wang, Chunhe Yang, Xiangzhen Yan, Yinping Li, Wei Liu, Cheng Liang, Jie Li. Dynamic response of underground gas storage salt cavern under seismic loads. Tunnel Undergr Space Technol. 2014;43:241–252. 19 Deng JQ, Yang Q, Liu YR. Time-dependent behaviour and stability evaluation of gas storage caverns in salt rock based on deformation reinforcement theory. Tunnel Undergr Space Technol. 2014;43:241–252. 20 Khaledi K, Mahmoudi E, Datcheva M, König D, Schanz T. Sensitivity analysis and parameter identification of a time dependent constitutive model for rock salt. J Comput Applied Math. 2015 in press, 10.1016/j.cam.2015.03.049. 21 Cristescu ND. A general constitutive equation for transient and stationary creep of rock salt. Int J Rock Mech Min Sci. 1993;30:125–140. 22 Cristescu ND, Hunsche U. Determination of a nonassociated constitutive equation for rock salt from experiments. In: Proceeding of IUTAM Symposium, Berlin: Springer Verlag; 1992: 511–523. 23 Cristescu ND, Hunsche U. A constitutive equation for salt. In: Proceeding of 7th International Congress on Rock Mechanics. Rotterdam: Balkema; 1993: 18211830. 24 Cristescu ND, Hunsche U. Time effects in rock mechanics. New York: John Wiley & Sons; 1998. 25 Jin J, Cristescu ND. An elastic viscoplastic model for transient creep of rock salt. Int J Plast. 1998;14:85–107. 26 Paraschiv-Munteanu I, Cristescu ND. Stress relaxation during creep of rocks around deep boreholes. Int J Eng. Sci. 2001;39:737–754. 27 Chan KS, Bodner SR, Fossum AF, Munson DE. A constitutive model for inelastic flow and damage evolution in solids under triaxial compression. Mech Mater. 1992;14:1–14. 28 Chan KS, Bodner SR, Fossum AF, Munson DE. A damage mechanics treatment of creep failure in rock salt. Int J Damage Mech. 1997;6:121–152. 29 Chen Z, Wang ML, Lu T. Study of tertiary creep of rock salt. J Eng Mech ASCE. 1997;123:77–82. 30 Hunsche Udo, Hampel Andreas. Rock salt-the mechanical properties of the host rock material for a radioactive waste repository. Eng Geol. 1999;52:271–291. 31 Peach CJ. Influence of Deformation on the Fluid Transport Properties of Salt Rocks. Utrecht: Proefschrift, RUU; 1991. 32 Stormont JC, Daemen JJK. Laboratory study of gas permebility changes in rock salt during deformation. Int J Rock Mech Min Sci Geomech Abstr. 1992;29:325– 342. 33 Popp T, Kern H. Ultrasonic wave velocities, gas permeability and porosity in natural and granular rock salt. Phys Chem Earth. 1998;23:373–378. 34 Hansen FD, Mellegard KD. Creep of 50-mm Diameter Specimens of Some Salt from Avery Island. Louisiana. RE/SPEC ONWI-104; 1980. 35 Mellegard KD, Senseny PE, Hansen FD. Quasi-static strength and creep characteristics of 100-mm diameter specimens of salt from Avery Island. Louisiana. RE/ SPEC ONWI-250; 1981. 36 Bathe KJ. Finite Element Procedures. New Jersey: Prentice-Hall; 1996. 37 Zienkiewicz OC, Cormeau IC. Visco-plastic solution by finite element process. Arch Mech. 1972;24(5–6)873–889. 38 Zienkiewicz OC, Cormeau IC. Visco-plasticity – plasticity and creep in elastic solids – a unified numerical solution approach. Int J Numer Methods Eng, 8; 1974. p. 821– 845. 39 Owen DRJ, Hinton E. Finite Elements in Plasticity: Theory and Practice. Swansea: Pineridge Press; 1980. 40 Crisfield MA. Non-linear Finite Element Analysis of Solids and Structures. 1: Essentials. 1991. 41 Boyle JT, Spence J. Stress analysis for creep. London: Butterworths; 1983. 42 Skrzypek JJ, Hetnarski RB. Plasticity and Creep–Theory, Examples and Problems. Boca Raton, Fla: CRC Press; 1993. 43 Pfeifle TW, Mellegard KD, Senseny PE. Preliminary Constitutive Properties for Salt and Nonsalt Rocks from Four Potential Repository Sites. Technical report RESPEC Inc.; 1983.