# Semianalytical Solution for Large Deformation of Salt Cavern with Strain-Softening Behavior

Authors: Lina Ran; Huabin Zhang; Qingqing Zhang
Journal: Advances in Civil Engineering
 (2020)

---

## Abstract

A semianalytical solution of stress and displacement in the strain-softening and plastic flow zones of a salt cavern is presented. The solution is derived by adopting the large deformation theory, considering the nonlinear Hoek–Brown (H-B) strength criterion. The Romberg method is used to carry out numerical calculation, and then, the large deformation law of displacement is analyzed. The results are compared with those obtained by former numerical methods, and the solutions are validated. The results indicate that the displacement of the plastic zone decreases with the increase in distance away from the salt cavern. Similarly, it decreases with an increase in the geological strength index or running pressure, with the running pressure having a more significant effect on the displacement. It increases with the dilation angle, and the impact degree gradually increases. Compared with the softening parameter,h, of the plastic zone, the flow parameter, f, has little impact on the displacement. The displacement of the plastic zone obviously increased when considering the strain-softening of salt rock. When considering the shear dilation and softening behaviors of salt rock, the analytical solution obtained by employing the experiential regression Hoek–Brown (H-B) criterion, which considers many factors such as the structural characteristics of the salt formation and the rock mass quality, is safer and closer to the actual situation. This study can provide reference for many applications, including but not confined to analyzing the deformation of the surrounding rock of an underground salt cavern storage facility during construction.


---

## Body

## 1. Introduction

Salt caverns are considered the best places in the world for energy reserves [1–5]. The reason is that salt cavern storage facilities, compared to other underground storage facilities, have many advantages. For example, it has simple hydrological conditions, a wide distribution range, and a complete structure and stratum status. Moreover, the mechanical properties of rock salt, compact structure, low porosity and permeability, and damage self-recovery capability can ensure the safety of salt cavern. An underground salt cavern storage facility can be constructed using a water solution to form underground adits in deep salt formations. However, the redistribution of the stress and strain around the cavity and the damage to the mechanical properties of the surrounding rock are closely linked with the solution process. In recent years, seal failure, volume shrinkage, excessive deformation or heavy collapse, serious ground subsidence, etc. have become among the most arduous challenges to the engineers. For example, in the Hutchison salt cavern gas storage in Kansas, USA, a ground subsidence with 90 m in diameter was observed. Internal pressure caused excessive volume shrinkage, and the Eminence salt cavern storage in Mississippi, USA, lost about 40% of its volume in 2 years [6]. Thus, analysis of the stability of underground salt cavern gas storage is very important for controlling this hazard.In order to ensure the safety and stability of underground adits, a large number of analytical methods for the stress and deformation of the surrounding rock have been suggested for ideal underground spherical caverns in the past. For stability analysis of deep underground cavity, the Hoek–Brown (H-B) failure criterion has been found to be more suitable than the commonly used Mohr–Coulomb (M-C) criterion. The H-B failure criterion was proposed for hard rock masses and very poor qualities of rock masses. It can reflect the inherent nonlinear damage characteristics of rock masses. Using the nonlinear H-B failure criterion, Brown et al. [7] derived a solution for the stress and displacement of the plastic zone of a tunnel. Using the nonassociated flow rules and considering three cases for the elastic strain of the plastic zone, Park and Kim [8] put forward a solution for the displacement of the plastic zone of a circular tunnel. Employing the H-B criterion, Lee and Pietruszczak [9] analyzed the stress and deformation field of circular adits with strain-softening behavior and presented a solution for the surrounding rock of a tunnel by assuming that the criterion parameters were linearly transformed with the shear strain. Introducing the nonassociated flow rules in the design of a hydropower tunnel, Rojat et al. [10] derived a closed-form solution for the surrounding rock of a tunnel using the H-B criterion. Park [11] dealt with a similarity solution for a spherical or circular adits excavated in elastic-strain-softening rock mass compatible with the M-C or H-B yield criterion. The partial differential equations from stress equilibrium, constitutive law, and consistency equations are simplified with first-order ordinary differential equations that solved via the Runge–Kutta (R-K) method. In China, Liu and Lin [12] derived the elastoplastic deformation of a tunnel by adopting the M-C criterion. Adopting the H-B strength theory and limit equilibrium equation, Yao et al. [13] and Pan et al. [14] conducted elastoplastic analyses of the surrounding rock in an idealized roadway with hydrostatic pressure. Wang et al. [15] proposed a brittle-plastic solution of a spherical cavity, and the solution is used in the strain-softening analysis. They also proposed numerical solutions and studied their application to determining stress–strain curve of strain-softening H-B rock masses. Wang et al. [16] proposed an elastic analytical solution for a spherical salt cavern with a lateral pressure coefficient (λ) and used the H-B failure criterion to derive the internal pressure limit when the plastic failure of the gas storage cavern did not appear. In order to study excavation problems for tunnels exhibiting strain-softening behavior, considering the effect of the confining stress, Wen and Yang [17] and Cai et al. [18] conducted numerical calculations by adopting the R-K method. They obtained the deformation around a tunnel and analyzed the displacement of the surrounding rock of a deep underground chamber by using the H-B criterion. Using the elastoplastic theory, Wang et al. [19], Fan et al. [20], and Hu and Zhu [21] presented analytical solutions to the expansion of spherical cavities and analyzed the effects of softening and shear dilation on the spherical expansion. Adopting the M-C criterion and nonassociated flow rule and introducing the large deformation theory, Wang and Zhu [22] obtained an elastoplastic analytical solution to the cavity expansion. Zhang et al. [23] established a multistep brittle-plastic model by simplifying the strain-softening process of the postfailure region, and they proposed analytical solutions of a spherical cavern with considering the deterioration of elastic parameters. Both M-C and H-B criteria are included in their analysis. Zhang et al. [24] analyzed the stress and deformation of the surrounding rock in a spherical salt cavern under ideal elastoplastic conditions, and they predicted the creep deformation of the surrounding rock via theoretical analysis.The analytical elastoplastic solution for a space problem accounting for large deformation requires a complicated mathematical deduction, and, even when it is combined with a numerical method, there are a limited number of studies that can be referenced. Wang and Qian [25] proposed a finite difference method. Adopting the large deformation theory, Guan et al. [26] proposed a numerical procedure that reflects the strain-softening behavior and studied the instability problem of a plane strain axially symmetric tunnel. The above studies mainly involved two-dimensional research on circular tunnels. Jiang et al. [27] provided key information for predicting deformational performance for a tunnel undergoing large deformation problem. Zhang et al. [28] presented a large strain numerical solution of a circular tunnel for the ground reaction curve in strain-softening rock masses. In consideration of the tunnelling-induced large deformation, Mo et al. [29] analyzed the cavity contraction in M-C and H-B media.In the previous results for salt cavern storage projects, there has been little information on the large deformation [30] properties of salt rock, and these results were obtained mainly by employing the M-C criterion. Accounting for the erosion effect of solution mining on the cavity wall, the H-B criterion, which is based on empirical regression, is especially practical with salt cavern gas storage. Hence, this paper considers the large deformation around a salt cavern by introducing the H-B criterion and taking the strain-softening property of salt rock into account. This study can have significant implications for the design and construction of a deep hollow spherical domain.


