Stress redistribution in a multilayer chamber for compressed air energy storage in abandoned coalmine: Elastic analytical insights and material choice

Compressed air energy storage (CAES) is attracting attention as one of large‐scale renewable energy storage systems. Its gas storage chamber is one of key components for its success. A successful utilization of an abandoned coalmine roadway depends on the stability of the gas storage chamber. The chamber is a multilayer structure and the redistribution of the stress and displacement in each layer is critical to the chamber stability. So far, this redistribution mechanism under any lateral pressure coefficient and internal air pressure has been unclear and thus a quick and easy evaluation on the CAES chamber stability becomes difficult. In this study, the redistributions of stress and displacement in each layer are analytically solved based on complex variable function theory. First, the stress and displacement in three CAES multilayer structures are obtained under any lateral pressure coefficient and internal air pressure. Then, a comprehensive parameter b1 is obtained to describe the redistribution of stress and displacement in rock mass, lining layer, and grout layer. Its linkage with lateral pressure coefficient, internal air pressure, and material properties is analytically expressed. Thirdly, an analytical expression is obtained for the influence range of internal air pressure on chamber stress. Finally, the role and selection of lining and grouting layers are explored at different lateral pressure coefficient and internal air pressure. It is found that the CAES chamber stability can be effectively and quickly evaluated by these analytical solutions and comprehensive parameters and enhanced by a material with low bulk modulus but high shear modulus.