## 2. Mechanics Model for Spherical Salt Cavern Gas Storage

In situ investigations indicated that salt caverns have all irregular shapes; they are different from each other. Most of the salt storage caverns are pear-like in the design stage. However, in order to facilitate the theoretical analysis and derivation, we can establish a spherical symmetry model for the stability analysis of a salt cavern under hydrostatic stress (side pressure coefficient is 1). Salt cavern gas storage is here regarded as a deep hollow spherical domain subjected to internal and external pressures, the radius of the salt cavern is designated asR0. The gas pressure (approximately homogeneous running pressure), pi, and hydrostatic pressure, p0, are the inner and outer boundaries of the model, respectively. Gravity and side pressure coefficient are not considered in the model. In order to meet the demands of spherical symmetry, the coordinate origin of the salt cavern mechanics model is designated in the center of the spherical adits. It agrees with the hypothesis of isotropy during the construction of a salt cavern. There exist elastic and plastic zones around the salt cavern. The far filed of the cavern (R2⟶∞) enters into an elastic state. The stress of the elastic zone is smaller than the strength of the salt rock; as a result, the deformations can recover after unloading. At r=R2, the radial stress and tangential strain are designated as σR2 and εR2, respectively. The near-zone of the cavern (R1⟶R2) enters into a strain-softening state. The rock in the strain-softening zone will undergo microcrack propagation and plastic deformation. At r=R1, the radial stress and tangential strain are designated as σR1 and εR1, respectively. Likewise, the near-zone of the cavern (R0⟶R1) enters into a plastic flow state. The rock in the plastic flow zone undergoes ultimate failure. What is more, the salt rock around the cavity satisfies the continuity conditions both in the elastoplastic interface, r=R2, and in the interface of the strain-softening and plastic flow zones, r=R1. Figure 1 shows the elastoplastic model of the rock surrounding spherical salt cavity.Figure 1
Elastoplastic model of spherical salt cavern. (a) Sketch of salt cavern deformation zones. (b) Simplified model of salt rock with strain-softening behavior.
(a)(b)The rock salt formation where salt cavern gas storage is located has a depth arranging from 1000 to 2000 m, so surrounding rock of the salt cavern is situated in high geostress levels. In situ investigations and compression tests indicated that the large plastic deformation played an unneglectable role during the analysis of the mechanical properties of rock salt under triaxial stress. Thus, the shear dilatation and strain-softening properties of rock salt [30–33] are considered in this study, and a strain-softening mechanical model is used.


## 3. Strength Criterion