| INTRODUCTION
The widespread use of renewable clean energy (such as hydropower, solar energy, and wind energy) requires a large-scale energy storage system to regulate the mismatch between energy demand and supply.Compressed air energy storage (CAES) technology as an emerging large-scale energy storage can solve the temporal and spatial mismatch in grid peak and energy use. 1,2The concept of using underground chamber as CAES was proposed by Stal Laval in 1949 3 and China now has the potential to develop large-scale and highquantity underground gas storage facilities. 4Until now, the dominant chamber options have been underground salt caverns, 5,6 abandoned mine chambers 7,8 or gas storage chambers in hard rock formations. 9,10The success of a CAES lies in successfully addressing the following issues: chamber stability, 11 compressed air leakage, 12 environmental impacts, 9 and some system factors. 9,13In terms of technical issues, the stability of gas storage chamber is caused by stress and displacement. 14Therefore, this paper will explore the issue of chamber stability through the analysis of stress and displacement.
The CAES chamber stability is a critical problem to the construction of CAES chamber with abandoned coalmine roadways.The roadways usually have soft surrounding rocks with excavation-induced damages, 14 thus a CAES chamber in coal mines has three features: First, the chamber is usually a three-layer structure with surrounding rock, lining reinforced concrete and grout layer.Second, it is usually subjected to nonhydrostatic external loads.Third, it is subjected to an internal air pressure due to the compressed air storage.The stress and displacement in the chamber are crucial to the chamber stability.Analytical solutions of stress and displacement in the storage chamber can quickly evaluate the chamber stability and thus favor to feasibility design.
The distribution of stress and displacement around the chamber can be obtained by either analytical solutions or numerical simulations.Finite element method (FEM) is the most widely used method for numerical simulations.Kim et al. 15 showed that tight concrete lining and enclosure structures with small damage areas are favor to chamber stability.Perazzelli et al. 16,17 proposed a limit equilibrium method for cavern stability analysis, analyzed the influence of geometric and geotechnical parameters on the uplift damage pressure, and explored the feasibility of CAES construction under different rock conditions.Prado et al. 18 analyzed the thermodynamic properties of a chamber with a seal layer and a 35 cm concrete lining when the working air pressure of chamber was between 5 and 8 MPa.Menéndez et al. 19 considered a 15 cm thickness of lining and a 5 m thickness of rock and analyzed the variation of gas temperature and pressure within the chamber under different operation conditions.The above results obtained through numerical simulations can evaluate the stability of a CAES chamber.However, an analytical solution is still being explored for a multilayer structure.A quick estimation of the distribution, regulation, and allocation of stress and displacement in composite structure chamber is extremely important to the stability evaluation of a CAES cavern, particularly in the stage of feasibility design.Some analytical solutions have been developed for the distribution of stress and displacement in a singlelayer chamber.For a circular chamber with a single layer, Kirsch 20 developed an analytical solution based on elasticity.Verruijt 21 applied the Laurent series to the complex function expansion and obtained the distribution of stress.Strack and Verruijt 22 further obtained the solution of the chamber stress under the gravity load in a half-plane elastomer.Lu et al. 23 obtained the stress solution of a single chamber under general external loads.Their solution took the effects of gravity and an arbitrary lateral stress into account, but they did not consider the presence of internal loads.Fahimifar and Zareifard 24 analyzed the tunnel stress distribution below water table under axisymmetric conditions, but they did not obtain a closed-form solution for this problem.Zhang et al. 25 used the unified strength theory to provide an analytical solution of the stress in a singlelayer structure, but the boundary conditions in their calculation model are symmetric and too simple.For single-layer structures, most of the analytical solutions are obtained at symmetric boundaries but without internal loads.
For a two-layer chamber, Wang and Li 26 obtained an analytical solution of stress through the liner-rock interaction with complex variable function.Lin et al. 27 used a complex variable approach and discrete Fourier transform theory to investigate the shallow buried chambers supported by segment linings.Zhou et al. 28 developed an axisymmetric model for multilayer chambers under internal loads and calculated the hydraulic-mechanical interaction in hydraulic chambers using the equivalent coupling method.These analytical solutions considered internal loads, but their structure has only two layers of rock mass and lining and their external loads were axially symmetrical.Thus, these solutions confront some difficulty in their application to the case under any lateral pressure coefficient.
A three-layer structure without internal gas pressure was analytically solved by EI et al. 29 with elasticity.They found that internal air pressure cannot be treated as a simple boundary traction.Mason and Abelman 30 also considered a three-layer structure, but only considered external boundary loads and ignored internal load.The analytical solutions of stress and displacement for three-or more-layer chambers under any lateral pressure coefficient and internal air pressure are still in development.The target is to understand the load transfer mechanism in different layers of the threelayer chamber under complex loads.
The multilayer chamber problems under complex loads can be solved in complex elasticity. 31Complex elasticity transforms an analysis problem through complex functions and uses complex variables to accurately describe the mechanical responses of the multilayer structures.Classic elasticity requires an inverse or semiinverse analysis for solutions, and the inverse process is complex. 32However, complex elasticity only requires analytical functions to conform to the Cauchy-Riemann equation.This greatly simplifies the calculation with analytical expression.In this paper, the complex elasticity is used to develop analytical solutions for multilayer chambers under both internal air pressure and nonhydrostatic external loads.Compared to previous simplified analytical solutions, this study tries to make three main improvements.First, the new analytical solutions can not only calculate the mechanical response of single-layer chambers, but also calculate multilayer chambers (up to three layers).The analytical solutions can be quickly calculated through our self-developed MATLAB code program.Second, the complex variable functions are combined with elasticity to avoid the difficulties brought by inverse and semi-inverse methods.Finally, the complex loads are included in the chamber and thus their impacts on each layer are clearly expressed.Particularly, the distributions of stress and displacement are quantitatively analyzed through analytical methods and a comprehensive parameter (b 1 ) is obtained for the redistribution of stress and displacement in rock mass, lining layer and grout layer under any lateral pressure coefficient, internal air pressure, and material properties.
The paper is organized as follows: In Section 2, the analytical solutions for multilayer chambers are derived using complex variable function and elasticity.The selfconsistency of these analytical solutions is checked through solution downgrade.In Section 3, these analytical solutions are verified with the analytical results in the literature and our FEM simulations.In Section 4, the difference between these analytical solutions and FEM solutions in the vertical and horizontal directions is further verified by stress concentration factors.The choice of different layer materials for the multilayer chamber is discussed.The influence range of internal air pressure on the chamber stress is analytically obtained.This demonstrates the application of these analytical solutions as theoretical support to CAES construction practice.Conclusions are drawn in Section 5.