The generalized H-B criterion is an experiential regression strength criterion. The H-B failure criterion was proposed in 1980s originally as follows [34]:(1a)σ1=σ3+σcmiσ3σc+11/2,where σc is the uniaxial compressive strength of the intact rock (MPa) and mi is a dimensionless material constant, with the value of mi varying from 0.001 to 25.0.To eliminate some of the deficiencies with the application of the original H-B failure criterion, a modified rock mass failure criterion has been developed and is presented by Hoek et al. [35]. Hoek et al. suggested a specific form of the strength criterion:(1b)σ1=σ3+σcimbσ3σci+sa,where σci is the uniaxial compressive strength of the intact rock (MPa) and mb, s, and a are dimensionless material constants.The H-B failure criterion was updated by discussing the use of the geological strength index (GSI) to describe the structure and condition of a wide range of types of rock masses and the selection of the disturbance factorD [36, 37]. The dimensionless material constants mb, s, and a can be calculated by (2) [36, 37]. The parameter a varies over a range from 0.5 to 0.666, where a=1/2+1/6e−GSI/15−e−20/3. Considering the difficulty to obtain theoretical solutions of nonlinear numerical problems, to simplify the calculation, the commonly used H-B strength criterion equation (1b) is used and the parameter a is given a value of 0.5 in this paper:(2)mbpeak=mieGSIpeak−100/28−14Dspeak=eGSIpeak−100/9−3D,mbres=mieGSIres−100/28sres=eGSIres−100/9,where the superscripts “peak” and “res” represent the peak parameter and residual parameter, respectively. Hoek and Brown proposed a nonlinear relation between the major and minor principal stresses at failure, assuming independence of the intermediate principal stress [29]. In recent years, many studies on generalized three-dimensional H-B strength criterion have been carried out [38–41]. The commonly used H-B criterion is modified continuously, and the generalized three-dimensional model is the future direction in the research and application.


## 4. Solutions in the Elastic Zone

For the symmetry of a spherical salt cavern subjected to homogeneous internal and external pressures, the corresponding relationships between both the stresses and the strains can be written asσ1=σ2=σθ=σφ,σ3=σr and ε1=ε2=εθ=εφ,ε3=εr. As a result, the salt rock around the cavern satisfies the force equilibrium equation:(3)dσrdr+2rσr−σθ=0.The compressive stress is assumed to be positive, the initial in situ stressp0 is not considered, and the stresses and strains in the elastic zone can be obtained, respectively, as follows:(4)σre=1−R2r3p0+R2r3σR2,σθe=1+12R2r3p0−12R2r3σR2,(5)εre=1Eσre−p0−2μσθe−p0,εθe=1Eσθe−p0−μσre+σθe−2p0.Then, the elastic strains within the elastic zone atr=R2 can be written as follows:(6)εre=δ,εθe=−12δ,where δ=1/2GσR2−p0 and G=E/21+μ.What is more, the displacement caused by the salt cavern’s construction can be derived as follows:(7)uer=14GσR2−p0R23r2.Assuming thatσr=R2=σR2, σR2 represents the radial stress at r=R2. The values of R2 and σR2 are not yet known. When elastic–perfectly plastic model is used, and strength parameters m and s in the H-B yield criterion are constant values, the values of R2 and σR2 can be obtained [24]. However, if strain-softening occurs in the plastic zone and if m and s in the H-B criterion change after peak strength, numerical methods have to be used to get the values of R2 and/or σR2. Besides, the stresses and strains, on the interface between elastic zone and softening zone and on the interface between softening zone and residual zone, have to be used in the solving process; the detail determination method of the elastoplastic interface and its corresponding stress will be analyzed later. Meanwhile, the stresses and displacements at the elasto-softening interface have to be used to obtain the stresses and displacements on the softening zone and the residual zone.


## 5. Stress Analysis in the Plastic Zone

### 5.1. Stress and Radius of the Strain-Softening Zone

Within the near-zone of a salt cavern,R1≤r≤R2, to build the H-B softening constitutive model, it is of critical importance to define a variable as the softening coefficient, which can reflect the linear degradation law of each softening model parameter during the strain-softening process (Figure 2). The coefficient is given as η=εθsp, where εθsp is the tangential plastic strain in the strain-softening zone. The functional relation between ωη and η is established as follows:(8)ωη=ωres−ωpeakη∗η+ωpeak,0≤η<η∗,ωres,η≥η∗,where ωη can be used to describe the H-B parameter mbs or ss; ωpeak and ωres are the corresponding peak parameter and residue parameter, respectively, and can be obtained by substituting GSIpeak and GSIres into (2), respectively; and η∗ represents the residue value of η.Figure 2
Relationships between softening coefficient and H-B parameters.According to the definition of the softening coefficient, we can obtain the following:(9)K⟶∞elastic−brittle−plasticity,η∗=0,K=0ideal elasto−plastic,η∗⟶∞,where, the parameter K is the curve slope for the strain-softening process. When K = 0, the elastic-brittle-plasticity model degenerated to an ideal elastoplastic model [42].After the rock around the cavern enters the strain-softening state,σθsη and σrsη vary with the softening coefficient and conform to the H-B criterion. To determine the stress within the strain-softening zone, the classical H-B criterion can be rewritten as follows:(10)σθsη=σrsη+σcimbsσrsησci+ss1/2.Combining (3) and (10), the radial stress can be derived as follows:(11)1mbsmbsσrsησci+ss1/2=lnr+C.Substitutingσrs=σR2 into (11) at r=R2, we can obtain the radial stress of the strain-softening zone as follows:(12)σrsη=σR2+2A1lnrR2+A2ln2rR2,where A1=mbsσciσR2+ssσci2; A2=mbsσci; mbs and ss vary with the variation of softening coefficient η; and 0≤η≤η∗.According to (10), the tangential stress can be determined as follows:(13)σθsη=σrsη+A1+A2lnrR2.Accounting for the boundary condition atr=R2, we can obtain the following:(14)3σR2+2mbsσciσR2+ssσci2−3p0=0,where the parameters mbs and ss are given the residual values.Because (14) is a nonlinear equation, fσR2=0, the iteration method can be employed to determine the value of σR2 in the strain-softening zone at r=R2. To determine the radius of the strain-softening zone, adding the stress–strain curve in Figure 1(b) and assuming that the elastic strain in the plastic zone abides by Hooke’s law, we can express the stress in the strain-softening zone as follows:(15)σθs−σrs=32p0−σR2−λεθs−σθs−σrsE,where λ=GSIpeak−GSIres/εR1−εR2.Considering the boundary conditions atr=R1 and combining (13) and (15), assuming R1/R2=t, we can obtain the following:(16)A1R1+A2R1lnt−32p0−σR2+λεR1−σθs−σrsE=0,where A1R1 and A2R1 are obtained by mbs and ss at r=R1, respectively.