| Analytic functions and general solution of stress and displacement
The distribution of stress and displacement in a chamber depends on the initial stress, the shape and location of the chamber, and the mechanical properties of rock mass. 33Complex functions can be used to solve this chamber problem if they satisfy both deformation compatibility and boundary conditions. 28If z = re iθ , the stress function should satisfy the equation as 34 where φ(z) and ψ(z) are complex analytic functions, and Re indicates the real part of these functions.The stress and displacement are obtained in the complex plane as where G is the shear modulus, μ is the Poisson's ratio, and k = (3 − 4μ) for the plane strain condition.The σ r , σ θ , τ rθ , u, and v are the radial stress, circumferential stress, shear stress, radial displacement, and circumferential displacement, respectively.It is noted that the methods described in Equations ( 1) and ( 2) are applicable to all geometric shapes. 35However, conformal transformations may be applied to different geometric shapes for simpler analytical functions φ(z) and ψ(z). 36ased on the Saint-Venant principle, 37 the selfweight of rock mass can be equivalent to a surface vertical load on boundary, while the horizontal load is generated according to the lateral pressure coefficient for different geological conditions.In this sense, Equation ( 2) is still applicable in dealing with deeply buried chambers.This simplification-induced error in the circumferential stress and displacement of the hole is insignificant. 37The most critical issue to the solution of a plane strain problem is to find two analytic functions φ(z) and ψ(z).They should satisfy both deformation compatibility and boundary conditions.The analytical solutions for three individual composite structures are obtained based on the following basic assumptions: (1) The material of each individual layer is of continuity, uniformity, isotropy, and elasticity.(2) The deformation of the composite structure is small.
(3) The boundaries of each layer are completely contacted in the initial and deformation processes.(4) A plane strain problem for the multilayer structure is considered.(5) Both external and internal surface loads are uniformly applied to the computational domain.

| A single-layer chamber
The single-layer chamber is shown in Figure 1.The analytic functions for the stress and displacement in rock mass are taken as where γ is the specific weight of rock, H is the vertical depth of circular chamber center from the ground, r 1 is the radius of circular chamber and K 0 is the lateral pressure coefficient.a 1 , b 1 , b 2 are coefficients to be determined.
After the analytic function of Equation ( 3) is substituted into Equation (2) and the real and imaginary parts are separated, the stress and displacement of rock mass are obtained as In the subscript labels here, the first letter stands for the material layer.For example, the first letter in σ rρ is r which stands for rock and the following letter is for the polar coordinates.u rρ stands for the radial displacement, and v rθ for the circumferential displacement.G r is the shear modulus of rock, μ r is the Poisson's ratio, and k r = (3 − 4μ r ).
The coefficients in Equation ( 4) are presented below, and their derivation process is given in Appendix I (Part 1).
where q a is the air pressure inside the chamber (called internal air pressure later).
The analytical solution of the stress and displacement in a single-layer chamber is obtained through Equation ( 4).The stress has two components: one is induced by the initial geostress (K 0 γH and γH) and the other is induced by external loads such as circumferential excavation unloading and internal air pressure.When the internal air pressure was 0 MPa, Kirsch 20 obtained the analytical solution of stress and displacement as This is the special case of Equation (4) without internal air pressure or q a = 0.

| A two-layer chamber with lining
A lining layer is usually added to a chamber for chamber stabilization. 38This practice effectively distributes the stress around chamber and relieves ground settlement, thus being widely applied to engineering practice.In this study, the lining is assumed to be completely contacted with the surrounding rock.That is, the outer radius of lining is geometrically equal to the inner radius of surrounding rock.r 1 is the inner radius of lining, r 2 is the outer radius of lining and is also the inner radius of surrounding rock.Figure 2 presents a two-layer chamber and its loads.This two-layer chamber has two loads: the geostress as the extend load at its external boundary and the air pressure at its internal boundary (cavern boundary).The analytic functions in the lining are The analytic function for the stress in surrounding rock is similar to that in Equation (3).The only replacement is r 1 with r 2 .The subscript l represents the lining and r represents the rock.
Equations ( 3) and ( 7) are used to develop the analytical solution of stress and displacement as 3 are nine coefficients to be determined by boundary conditions.Their final results are presented in Equation ( 9) and their derivation process is shown in Appendix I (Part 2).
F I G U R E 2 A two-layer chamber and its loads. where and I r r = / 1 2 1 .If G l , k l , and r 2 are all zero, this means that the lining layer is not present and the structure becomes a single-layer chamber, Equation ( 9) is degraded to Equation (5).When the internal air pressure is zero, the Kirsch solution in Equation ( 6) is then obtained.Therefore, the solution for two-layer chambers theoretically contains that for single-layer chambers.

| A three-layer chamber with lining and grout layers
A three-layer tunnel is a composite structure composed of a lining layer, a grouting layer, and surrounding rock mass.The surrounding rock mass usually has a lot of defects.These defects include soft rock, voids or numerous natural cracks, and excavation-induced damage.They are not favorable for chamber stability and thus a grouting reinforcement is often required and a grout layer is formed. 39,40Grouting technology is widely used in chamber construction to reinforce the defects rock mass near the lining layer.r 1 is the inner radius of lining, r 2 is the outer radius of the lining but also the inner radius of grout layer, and r 3 is the outer radius of grout layer but also the inner radius of the surrounding rock.Figure 3 is a typical threelayer chamber and its loads.
The analytic function in Equation ( 3) remains for the surrounding rock layer and in Equation ( 7) for the lining layer.The analytic function for the grout layer is taken as The stress and displacement in the grout layer are then obtained as F I G U R E 3 A three-layer chamber and its loads.(12a) The subscript g denotes the grout layer.a b , ,  (12b) (12c)  (12d)   | 4205 where   .
If G g , k g , and r 3 are all zero, the grout layer is not present and the three-layer structure becomes a twolayer chamber with only lining layer and rock mass, Equation ( 12) is degraded to Equation (9).If G g , k g , G l , k l , r 3 , and r 2 are all zero, both the lining layer and grout layer are not present, and the structure becomes a single-layer chamber, Equation ( 12) is degraded to Equation (5).When the internal load is zero, the Kirsch solution in Equation ( 6) is still obtained.These show that the solutions for the three-layer chamber are self-consistent with those for the one-layer and two-layer chambers.

| Comparison and verification of analytical solution
The analytical solution derived in Section 2 is compared with Dong et al. 41 and Du et al. 42 in Figures 4 and 5, respectively.Their solutions also included analytical solutions and numerical simulations.Therefore, we use our analytical solutions as a standard (denoted as i = 3) to evaluate the relative error with their analytical solutions, their numerical simulations, and our simulations in this article.This approach can effectively analyze the relative differences between the analytical solution proposed in this article and the results calculated by other methods.The relative error is calculated by where v i is the calculation result with the ith method and a denotes the analytical solution.i = 1 denotes the calculation based on the analytical solution obtained by others; i = 2 denotes the calculation based on the simulation results obtained others; i = 4 denotes the calculation based on our simulation results.Generally, our analytical solution is in good agreement with their results.Hence, our analytical solution can correctly describe both radial and circumferential stresses.

| Establishment of the numerical model
To verify these analytical solutions, a sketch model is presented in Figure 6 for numerical simulations.The FEM solution is obtained by COMSOL Multiphysics.In this section, the analytical solutions are compared with the FEM results in cloud diagrams to check the rationality of analytical solutions in calculating stress distribution.In the numerical simulations, r r = 2 m and N = 20 m.The lateral pressure coefficient (K 0 ) is taken as 0.5, 1.0, and 1.5, respectively.The lateral pressure coefficient is taken based on different geological conditions.The measurements from Chinese coal mines show that the lateral pressure coefficient is generally distributed between 0.5 and 2.0. 43When the lateral pressure coefficient K 0 = 0.5, the vertical load is twice the horizontal load.When K 0 = 1.0, the vertical load is equal to the horizontal load, and the in situ stress is hydrostatic.When K 0 = 1.5, the horizontal load is 1.5 times the numerical load.The crust stress is stronger.Therefore, these three values of lateral pressure coefficient correspond to three geologoical conditions.The mechanical responses of the chamber under these three different geological conditions should be considered, respectively.The air pressure is taken as 0 and 6 MPa, respectively.No gravity is considered in the domain but vertical and horizontal loads are applied at the external boundary.The variation of stress out of the 20 m range is not significant, thus the lateral pressure is uniformly distributed.