### 5.2. Stress and Radius in the Plastic Flow Zone

Within the near-zone of a salt cavern,R0≤r≤R1, the stress satisfies (1b) and (3). The H-B parameters mbr and sr are residual values, which can be derived by combining (2) and parameter GSIres. Combining (1b) and (3) results, we obtain the radial stress of the plastic flow zone as(17)1mbrmbrσrrσci+sr1/2=lnr+C.Substitutingσrr=pi into (17) at r=R0 results in the radial stress:(18)σrr=pi+2B1lnrR0+B2ln2rR0,where B1=mbrσcipi+srσci2 and B2=mbrσci.Adding (1b), we can determine the stress of the plastic flow zone, σθr, as follows:(19)σθr=σrr+B1+B2lnrR0.Accounting for the continuous conditions, atr=R1, the stress of the plastic flow zone σrr should equal that of the strain-softening zone σrs. Combining (12) and (18), we can obtain the radial stress as(20)σR1=pi+2B1lnR1R0+B2ln2R1R0=σR2+2A1lnt+A2ln2t.Adding (20), we can determine the range of the plastic flow zone as(21)R1=R0eB12+B2B3−B1/B2,where B3=σR1−pi.


## 5.1. Stress and Radius of the Strain-Softening Zone

Within the near-zone of a salt cavern,R1≤r≤R2, to build the H-B softening constitutive model, it is of critical importance to define a variable as the softening coefficient, which can reflect the linear degradation law of each softening model parameter during the strain-softening process (Figure 2). The coefficient is given as η=εθsp, where εθsp is the tangential plastic strain in the strain-softening zone. The functional relation between ωη and η is established as follows:(8)ωη=ωres−ωpeakη∗η+ωpeak,0≤η<η∗,ωres,η≥η∗,where ωη can be used to describe the H-B parameter mbs or ss; ωpeak and ωres are the corresponding peak parameter and residue parameter, respectively, and can be obtained by substituting GSIpeak and GSIres into (2), respectively; and η∗ represents the residue value of η.Figure 2
Relationships between softening coefficient and H-B parameters.According to the definition of the softening coefficient, we can obtain the following:(9)K⟶∞elastic−brittle−plasticity,η∗=0,K=0ideal elasto−plastic,η∗⟶∞,where, the parameter K is the curve slope for the strain-softening process. When K = 0, the elastic-brittle-plasticity model degenerated to an ideal elastoplastic model [42].After the rock around the cavern enters the strain-softening state,σθsη and σrsη vary with the softening coefficient and conform to the H-B criterion. To determine the stress within the strain-softening zone, the classical H-B criterion can be rewritten as follows:(10)σθsη=σrsη+σcimbsσrsησci+ss1/2.Combining (3) and (10), the radial stress can be derived as follows:(11)1mbsmbsσrsησci+ss1/2=lnr+C.Substitutingσrs=σR2 into (11) at r=R2, we can obtain the radial stress of the strain-softening zone as follows:(12)σrsη=σR2+2A1lnrR2+A2ln2rR2,where A1=mbsσciσR2+ssσci2; A2=mbsσci; mbs and ss vary with the variation of softening coefficient η; and 0≤η≤η∗.According to (10), the tangential stress can be determined as follows:(13)σθsη=σrsη+A1+A2lnrR2.Accounting for the boundary condition atr=R2, we can obtain the following:(14)3σR2+2mbsσciσR2+ssσci2−3p0=0,where the parameters mbs and ss are given the residual values.Because (14) is a nonlinear equation, fσR2=0, the iteration method can be employed to determine the value of σR2 in the strain-softening zone at r=R2. To determine the radius of the strain-softening zone, adding the stress–strain curve in Figure 1(b) and assuming that the elastic strain in the plastic zone abides by Hooke’s law, we can express the stress in the strain-softening zone as follows:(15)σθs−σrs=32p0−σR2−λεθs−σθs−σrsE,where λ=GSIpeak−GSIres/εR1−εR2.Considering the boundary conditions atr=R1 and combining (13) and (15), assuming R1/R2=t, we can obtain the following:(16)A1R1+A2R1lnt−32p0−σR2+λεR1−σθs−σrsE=0,where A1R1 and A2R1 are obtained by mbs and ss at r=R1, respectively.


## 5.2. Stress and Radius in the Plastic Flow Zone