| Validation of analytical solutions for single-layer chamber
The single-layer chamber has following parameters: a circular cavern surface with a radius r = 2 m, the elastic modulus of 13.7 GPa and the Poisson's ratio of 0.24.The insitu rock cores were obtained by using the hollow inclusion coring method.The core samples were subsequently tested in laboratory and the stress-strain curve was obtained and presented in Figure 7.The elastic modulus is obtained from the slope of the linear stress-strain curve, and the Poisson's ratio is obtained from the ratio of radial strain to axial strain.The far-field vertical stress is 10 MPa and the horizontal stress is determined by lateral pressure coefficient.Von Mises stress expresses the complex stress state.This simplification can more directly evaluate the current stress state and quickly compare the analytical solution with finite element results.The cloud results for only a lateral pressure coefficient of 1.5 are presented here.Figure 8 shows the stress cloud for the single-layer chamber.The symbol has the following meanings.R denotes the rock layer, A denotes the analytical solution, and S denotes the numerical solution.The final number indicates the air pressure, where 0 means no air pressure and 1 means air pressure.The results for other two lateral pressure coefficients are shown in Appendix II-1 for comparison.When K 0 < 1.0, the highest stress appears at the left and right ends of the chamber.This implies that the chamber has the most risk of stability at its horizontal end.With the increase of lateral pressure coefficient, the risk point gradually moves towards the vertical end.The magnitude and variation of peak stress are compared in Table 1.These results show that the effect of internal air pressure causes a significant reduction of stress around the chamber.The change of peak stress refers to the ratio of peak stress with and without internal air pressure if other parameters are kept the same, thus being a dimensionless number.

| Validation of analytical solutions for a two-layer chamber
The two-layer chamber still has the radius of 2 m for a circular cavern.The lining layer is assumed to use the concrete C60.The layer thickness is 0.5 m.The layer material has the elastic modulus of 40 GPa and a Poisson's ratio of 0.2 according to the reference values given by Wang et al. 44 A range of 20 m is chosen for simulation.Figure 9 shows the stress cloud for this two-layer chamber with lining.Symbol L represents the lining layer.Table 2 summarizes the distribution of peak stress with and without internal air pressure.Appendix II-2 gives those figures at other two lateral pressure coefficients.
The role of lining layer is observed by comparing the stresses in Figure 8 for single-layer chamber and in Figure 9 for two-layer chamber.They show that the lining layer bears the main load from the internal air pressure and thus maintains the chamber stability.a In this column, ✖ means without internal load and ✔means with internal load.

| Validation of analytical solutions for a three-layer chamber
The same geometry and parameters as the two-layer chamber are used here, but the lining thickness is 0.2 m and the grout layer thickness is 0.5 m.The grout layer has an elastic modulus of 25 GPa and a Poisson's ratio of 0.25. 45Figure 10 shows its stress cloud, where symbol G denotes the grout layer.Appendix II-3 presents more figures at the other two lateral pressure coefficients for detailed comparison.
The maximum stress distribution is in X-shape when the number of layers is even, while the maximum stress distribution is basically the same as that of a single-layer structure when the number of layers is odd.This is because the ends of the odd-numbered layers are | 4211 subjected to higher stress.The maximum stress is at the left and right for the lateral pressure coefficient less than 1.0, and at the top and bottom for the lateral pressure coefficient greater than 1.0.The stress in the evennumbered layers moves in the corresponding direction and induces this X-shaped stress distribution.The magnitude and variation of peak stress are summarized in Table 3.It is observed that the grout layer does reduce the stress in the rock mass, but this reduction of stress is not significant compared to the role of lining layer.Therefore, the lateral pressure coefficient and internal gas pressure have a significant impact on chamber stability.

| THE ROLES OF LINING AND GROUT LAYERS
Previous sections investigated the possible influences of the lining and grout layers on the chamber through the analytical solutions of stress and displacement for singlelayer, two-layer, and three-layer chambers.The correctness of stress distribution and peak stress were compared with the stress clouds obtained in this study and other scholars, but the magnitude of stress at each point is not verified yet.This section further explores the feasibility of these analytical solutions by comparing the point stress in the horizontal direction with that in the vertical direction.Their radial and circumferential stresses in each layer are compared at different lateral pressure coefficients.

| Single-layer chamber
For the dimensionless analysis of all data, the stress concentration factor is defined as σ r /γH for radial stress and as σ θ /γH for circumferential stress.The stress concentration factor is positive for compressive stress and negative for tensile stress.Figure 11 illustrates the radial and circumferential stresses in the horizontal and vertical directions for a single-layer chamber.0°denotes the horizontal direction and 90°refers to the vertical direction.RS is the ratio of radial stress (σ r /γH) and CS is the ratio of circumferential stress (σ θ /γH).The angular scale 0 means no air pressure and 1 means the existence The magnitude and variation of peak stress in the three-layer chamber.| 4213 of air pressure.This figure shows that the analytical solution generally agrees well with the FEM simulation, but small differences are observed.The details are summarized as follows:

Lateral
(1) The radial stress concentration factor gradually approaches to the lateral pressure coefficient in the horizontal direction and to 1.0 in the vertical direction.The circumferential stress concentration factor has an opposite change pattern.The radial stress concentration factor grows faster at higher lateral pressure coefficient.(2) The air pressure in the chamber increases the radial stress but reduces the circumferential stress even to tensile stress.The maximum radial stress does not occur at the edge of chamber, but moves inwards (see Figure 11 RS 1 -0°).

| Two-layer chamber with lining layer
The radial and circumferential stresses of the two-layer chamber are presented in Figure 12.The following are observed from these curves: (1) The radial stress concentration factor still eventually approaches to the lateral pressure coefficient in the horizontal direction and to 1.0 in the vertical direction.The lining layer leads it to a further increase of convergence rate.
(2) The circumferential stress is in a "bow and arrow" shape.When the lateral pressure coefficient is less than 1.0, the circumferential stress is horizontally located at the upper part of the "bow and arrow" and vertically located at the lower part.The possible mechanism on the "bow and arrow" shape is analyzed below.When the lateral pressure coefficient is 0.5, the vertical stress is greater than the horizontal stress, and the circular cavity deforms into an elliptical cavity.At this time, the ellipse has the major axis in the horizontal direction and the minor axis in the vertical direction.The horizontal circumferential direction is mostly in a tensile state, while the vertical direction is converted into circumferential compressive stress.The horizontal outermost circumferential stress of the rock layer is basically equal to the vertical load acting on the outer side, converted into a circumferential stress set of only 0.5.When the lateral pressure coefficient is 1.5, the axial stress in the horizontal and vertical directions is exactly opposite to the result at 0.5.
(3) The air pressure in the chamber reduces the circumferential compressive stress in the lining layer but increases the circumferential tensile stress in the rock layer.This increase or decrease in circumferential stress will significantly increase the local stress difference, thereby threatening the stability of the chamber.