Within the near-zone of a salt cavern,R0≤r≤R1, the stress satisfies (1b) and (3). The H-B parameters mbr and sr are residual values, which can be derived by combining (2) and parameter GSIres. Combining (1b) and (3) results, we obtain the radial stress of the plastic flow zone as(17)1mbrmbrσrrσci+sr1/2=lnr+C.Substitutingσrr=pi into (17) at r=R0 results in the radial stress:(18)σrr=pi+2B1lnrR0+B2ln2rR0,where B1=mbrσcipi+srσci2 and B2=mbrσci.Adding (1b), we can determine the stress of the plastic flow zone, σθr, as follows:(19)σθr=σrr+B1+B2lnrR0.Accounting for the continuous conditions, atr=R1, the stress of the plastic flow zone σrr should equal that of the strain-softening zone σrs. Combining (12) and (18), we can obtain the radial stress as(20)σR1=pi+2B1lnR1R0+B2ln2R1R0=σR2+2A1lnt+A2ln2t.Adding (20), we can determine the range of the plastic flow zone as(21)R1=R0eB12+B2B3−B1/B2,where B3=σR1−pi.


## 6. Deformation Analysis in the Plastic Zone

According to the plastic theory, the strain in the plastic zone can be divided into the elastic and plastic strains. For salt rock, the large strain theory is more suitable for the plastic strain. Generally when the degree of deformation is less than 10%, it can be expressed by engineering strain; otherwise, logarithmic strain should be used. The deformation of salt rock under triaxial stress state indicates the large deformation characteristics of salt rock under high confining pressure. Using the logarithmic strain to analyze the large deformation characteristics of salt rock is more suitable [30]. Therefore, the strain can be expressed in the logarithmic form as ε=∫r0rdl/l=lnr/r0, where r0 and r are variables in the calculations. Then, the geometric equation can be rewritten as (22). Although the elastic strain in the plastic zone has various definitions [43], it is assumed that it conforms to Hooke’s law in this study:(22)εr=lndrdr0,εθ=εφ=lnrr0.


### 6.1. Displacement Analysis in the Strain-Softening Zone

In the zone ofR1≤r≤R2, the strains can be stated as follows:(23)εrs=εrse+εrsp,εθs=εθse+εθsp,where the superscripts “e” and “p” denote the elastic and plastic strains of the strain-softening zone, respectively.Employing the flow rule establishes the relationship between the components of the plastic strain. According to the nonassociated flow law, we can get the relationship between the plastic strain increments as follows:(24)dεrspdεθsp=−2h,where h is a function of ψ, h=1−sinψ/1+sinψ, and ψ is the angle of dilation obtained by a triaxial compressive test.Combining (23) and (24) results in the following:(25)hdεrs+2dεθs=hdεrse+2dεθse.Based on the integration of (25) and by adding (5), we can derive the following:(26)hεrs+2εθs=fsr,where fsr can be given as fsr=1/2G1+μh−2μσrs+21−μ−hμσθs−h+21−2μp0.After substituting the radial stress of (12) and tangential stress of (13) into the above form for fsr, we can obtain(27)fsr=12G1+μC1+C2lnrR2+C3ln2rR2,where C1=h+21−2μσR2−p0+21−μ−hμA1, C2=2h+21−2μA1+21−μ−hμA2, and C3=h+21−2μA2.Combining (22) and (26) leads to the following equations:(28)hlndrdr0+2lnrr0=fsr,lnrr02/hdrdr0=fsrh.Definite integral calculation can be performed on the intervalsr0∈r0,R2−uR2 and r∈r,R2, where uR2 represents the displacement around a salt cavern in the elastic zone at r=R2. Based on (7), we can obtain uR2=1/2δR2.The integration of (28) leads to the displacement for a large deformation in the strain-softening zone:(29)usr=r−r0=r−1−12δR2h+2/h−h+2h∫rR2r2/he−fsr/hdrh/h+2.Using the method of definite integration by substitution, we can deduce (29) as follows:(30)usr=r−r0=r−1−12δR2h+2/h−h+2h∫01R2−rx+r2/hR2−re−1/2G1+μC1+C2lnR2−rx+r/R2+C3ln2R2−rx+r/R2/hdxh/h+2.Furthermore, the solution of (30) can be derived based on Romberg integration. At r=R1, according to (28), the tangential strain at the interface of the strain-softening zone and plastic flow zone can be given as follows:(31)εR1=lnrr0=ln1−12δ1th+2/h−h+2hT⇒T=hh+21−12δ1th+2/h−eεR1,where T=∫011/t−1x+12/h1/t−1e1/hEC1+C2ln1−tx+t+C3ln21−tx+tdx.Combining (16) and (31), we can obtain the nonlinear differential equation for t. However, this equation has no closed-form solution. To find the values of t and εR1, the numerical iteration method is employed with the aid of MATLAB. Next, substituting lnt into (20) produces σR1. Then, substituting σR1 into (21), we can obtain R1 and R2.


### 6.2. Displacement Analysis in the Plastic Flow Zone