| The role of lining
When the lateral pressure coefficient is 1.0, the chamber is under hydrostatic pressure and its stress is Equation ( 14) has the only variable of b 1 .The effects of lining and grout layers on the rock layer can be easily studied through this parameter.For a two-layer chamber, b 1 is still given in Equation ( 9) and can be reduced to the following: where 2 > 0, thus, 0 < b 1 < 2. For a single-layer chamber, Equation ( 5) can be reduced to a 1 (16)   When the chamber has no internal air pressure, b 1 is still smaller than 2 from Equation ( 14), but is 2 from Equation (16).The lining layer makes the value of Equation ( 15) smaller than that of Equation ( 16), and thus enhances the radial stress but reduces the circumferential stress.Let γH be 10 MPa and n = (r 2 − r 1 )/r 1 .The stress difference (σ ρ − σ θ ) is closely related to the Von-Mises stress, thus directly affects the chamber stability through shear failure.In general, the chamber stability is worse for bigger stress difference.The | 4215 relationship among stress difference (σ ρ − σ θ ) and shear modulus ratio (G l /G r ), lining thickness, k l , is shown in Figure 13.
Both Equation (15) and Figure 13 show that (1) b 1 in a multilayer chamber is always smaller than that in single-layer chamber and the stress of twolayer chamber is always smaller than that of singlelayer chamber.On this sense, a lining can greatly reduce stress difference (σ ρ − σ θ ) and improve chamber stability.(2) The stress difference decreases as the shear modulus ratio (G l /G r ) and lining thickness ratio (n = (r 2 − r 1 )/ r 1 ) increase.When G l /G r is small, the stress difference decreases faster.The relationship between stress difference and G l /G r follows a hyperbolic relationship.(3) Figure 12A shows that after a certain ratio of G l /G r and n is reached, continuously increasing the shear modulus of material and lining thickness is not effective in reducing the displacement and stress difference (σ ρ − σ θ ) in chamber.Figure 13B shows that G l /G r has more pronounced effect on the stress difference (σ ρ − σ θ ) than k l .For a given shear modulus ratio, smaller k l implies smaller stress difference (σ ρ − σ θ ), although this variation is not significant.Therefore, any attempt is trivial to maintain chamber stability by using a lining layer with a high shear modulus or by increasing the lining thickness without restriction or by choosing a material with low k l .
As summary, a 1 , b 1 , and b 2 are the parameters in the analytical solutions of stress and displacement in the rock layer.b 1 is related to external loads and material properties, whereas a 1 and b 2 are only related to material properties.Thus, b 1 can be regarded as a universal parameter to evaluate chamber stability.

| Three-layer chamber with lining and grout layers
The radial and circumferential stresses in the three-layer chamber are illustrated in Figure 14 in the horizontal and vertical directions, respectively.For a three-layer chamber, the accuracy of the analytical solution decreases relative to that in the two-layer chamber.However, overall results are still in good agreement with numerical simulations.The conclusions are mainly drawn as follows: (1) The radial stress concentration factor still eventually approaches to the lateral pressure coefficient in the horizontal direction and to 1.0 in the vertical direction.Compared to the two-layer chamber, the three-layer chamber approaches more quickly.(2) The shape and pattern of the circumferential stress concentration factor are consistent to those in the two-layer chamber.The curve in (a 1 -0°) tends to decrease and then increase when the lateral pressure coefficient is small, whereas the curve increases and then decreases in the two-layer chamber.The curve in (a 1 -90°) increases and then decreases for smaller lateral pressure coefficient, but the result in two layers decreases and then increases.These results indicate that the grout layer somewhat reduces the radial stress in the horizontal direction and thus reduces the horizontal displacement.| 4217 (3) The internal air pressure weakens the circumferential compressive stress in the lining and grout layers but increases the circumferential tensile stress in the lining layer only in the horizontal direction and inhibits the growth of tensile stresses in the vertical direction.The circumferential tensile stress is mainly induced by internal air pressure.

| The roles of lining and grout layers
The stress distribution in the rock under hydrostatic pressure was presented in Section 4.3 and b 1 in the threelayer chamber is given by A function f is defined as the difference of b 1 between the three-layer chamber and the two-layer chamber and is where and f < 0. This means that b 1 is smaller in a three-layer chamber than in a two-layer chamber.
Let γH be 10 MPa, m = (r 3 − r 1 )/r 1 , and l = (r 3 − r 2 )/r 2 .The relationship between stress difference (σ ρ − σ θ ) and G l /G r , G g /G r , k l , k g is shown in Figures 14-17.Figure 14 presents the relationship among stress difference (σ ρ − σ θ ), shear modulus ratio and chamber thickness.Figure 16 shows the effect of G l /G r versus k l , k g on the stress difference (σ ρ − σ θ ) and Figure 17 shows the effects of G g /G r plus k l , k g on the stress difference (σ ρ − σ θ ).
Comparing Figures 12 and 14 can draw similar conclusions but some differences are highlighted here: (1) Figure 14B shows that the stress difference decreases faster when G g /G r is higher.As the thickness of grout  layer increases, the effect of grout layer thickness on the stress difference (σ ρ − σ θ ) becomes less and less.The grout layer plays more role in maintaining chamber stability than the lining layer.(2) The effects of G l /G r , k l , and k g on the stress difference are significant as shown in Figure 16. Figure 17 further shows that G g /G r has a significant effect on the stress difference than k l and k g .When a material with low k l and k g is chosen, the stress difference (σ ρ − σ θ ) is lower and the chamber is more stable.