Likewise, in the zone ofR0≤r≤R1, the strains can be expressed as follows:(32)εrr=εrsr=R1+εrrp,εθr=εθsr=R1+εθrp.Based on the nonassociated flow law, we can determine the plastic strain increments as follows:(33)dεrrpdεθrp=−2f,where f is a function of ψ′; f=1−sinψ′/1+sinψ′; and ψ′ can be got through a triaxial compressive test.Combining (32) and (33), we will have the following:(34)fdεrr+2dεθr=fdεrre+2dεθre.Based on the integration of (34) and by adding (5), we can derive the following:(35)fεrr+2εθr=frr,where frr=1/2G1+μσrr−p0−2μσθr−p0f−h+fsR1.After substituting the radial stress of (18) and tangential stress of (19) into frr, we can obtain(36)frr=f−h2G1+μD1+D2lnrR0+D3ln2rR0+fsR1,where D1=1−2μpi−p0−2μB1; D2=21−2μB1−2μB2; D3=1−2μB2; and fsR1=1/2G1+μC1+C2lnR1/R2+C3ln2R1/R2.Combining (32) and (34) results in the following:(37)flndrdr0+2lnrr0=frr,lnrr02/fdrdr0=frrf.Definite integral calculation can be performed onr0∈r0,R1−uR1 and r∈r,R1, where uR1 is based on (30) and represents the displacement in the strain-softening zone at r=R1. The integration of (37) leads to the displacement for a large deformation in the plastic flow zone:(38)urr=r−r0=r−R1−uR1f+2/f−f+2f∫rR1r2/fe−frr/fdrf/f+2.Using the definite integration by substitution, we can deduce (38) as follows:(39)urr=r−r0=r−R1−uR1f+2/f−f+2f∫01R1−rx+r2/fR1−re−f−h/fED1+D2lnR1−rx+r/R0+D3ln2R1−rx+r/R0−1/ffsR1dxf/f+2.Then, (39) can be solved by the method of Romberg integration.


## 6.1. Displacement Analysis in the Strain-Softening Zone

In the zone ofR1≤r≤R2, the strains can be stated as follows:(23)εrs=εrse+εrsp,εθs=εθse+εθsp,where the superscripts “e” and “p” denote the elastic and plastic strains of the strain-softening zone, respectively.Employing the flow rule establishes the relationship between the components of the plastic strain. According to the nonassociated flow law, we can get the relationship between the plastic strain increments as follows:(24)dεrspdεθsp=−2h,where h is a function of ψ, h=1−sinψ/1+sinψ, and ψ is the angle of dilation obtained by a triaxial compressive test.Combining (23) and (24) results in the following:(25)hdεrs+2dεθs=hdεrse+2dεθse.Based on the integration of (25) and by adding (5), we can derive the following:(26)hεrs+2εθs=fsr,where fsr can be given as fsr=1/2G1+μh−2μσrs+21−μ−hμσθs−h+21−2μp0.After substituting the radial stress of (12) and tangential stress of (13) into the above form for fsr, we can obtain(27)fsr=12G1+μC1+C2lnrR2+C3ln2rR2,where C1=h+21−2μσR2−p0+21−μ−hμA1, C2=2h+21−2μA1+21−μ−hμA2, and C3=h+21−2μA2.Combining (22) and (26) leads to the following equations:(28)hlndrdr0+2lnrr0=fsr,lnrr02/hdrdr0=fsrh.Definite integral calculation can be performed on the intervalsr0∈r0,R2−uR2 and r∈r,R2, where uR2 represents the displacement around a salt cavern in the elastic zone at r=R2. Based on (7), we can obtain uR2=1/2δR2.The integration of (28) leads to the displacement for a large deformation in the strain-softening zone:(29)usr=r−r0=r−1−12δR2h+2/h−h+2h∫rR2r2/he−fsr/hdrh/h+2.Using the method of definite integration by substitution, we can deduce (29) as follows:(30)usr=r−r0=r−1−12δR2h+2/h−h+2h∫01R2−rx+r2/hR2−re−1/2G1+μC1+C2lnR2−rx+r/R2+C3ln2R2−rx+r/R2/hdxh/h+2.Furthermore, the solution of (30) can be derived based on Romberg integration. At r=R1, according to (28), the tangential strain at the interface of the strain-softening zone and plastic flow zone can be given as follows:(31)εR1=lnrr0=ln1−12δ1th+2/h−h+2hT⇒T=hh+21−12δ1th+2/h−eεR1,where T=∫011/t−1x+12/h1/t−1e1/hEC1+C2ln1−tx+t+C3ln21−tx+tdx.Combining (16) and (31), we can obtain the nonlinear differential equation for t. However, this equation has no closed-form solution. To find the values of t and εR1, the numerical iteration method is employed with the aid of MATLAB. Next, substituting lnt into (20) produces σR1. Then, substituting σR1 into (21), we can obtain R1 and R2.


## 6.2. Displacement Analysis in the Plastic Flow Zone

Likewise, in the zone ofR0≤r≤R1, the strains can be expressed as follows:(32)εrr=εrsr=R1+εrrp,εθr=εθsr=R1+εθrp.Based on the nonassociated flow law, we can determine the plastic strain increments as follows:(33)dεrrpdεθrp=−2f,where f is a function of ψ′; f=1−sinψ′/1+sinψ′; and ψ′ can be got through a triaxial compressive test.Combining (32) and (33), we will have the following:(34)fdεrr+2dεθr=fdεrre+2dεθre.Based on the integration of (34) and by adding (5), we can derive the following:(35)fεrr+2εθr=frr,where frr=1/2G1+μσrr−p0−2μσθr−p0f−h+fsR1.After substituting the radial stress of (18) and tangential stress of (19) into frr, we can obtain(36)frr=f−h2G1+μD1+D2lnrR0+D3ln2rR0+fsR1,where D1=1−2μpi−p0−2μB1; D2=21−2μB1−2μB2; D3=1−2μB2; and fsR1=1/2G1+μC1+C2lnR1/R2+C3ln2R1/R2.Combining (32) and (34) results in the following:(37)flndrdr0+2lnrr0=frr,lnrr02/fdrdr0=frrf.Definite integral calculation can be performed onr0∈r0,R1−uR1 and r∈r,R1, where uR1 is based on (30) and represents the displacement in the strain-softening zone at r=R1. The integration of (37) leads to the displacement for a large deformation in the plastic flow zone:(38)urr=r−r0=r−R1−uR1f+2/f−f+2f∫rR1r2/fe−frr/fdrf/f+2.Using the definite integration by substitution, we can deduce (38) as follows:(39)urr=r−r0=r−R1−uR1f+2/f−f+2f∫01R1−rx+r2/fR1−re−f−h/fED1+D2lnR1−rx+r/R0+D3ln2R1−rx+r/R0−1/ffsR1dxf/f+2.Then, (39) can be solved by the method of Romberg integration.


## 7. Discussion

Whenη∗⟶∞, there is no softening behavior, and the zoned phenomenon does not remain in the plastic zone. In addition, the constitutive relationship of the salt rock is the ideal elastoplastic one. The radial stress and tangential stress in the plastic zone can be written as (40) and (41), respectively:(40)σrp=pi+2A1lnrR0+A2ln2rR0,where A1=mbσcipi+sσci2, A2=mbσci, and mb and s are stationary values obtained by substituting the parameter GSIpeak into (2):(41)σθp=σrp+A1+A2lnrR0.Considering the boundary condition atr=R2, we can express (14) as follows:(42)3σR2+2mbσciσR2+sσci2−3p0=0.The radial stress atr=R2, σR2, can be derived through numerical iteration. Then, plastic radius R2 can be determined according to (40):(43)R2=R0eA12+A2σR2−pi−A1/A2.When calculating the displacement in the plastic zone, (26) is rewritten as follows:(44)hεrp+2εθp=fpr,where fpr=1/2G1+μK1+K2lnr/R0+K3ln2r/R0, K1=h+21−2μpi−p0+21−μ−hμA1, K2=2h+21−2μA1+21−μ−hμA2, and K3=h+21−2μA2.Combining (22) and (44), we can obtain the following:(45)lnrr02/hdrdr0=fprh.Moreover, definite integral calculation can be performed onr0∈r0,R2−uR2 and r∈r,R2, where uR2 is based on (7) and represents the displacement in the elastic zone at the interface of the plastic zone and elastic zone. The integration of (45) leads to the displacement for a large deformation in the plastic zone with an ideal elastoplastic rock mass:(46)upr=r−r0=r−R2−uR2h+2/h−h+2h∫rR2e−fpr/hr2/hdrh/h+2.Using the method of the definite integration by substitution, we can deduce (46) as follows:(47)upr=r−r0=r−R2−uR2h+2/h−h+2h∫01e−1/2G1+μK1+K2lnR2−rx+r/R0+K3ln2R2−rx+r/R0/hR2−rx+r2/hR2−rdxh/h+2.Furthermore, the solution of (47) can be determined based on Romberg integration.


## 8. Case Study