| Influence range of internal gas pressure on chamber stress
When the lateral pressure coefficient is 1.0, the chamber is under hydrostatic pressure and its stress is described by Equation ( 14).The influence range of internal air pressure on the chamber stress can be calculated by the decay of additional radial stress.This is because the internal air pressure has a very small impact at the infinite distance.The influence range of internal air pressure on the chamber stress can be obtained by setting a tolerance threshold.The radial stress error rate function is defined as SRE and the tolerance threshold is e%.
where b 1 0 denotes the value at no internal air pressure and b 1 1 represents the value at internal air pressure.
For a single-layer chamber, b 1 0 and b 1 1 can be written as For a two-layer chamber, b 1 0 and b 1 1 become    (1) The analytical solutions of stress and displacement were derived for single-layer, two-layer and threelayer chambers under any lateral pressure coefficient and any internal air pressure.These analytical solutions are mathematically self-consistent and can be used to quickly analyze and evaluate the distribution of stress and displacement in each layer of a multilayer chamber.(2) Internal air pressure significantly modifies the location and magnitude of risk points in the CAES chamber.Generally, internal air pressure increases radial compressive stress and decreases circumferential compressive stress even to tensile stress.Under internal air pressure, circumferential stress changes much greater than radial stress and becomes the dominant factor to chamber stability.(3) At different lateral pressure coefficients, the rock layer in the three-layer chamber is the least stressed, least displaced, and most stable compared with single-and two-layer chambers.The stress difference in chamber is inversely proportional to their thickness and correlated with the shear modulus and Poisson's ratio of these layers.The influence range of internal air pressure on the chamber stress is bigger for a more complex chamber structure.(4) The comprehensive parameter b 1 obtained in this study can quickly evaluate the redistribution of stress and displacement in rock mass, lining layer and grout layer under any lateral pressure coefficient, internal air pressure and material properties.The influence range of internal air pressure on chamber stress is analytically obatined.These two parameters can be used to evaluate and enhance the CAES chamber stability through the selection of lining and grouting layers.

5
Comparison of our results with those of Dong et al.41 (A) radial stress concentration σ r /γH, (B) circumferential stress concentration σ θ /γH along horizontal direction, (C) radial stress concentration σ r /γH, and (D) circumferential stress concentration σ θ /γH along vertical direction.SUN ET AL.Comparison of our results with those of Du et al. 42 (A) radial stress concentration σ r /γH, (B) circumferential stress concentration σ θ /γH along horizontal direction, (C) radial stress concentration σ r /γH, and (D) circumferential stress concentration σ θ /γH along vertical direction.F I G U R E 6 A sketch model for numerical simulations on a three-layer chamber.F I G U R E 7 Stress-strain curve of rock in a coal mine.

F I G U R E 8
Stress cloud of single layer model.T A B L E 1The magnitude and variation of peak stress in the single-layer chamber.

10
Stress cloud in the three-layer chamber model.SUN ET AL.

F
I G U R E 12 Comparison of stress concentration factors in two-layer chambers.SUN ET AL.

F
I G U R E 13 Relationship between stress difference and related parameters.(A) Stress difference, shear modulus ratio and thickness and (B) stress difference, shear modulus ratio and k l .F I G U R E 14 Comparison of stress concentration factors for three-layer chambers.SUN ET AL.