Large-scale underground salt cavern gas storage construction has been carried out in Huaian, Jiangsu, Pingdingshan, Henan, Wanzhou, Chongqing, Yunying, Hubei, etc., those storage groups will be used as oil and gas reserves. At present, there is just one underground gas storage salt cavern (Jintan gas storage) in China at running condition [44]. To show the accuracy of the proposed solutions, the results of the displacement, the radial and circumferential stresses, and the variation of r/R0 are compared using the dataset as shown in Table 1. Table 1 lists the salt cavity structural and rock mechanical parameters. E, μ, and ψ denote the elastic modulus, poison ratio, and dilation angle of salt rock, respectively. Table 2 lists the H-B parameters of the salt rock. After the water encroachment of salt cavity storage, the damaged quality and structure of the rock mass decrease the GSI value. Accounting for the empirical regression H-B parameter, GSI affects the displacement, with different GSI values of 60 and 70.Table 1
Parameters.
R0 (m)p0 (MPa)pi (MPa)E (GPa)μψ304010100.30,φ/4,φ/2Table 2
H-B parameters of salt rock.
GSISGSIrσci (MPa)miDmbsssmbrsr603424.44.10.20.840.00850.300.0004703624.44.10.21.250.02810.320.0005Primarily, through numerical iteration,σR2 is calculated on the basis of (14). Furthermore, the nonlinear differential equations (48) are utilized according to (16) and (31). Then, the unknown variables t and εR1 are determined with the aid of the MATLAB software. Next, σR1,R1, and R2 are calculated one by one. Table 3 lists the related data:(48)A1R1+A2R1lnt−32p0−σR2+λεR1−σθs−σrsE=0,∫01θtdx=hh+21−12δ1th+2/h−eεR1,where θt=1/t−1x+12/h1/t−1e1/hEC1+C2ln1−tx+t+C3ln21−tx+t.Table 3
Dataset.
GSIδσR2εR1tσR1R1R2A1A2B1B2B360−0.00224.8840.00110.642714.197137.514958.3708A1 (0) = 22.674A1 (η∗) = 13.5051A2 (0) = 20.496A2 (η∗) = 7.328.56967.324.197170−0.002322.390.00120.691313.154035.416551.23A1 (0) = 27.815A1 (η∗) = 13.950A2 (0) = 30.420A2 (η∗) = 7.8088.8537.8083.1540First, it is particularly necessary to validate the analytical result. Therefore, the elastic-brittle-plasticity model is simplified to an ideal elastoplastic model. The solution of (47) is derived and compared with the numerical solutions that are carried out with the use of the H-B constitutive model in FLAC3D. Figure 3 shows the principle stress and displacement of the rock masses around the cavity, respectively. It can be seen from Figure 3 that the solutions obtained with different H-B criterion parameters are in accord with the distribution law of deep spherical adits and agree with the published outcomes [24]. The variation in GSI affects the principle stress but not obviously. Without considering softening behavior, the analytical solution is closer to the numerical solution. However, when taking the softening behavior into account, the plastic zone is enlarged, and a great difference exists between the analytical and numerical results, particularly at the cavity wall. The numerical solution for the tangential stress is significantly larger than the analytical solution. It can be seen from Figure 3(b) that GSI markedly affects both the displacement at the cavity wall and the plastic radius of the surrounding rock. For example, the GSI increases by 10 and, as a result, the displacement at the cavity wall decreases from 0.72 m to 0.51 m, while the plastic radii R1 and R2 decrease by 6% and 12%, respectively. The influence of GSI on R2 is more obvious. The numerical result for the displacement obtained when adopting the ideal elastoplastic model is evidently smaller than the analytical result obtained when considering the softening behavior, but it should be noted that the numerical result is fundamentally in accord with the analytical result obtained without taking the softening behavior into account. It can be seen from the results that the variation of GSI and/or the softening behavior greatly affects the results. Thus, if the effects of the two factors are not considered, the conclusions will have large deviations.Figure 3
Validation of H-B solution against numerical results for spherical domain: (a) distributions of principal stresses; (b) distributions of displacements.
(a)(b)Alejano and Alonso [45] suggested a definition of dilatancy; it is a volume change caused by the shear distortion which could not recover. Hoek and Brown [46] proposed that the dilation angle of rocks may be extremely small. A dilation angle of zero shows that the volume is preserved. Three cases, ψ=0, ψ=0.25φ, and ψ=0.5φ, are considered here. The relationships between the dilation angles and the displacements within the plastic zone of a salt cavern are illustrated in Figure 4. The figure shows that the displacement at the same point increases with the dilation angle, but the degree of this effect decreases with an increase in the distance away from the salt cavern. The influence of the dilation angle on the displacement at the cavity wall is the most striking. Taking the softening behavior into account, when the angle of dilation is ψ=0, ψ=0.25φ, or ψ=0.5φ, the displacement at the cavity wall is 0.72 m, 0.84 m, or 0.97 m, respectively. Thus, the extent of the influence increases with the dilation angle. We next consider the displacements of the plastic flow zone for cases of f=0.7699, f=0.8687, and f=1 when h = 0.7699 (ψ=0.25φ) and f≠h. Figure 5 illustrates that the displacement of the plastic flow zone decreases with an increase in f. In other words, the displacement increases with dilation angle ψ′, but the variations are minor. Compared to h, the effect of f on the displacement is relatively small.Figure 4
Relationships between dilation angle and displacements of salt cavern in plastic zone.Figure 5
Displacements of surrounding rock in strain-softening zone.Figure6 indicates the relationships between the operating pressure and displacement at the cavity wall. As we can see from Figure 6, the displacement decreases with an increase in the operating pressure. When the operating pressure is increased six times, both the analytical results and numerical results are decreased by 80%, and the reduction rates are essentially the same. Under a pressure of 2 MPa, an analytical solution of 2.4 m is attained, but when the softening behavior is not considered, the analytical solution and numerical solution are merely 0.55 m and 0.75 m, respectively. Thus, the difference is evident. The variations are caused by the fact that the numerical solution does not take the softening behavior into account, and the definition of the elastic strain within the plastic zone is different. The results demonstrate that the operating pressure has a significant effect on the displacement of the salt cavern storage at the cavity wall. In order to ensure safe operation of the reservoir, the minimum running pressure should not be too low. Thus, it is safer to design the operating pressure based on the analytical solution.Figure 6
Relationships between displacement activity wall and operating pressure.


## 9. Conclusions

(1)
Under the condition of three-dimensional model, considering the large strain characteristics of salt rock, the logarithmic strain is introduced to analyze the theoretical solution for the large deformation of the surrounding rock of salt cavity. And a semianalytical solution is obtained via Romberg numerical method. The stress, strain, and large deformation displacement of the strain-softening and plastic flow zones of a salt cavern in H-B rock masses were obtained.(2)
The H-B parameter, GSI, was found to have a significant impact on the radius and displacement of the plastic zone. Both the plastic radius and displacement decreased with an increase in the GSI parameter, which indicated that the plastic zone and displacement decreased with an increase in the rock mass.(3)
The dilation angle also affected the distribution of the displacement in the salt cavern. The displacement was positively correlated with the angle of dilation. Compared with the impact of the strain-softening parameter,h, the flow parameter, f, had little effect on the displacement. The dilation angle and GSI affected the displacement of the salt cavern most significantly at the cavity wall.(4)
The displacement and plastic zone became obviously greater when the strain-softening behavior was considered, especially for a poor rock mass and lower operating pressure. Therefore, the strain-softening properties cannot be ignored when designing an actual project.



---

*Source: 8851177-2020-08-31.xml*