Abstract
The contradiction between supply and demand of energy leads to more and more attention on the large-scale energy storage technology; Compressed Air Energy Storage (CAES) technology is a new energy storage technology that is widely concerned in the world. The research of coupled heat transfer and seepage in fractured surrounding rocks is the necessary basis to evaluate the operation safety and effectiveness of CAES. Current studies point to the possibility of cracking in concrete liner seals, but the thermodynamic processes and leakage characteristics of compressed air in the presence of cracking and the heat transfer characteristics of seepage have not been addressed and reported. In order to investigate the leakage, the gas seepage and heat transfer law in fractured rock when the hard rock CAES gas reservoir seal cracks, the COMSOL fracture Darcy module, and the non-Darcy Forchheimer model are used as the constitutive seepage. The global ODE is used to calculate the thermodynamic process of compressed air in gas storage with coupled seepage and heat transfer process. The pressure and temperature of compressed air are obtained as the unsteady boundary of the seepage heat transfer model. A program for calculating the seepage and heat transfer characteristics of fractured surrounding rock in the CAES gas reservoir is established. On this basis, with the proposed Suichang CAES cavern as the background, the seepage and heat transfer characteristics of the fractured surrounding rock of the gas storage are studied. The results showed that when there are fewer cracks in the lining and surrounding rock of the air reservoir, the air pressure decreases due to a small amount of air leakage after 30 operation cycles, and the leakage rate of each cycle is 0.7% of the gas storage capacity, but it still meets the engineering requirements. If the plant is operating under these conditions, the charging rate will need to be increased by 1.2 kg/s per cycle charging stage. In the discharging and power generation phase, the high-pressure air that previously percolated into the rock mass cracks could flow back into the air storage through the lining cracks. Therefore, it is incorrect and unreliable to consider the gas which flows out from the inner surface of the lining as unusable. When the lining crack width is less than 0.3 mm, the seepage flow is Darcy flow and the non-Darcy effect can be ignored; when the lining crack width is greater than 0.5 mm, the non-Darcy effect of seepage cannot be ignored. The gas velocity in the surrounding rock fracture medium is on the order of 0.01 m/s with an influence range of over 100 m, and the gas velocity in the pore medium is on the order of 10-6 m/s with an influence range of 50 m. The findings of this study contribute to a better understanding of the interaction between the thermodynamic properties of compressed air and the seepage heat transfer process in compressed air storage underground reservoirs, as well as the gas leakage process in the event of liner seal cracking.
1. Introduction
Compressed Air Energy Storage (CAES) is a new technology suitable for large-scale power energy storage; it is one of the key technologies to solve the instability and intermittency of large-scale clean energy such as wind and solar energy [1–4]. At present, the large-scale commercial Compressed Air Energy Storage in operation uses salt caverns as a gas reservoir; however, the site selection of salt rock underground gas storage has great geographical limitations; manual excavation of hard rock gas storage can solve the problem of location limitation of underground gas storage ideally [5]. In the operation of the CAES power station, the gas reservoir is in the state of high temperature and high pressure with frequent cycle changes. Chen et al. [6] pointed that granite, as one of the alternative surrounding rocks for underground gas storage, has a cracking threshold temperature of 60-70°C under temperature damage; at the same time, minor cracks also exist in the concrete lining, which is used for the sealing base of the underground gas reservoir, due to its own defects and the limitation of construction technology. The gas will leak along the cracks in the concrete lining, and cracks in the rock mass after the sealing layer cracks under the high-pressure action of the internal compressed air. When the variable temperature compressed air enters the concrete lining and surrounding rock, it will lead to the uneven change of the internal temperature and then cause the temperature stress. After a long-term operation, it may endanger the normal performance of the sealing function of the gas storage and even affect the safe operation of the underground gas storage.
The temperature and pressure variations of compressed air in the underground storage of compressed air energy have been a hot topic of interest and focus for scholars today; Raju and Khaitan [7], Kushnir et al. [8], Zhou et al. [9], and Xia et al. [10] have established the analytical equations for temperature and pressure by means of mathematical analytical solutions. However, they focus more on the effect of different thermal treatments (adiabatic, thermal conductive, or thermostatic) of the cavity wall on the energy storage in the cavity. Kim et al. [11] pointed out that during the frequent cycle of charging and discharging, compressed air can leak from the initial defects of the structure or weak construction joints. Pornkasem et al. [12] studied the fracture initiation and extension direction of different high internal pressure underground gas storage reservoir rock bodies. The seepage and heat transfer characters in the surrounding rock under liner seal cracking failure and the effect of leakage on the thermodynamic process of compressed air in the reservoir have not been studied and reported yet.
The temperature variation of the surrounding rock and concrete lining of the air reservoir is controlled by the heat conduction of the solid medium and the convective heat transfer between the gas and solid in the fractures and pores. Due to the low viscosity of gas, non-Darcy flow may occur in the flow through fractures. However, there are few reports on the nonseepage heat transfer characteristics of fractured rock mass. In the existing literature, researches on seepage and heat transfer of fractured rock mass mostly focus on Darcy seepage and heat transfer characteristics of water in fractures [13–17]. The non-Darcy seepage characteristics of gas have been widely reported in the research of coalbed methane [18–20]. However, these studies seldom consider the influence of seepage heat transfer.
In terms of numerical simulation methods of gas seepage in fractured rock mass, most commercial software at present have the Darcy seepage model but do not have the non-Darcy seepage function. The results show that using the non-Darcy correction factor to modify the Darcy model is a simple and reliable method to calculate the non-Darcy seepage [21].
The physical process of the research content of this paper is shown in Figure 1. The research process was as follows: Firstly, the non-Darcy seepage theory was introduced, and the expression of non-Darcy permeability corrected by the non-Darcy correction factor without the flow rate term was given. Then, the thermodynamic process of compressed air in compressed air storage reservoir considering leakage and heat transfer process was introduced. Then, the non-Darcy correction factor was used to modify the Darcy model of fracture in the groundwater module of COMSOL numerical analysis software in this study, so as to realize the non-Darcy seepage analysis function of fracture in COMSOL. The seepage heat transfer model with a sand-filled fissure was used to verify the model of this paper and prove the correctness of the model of this paper. Finally, the background of a compressed air storage power station in Suichang, Zhejiang Province, was studied. The global ODE coupled with seepage heat transfer was used to calculate the variation of gas temperature and pressure in the gas storage under the state of seepage. The leakage process and seepage heat transfer characteristics of the liner cracking were studied, and the influence of the relative position of the cracking location and the surrounding rock fracture on the seepage heat transfer process was analyzed.
The advantages of this study are the elucidation of leakage when the liner seal cracks and fails and the effect of cracked leakage on the thermodynamic processes of compressed air in the storage reservoir, as well as the study of the seepage heat transfer characteristics of compressed air after cracking. At the same time, a non-Darcy correction factor without the flow velocity term is proposed in the theory, which enables any numerical simulation software containing the Darcy model to calculate the non-Darcy seepage process by this correction factor, a mathematical model of compressed air thermodynamics considering seepage and heat transfer is established, and the results of the analytical solution are used as the unsteady state boundary to calculate the seepage heat transfer process in the gas storage reservoir. The disadvantage of this study is that the lining cracking conditions (cracking location, number of cracks, cracking width) and the distribution of cracks in the surrounding rock are autonomous random assumptions, which may not be the same as the actual situation. And the dynamic change of fracture width with the process of gas filling and deflation is not considered.
2. Theory of Non-Darcy Flow and Heat Transfer in Fractured Rock Mass
2.1. Non-Darcy Seepage Theory
In the case of gas seepage in fractures or porous media, the flow rate is proportional to the pressure gradient and inversely proportional to the ratio of permeability and viscosity, which is Darcy’s law. When the gas flow rate is high, the inertia force of the gas flow is large, which will produce the phenomenon of deviation from Darcy’s law. The Forchheimer equation describes the non-Darcy seepage by adopting the square term of velocity to represent the inertia force of the fluid, and the formula is as follows: where is aerodynamic dynamic viscosity (Pa s), is the permeability of porous or fractured media (m2), is the gas seepage velocity (m/s), is the non-Darcy coefficient (dimensionless), , and is the density of the gas (kg/m3). Its physical significance is that the seepage resistance is composed of two parts: the first part is viscous resistance, which is proportional to velocity, and the second part is nonlinear seepage inertia additional resistance.
It is pointed out that the Forchheimer equation could be modified by the Darcy law to reflect nonlinear seepage characteristics [21]. The modification can be written as where is the non-Darcy correction factor; its calculation expression is shown as follows: When the velocity is very small, is equal to 1; the above Equation (2) is still Darcy’s law.
Non-Darcy permeability is defined as follows:
By substituting the non-Darcy seepage rate into Equation (2), the Darcy seepage expression form reflecting the characteristics of the non-Darcy seepage can be obtained. Since in Equation (4) contains the velocity term, this calculation formula uses permeability with the velocity term to calculate velocity, which is not applicable in some numerical software, so Equation (4) needs to be further modified. After the following derivation and transformation, another expression formula of permeability of the non-Darcy model is proposed. Details of ’s derivation process can be found in the Appendix. The expression for without the flow rate term is shown as follows:
2.2. Analytical Equation of Seepage and Heat Transfer Boundary of CAES Gas Reservoir
The boundary of seepage and heat transfer analysis of surrounding rock of an underground gas reservoir is closely related to the temperature and pressure of compressed air. The equation of temperature change in the gas storage during the operation of a CAES power station and its derivation process has been introduced in detail in literature [22]. The differential equation of temperature change is directly quoted here as follows: where is the compressed air temperature in the gas reservoir (K), is the air charging mass rate (kg), is the specific heat capacity (J/(kg·K)), is the temperature at which the gas is charged (K), is the volume of gas cavern (m3), is the compressed gas pressure in the air reservoir (Pa), is the convective heat transfer coefficient between air and tunnel wall (W/(m2·K)), is the convective heat transfer area of the tunnel wall (m2), is the temperature of the cave wall (K), and is the quality of compressed air in the gas reservoir (kg).
The differential equation of compressed air mass variation in the gas reservoir is shown as follows: where is the discharging air mass flow rate (kg/s) and is the mass flow rate of gas leakage from the air reservoir (kg/s).
The differential equation of compressed air pressure variation in the gas storage can be obtained from the idea gas state equation as shown: where is the gas constant (J/(kg·K)).
The gas reservoir leakage rate is calculated by the following formula: where is the thickness of the 2D model, which is the length of the cave (m), and is the arc length of inner surface of circular lining sealing layer (m).
3. Secondary Development and Verification of Calculation Program
COMSOL has advantages in multifield coupling, such as good convergence, fast calculation speed, and a large number of physical and mechanical models. In the latest version of COMSOL, the porous media seepage analysis module has included the non-Darcy seepage model, but the fracture seepage analysis module does not have the calculation function of the non-Darcy seepage. Therefore, based on the non-Darcy seepage analysis method mentioned above in this paper, the non-Darcy seepage analysis function of fracture of the COMSOL software platform is realized.
3.1. Secondary Development Process of Calculation Program
The global ODE and the modification of COMSOL Darcy’s law module were used for the secondary development of the calculation program. The secondary development idea is as follows, and the calculation process of the development program is shown in Figure 2. (1)The constitutive model of fracture seepage under the Darcy’s law module was modified by , and the constitutive model was changed into the non-Darcy Forchheimer equation(2)The global ODE in the mathematics module was used to calculate the change process of the compressed gas mass (Equation (10)), temperature (Equation (9)), and pressure (Equation (11)) in the gas reservoir; the compressed air mass was coupled with the gas leakage on the inner surface of the lining of the seepage module, and the change of temperature was coupled with the heat transfer on the inner surface of the lining of the heat transfer module(3)The gas pressure and temperature in the gas storage calculated by the global ODE were, respectively, taken as the wall pressure boundary of the modified seepage field and the wall temperature boundary of the heat transfer field, and the gas seepage flow and heat transfer at the wall were calculated based on these boundary conditions
3.2. Program Verification
To verify the correctness of the program, the physical test of seepage and heat transfer in sand-filled fractures in literature [17] was adopted to verify the model. The verified model is shown in Figure 3 with a size of ; four fractures V1, V2, H1, and H2 are set in the model. The upper ends of V1 and V2 were water injection ports. The temperature of the injected water was 16.4°C; the lower end was the water outlet; the left side of the model was the temperature heat source 95°C; the upper boundary, the right boundary, and the lower boundary were the adiabatic boundaries; and the solid pore seepage was not considered here. The initial pressure condition was a uniform pressure of 0 MPa, the calculation parameters of the model are shown in Table 1.
The outlet temperatures of V1 and V2 fractures increased by 40.4 and 15.2°C, respectively, in the literature physical experiment; the outlet temperatures of V1 and V2 were, respectively, increased by 40.04°C and 16.4°C by using the redeveloped verification program in this paper; the relative error was 0.89% and 7.89% with the measured values obtained from the physical experiment, in which the numerical simulation results of the outlet temperature of V1 and V2 obtained in literature [17] and the physical experiment results showed an error of 6.3% and 10.1%, respectively. The error of the modified model in this paper was lower than that of the original literature, indicating that the accuracy of this modified model was higher than that of the original model. Figure 4 shows the comparison of the temperature of vertical fracture V1 at different heights calculated by the modified model in this paper with the measured results of the physical test and the calculated results of the numerical model in literature [17]. The temperature distribution of fracture V1 at different heights calculated by the modified program in this paper was highly consistent with the results of a physical experiment and mathematical model calculation in literature [17]. Moreover, the modified model in this paper had a better fit with the physical experiment in the literature. According to the comparison results of the temperature distribution at different heights of fracture V1 and outlet values of V1 and V2 between the revised model in this paper and the physical tests and numerical models in the original literature, it can be seen that the revised model in this paper was correct and reasonable in calculating the process of seepage and heat transfer with high accuracy.
4. Seepage and Heat Transfer Characteristics of Fractured Surrounding Rock in CAES Air Reservoir
4.1. Numerical Model and Calculation Conditions
In order to study the seepage and heat transfer characteristics of fractured surrounding rock of CAES gas storage in the operation of a power station under the condition of high-pressure air leakage, a numerical analysis model for gas permeability of the air reservoir was established based on the underground air reservoir project of a CAES power station planned to be built in Suichang, Zhejiang Province. The underground air reservoir was a tunnel-type with circular cross-section and located in the granite strata, the cross-section diameter of the air reservoir was 15 m, and the size of the numerical model was . The research showed that concrete can be used as sealing layer of gas in a high-pressure underground gas storage when its permeability was low enough [11]. Therefore, concrete was considered as the sealing layer in this calculation. Assume that there were 3 fractures in the lining sealing layer, the location of fractures was random, and the fracture lengths were 50 m. Considering that there will still be a certain degree of contact surface between the concrete lining and the surrounding rock after the completion of construction [23], the contact surface between the concrete lining and the crack was set as a fracture attribute. In order to investigate the seepage and heat transfer characteristics of fractures in surrounding rock, and according to the geological prospecting report, the structure of surrounding rock mass in the roadway of +240 m, +210 m, +160 m, +120 m, and +70 m in the five middle sections was complete, and the joints and fractures were, respectively, undeveloped, generally developed, and closed without filling; 15 long fractures were assumed in the cross-section and randomly generated by the Monte Carlo method. In order to study the influence of the relative position of surrounding fractures and lining fractures on the process of seepage and heat transfer, the random process ensures that one rock fracture was in alignment with the lining crack and another is misplaced with the lining crack. The model grid is shown in Figure 5.
The calculation process included the initial charging process, so the initial pressure of the seepage field was 0.1 MPa and the initial temperature was 15°C. In order to avoid the influence of the model size on the calculation results, the outer boundary of the model was set as an infinite element domain, and the outer boundary of the infinite element domain was set as the boundary of the room temperature and atmospheric pressure, namely, 0.1 MPa and 15°C. The inner surface of the concrete lining was a custom global ODE unsteady temperature pressure boundary coupled with seepage flow and heat transfer. According to the experience of the hard rock air reservoir in the Hunan Pingjiang test reservoir, the upper and lower limits of operating pressure were determined to be 7-10 MPa, and then, the charging and discharging rates were determined according to the time of the operation stage and the upper and lower limits of pressure. The calculation process included the initial charging stage and 30 charging and discharging cycles. The initial charging stage was 1 d, the charging rate was 53 kg/s, and the air pressure in the air reservoir reached 10 MPa. In the operation process, a complete cycle time was 1 d, which were, respectively, the charge stage (8 h), high pressure storage stage (4 h), discharge stage (4 h), and low-pressure storage stage (8 h). During stable operation, the charging rate was 47.5 kg/s, and the deflating rate was 95 kg/s. The parameters of the fractured surrounding rock model of the Suichang gas storage reservoir are shown in Table 2.
Calculation conditions: (1) air leakage of concrete lining sealing layer was not considered; (2) the cracking and air leakage of the concrete lining sealing layer were considered.
The model calculation was based on the following assumptions: (1) local thermal equilibrium assumption, (2) the fractures in the sealing layer of concrete lining was a linear crack, and (3) the influence of pressure variation on the gap width between the lining and the surrounding rock during operation was ignored.
4.2. Influence of Lining Cracking on Temperature and Pressure Variation in Gas Storage
Figure 6 shows the changes of compressed air temperature and pressure in 30 cycles of the air reservoir under the conditions of no gas leakage and cracking leakage of lining. When the concrete lining was used as the sealing layer, the pressure operating range was 7-10 MPa in 30 cycles, and the concrete lining cracked according to the assumed conditions; the air pressure in the gas storage was reduced by 1.5 MPa after 30 operation cycles. The leakage of concrete lining had little influence on the compressed air temperature. When the current charging and discharging rate was adopted, the maximum temperature in the initial charging process would reach 401.8 K. Then, the temperature gradually decreased with the process of charging and discharging alternately, and the temperature varied from 266.5 to 300.4 K in the stable operation stage. When the lining cracking gas leakage was considered, the minimum air temperature decreases by 3 K and the maximum air temperature increases by 3 K.
Figure 7 shows the leakage rate of compressed air at the tunnel wall when the concrete lining cracks. When the concrete lining was used as the sealing layer and cracks, in the process of circulating operation, the gas leakage rate varied within the range of 0.2-0.6 kg/s; the average gas leakage rate in the whole operation process was 0.4 kg/s, accounting for 0.8% of the charging rate; and the leakage amount of each cycle accounts for 0.7% of the total gas mass in the air reservoir. The leakage rate per cycle was less than 1%; this cracked concrete lining sealing layer still meets the engineering needs. The gas seepage flux distribution in the pore medium at the cavern boundary under four different operating states in the 15th cycle is shown in Figure 8. It was worth noting that at the end of the charging stage and the end of the high-pressure air storage stage, the gas seepage flux forms a trough around the fractures of the lining. The reason was that after the lining cracks, the air seepage pressure at the fractures was high, which led to the increase of the air pressure in the lining near the fractures and the decrease of the pressure gradient with the internal compressed air. The seepage velocity decreased, and then, the seepage flux decreased. Moreover, in the end of the discharge stage, the gas seepage flux appeared positive in a long range; it showed that the air in this stage was percolating into the cavern due to the low gas pressure inside the cavern. The results showed that when calculating the leakage of gas storage, it was necessary to consider whether there is supplementary effect of air backflow in the low-pressure stage during the circulating operation.
4.3. Spatial Distribution of Seepage and Heat Transfer Characteristics
Pressure and temperature distribution of gas seepage in fractured rock mass in the end of the charging stage of the initial charging and the 5th, 15th, and 30th cycles is shown in Figures 9 and 10. The distribution of fractures had a great influence on the distribution of seepage pressure. The fractures in the surrounding rock were the main channel of air seepage. The high-pressure air flowed outward from the fractures in the lining and surrounding rock and at the same time flowed into the surrounding rock medium from the gap between the lining and surrounding rock. The main high pressure seepage channels were opposite fractures and misplaced fractures connected with lining fractures, which form high pressure extension zones.
The temperature range of the surrounding rock was small, and the temperature at the lining fracture changed rapidly with the compressed air temperature in the cavern. A high temperature extension zone would be formed in the surrounding rock facing fracture and dislocation fracture, but the extension distance was very limited, which was 6 m and 5 m. With the circulation process, the temperature of lining and surrounding rock would form an isothermal circle of alternating high and low temperatures. After stable circulation, the temperature of lining would be lower than that of surrounding rock in the process of change.
Air seepage velocity distribution in the porous medium of fractured surrounding rock is shown in Figure 11; it showed that the higher air seepage velocity in the surrounding rock was the cracking point of the lining and the junction of the surrounding rock fissure and the surrounding rock. The maximum flow velocity decreased with the increase of operating time. The higher air seepage velocity in the initial stage was the junction of the cracking point and the surrounding rock; at the later stage of operation, the higher velocity was the end of the surrounding rock fractures, the maximum flow rate decreased from to , and the range of gas seepage also gradually expanded outward from the initial crack point of lining and surrounding rock fractures. After 30 cycles, the range of gas seepage reached 50 m along the porous media outside the lining, and the influence range along the fracture media was related to the fracture connectivity.
4.4. Seepage and Heat Transfer Evolution Process Analysis of Fractured Surrounding Rock of Air Reservoir
In order to investigate the influence of location of lining fractures and surrounding rock fractures on the seepage and heat transfer process of the air reservoir, monitoring points were selected from the lining, lining cracks, and the gap between the lining and surrounding rock. The distribution of monitoring points is shown in Figure 12. Lining monitoring points were marked L, monitoring points in the gap between lining and surrounding rock were marked g, surrounding rock monitoring points were marked R, and subscripts f and p denoted fractures and porous media, respectively.
Figures 13 and 14, respectively, show the variation of gas seepage velocity in lining fractures and surrounding rock fractures. The gas seepage velocity in the fracture was on the order of 0.01 m/s. The lining fracture directly connected with the surrounding rock fracture had the highest flow velocity, and the seepage velocity at Lf4 was 0.008 m/s higher than that at Lf1 in the initial cycle operation and 0.002 m/s higher in the stable operation. The gas seepage velocity at Lf6 and Lf5 was 0.01 m/s higher than that at Lf3 and Lf2. Therefore, the gas seepage velocity was higher at the intersection of the lining fracture and the gap of the lining and surrounding rock. The gas seepage velocity at Lf6 was higher than that at Lf5, and the gas seepage velocity at Lf3 was higher than that at Lf2. Therefore, the gas flow velocity in the lining fractures near the surrounding rock fractures was higher than that in the lining cracks far away from the surrounding rock fractures. The gas flow velocity at Rf1 and Rf3 was much higher than that at Rf2 and Rf4. It can be seen that the gas flow velocity in the surrounding rock fractures which connected with the lining fractures was higher than that in the dislocation fractures. The maximum gas seepage velocity in Rf1 was 0.003 m/s higher than that at Rf3; that is, the maximum gas seepage velocity decreased by 0.003 m/s in the distance of 6.24 m in the opposite fracture. Fractures in the surrounding rock that were not connected to the gap had very low gas flow rates.
Temperature changes at monitoring points at lining fractures and the gap between the lining and surrounding rock are shown in Figure 15. The temperature change at each measuring point of the lining was the same as the air temperature change law in the air reservoir. Since the thickness of the lining was only 0.5 m, the temperature difference at each point was not large, with 10 K maximum. The monitoring point with the highest temperature was Lf4, where the gas flow rate was the highest and brings more heat energy, so the temperature was higher, while gf1, gf2, and gf3 had the lowest temperature because of their distance and low flow rate.
The variation of gas seepage velocity at the monitoring point in the surrounding rock porous media of air reservoir is shown in Figure 16. The gas seepage velocity in the surrounding rock porous media was on the order of 10-7 m/s. The order of gas velocity at each point 8 m away from the outer wall of lining was . At the same distance, the gas seepage velocity in the porous media around the surrounding rock directly opposite the fracture was higher than that of the pure porous media, and the gas seepage in the proximal porous media was promoted by the surrounding rock directly opposite the fracture. The seepage velocity of Rp1 was higher than that of Rp2. It can be seen that, although the seepage pressure was prevented from being transferred to the far end by the distant fracture, its dispersion effect on the seepage pressure will lead to the gas seepage velocity in the porous media near the wall being higher than that of the pure porous media. As the gas seepage velocity of Rp2 was higher than that of Rp6, it can be seen that although the distal end of the dislocation fracture was the pressure extension zone to promote the seepage, it hindered the gas seepage at the proximal end. At a distance of 14.24 m from the outer wall of the lining, the order of seepage velocity at each monitoring point was . In a remote place, the existence of a directly opposite fracture and dislocation fracture would lead to the gas seepage velocity in the surrounding rock porous media lower than that in the pure porous media. The reason was that the same seepage pressure extends farther in the fracture, leading to higher pressure in the around rock, and the pressure gradient between the rock and the outer wall of the lining was lower, so the gas seepage velocity in the porous media far affected by the fracture was lower than that of the pure porous media.
Pressure changes at monitoring points of surrounding rock porous media are shown in Figure 17. The pressure at the monitoring points Rp3, Rp4, Rp8, and Rp9 which are around the directly opposite fracture was the highest in the surrounding rock, and the pressure at the monitoring point Rp1 which is located between the far away fracture and the lining was second, followed by the pressure at the monitoring point of the porous media around the dislocation fracture. The difference between the pressure of the porous media near the proximal dislocation fracture and that of the pure porous media mainly depended on whether it was on the side affected by the lining fracture, while the pressure at the far end was obviously higher than that at the pure porous media.
The pressure changes at each measuring point of the gap between the lining and surrounding rock and the cracks in the surrounding rock are shown in Figure 18. The pressure at Rf1 was roughly the same as that at gf1. The pressure at the intersection of the lining fracture and gap was equal to the seepage pressure 8 m away from the surrounding rock directly opposite the fracture; however, the amplitude of the pressure change at the intersection was smaller than that at the directly opposite fracture. According to the pressure variation of gf2 and Rf3, the pressure at the gap points not affected by lining cracks was approximately equal to the pressure at the distance of 16.24 m surrounding rock directly opposite the fracture, and the variation range was also approximately the same. The pressure at each point of the dislocation fracture and the far away fracture was low.
5. Discussion
When calculating the high-pressure seepage of fractured rock mass, the Forchheimer equation should be used as the seepage constitutive equation when it is uncertain whether there would be non-Darcy seepage in advance. When the seepage velocity was low, the inertial resistance term of the Forchheimer equation took a small part in the pressure loss, and the difference between the flow velocity calculation and the Darcy seepage constitutive calculation was very small. When the flow velocity was high, the nonlinear term of inertia resistance led to a large proportion of pressure loss, and the nonlinear effect cannot be ignored, so the calculation result of the Forchheimer equation was more accurate. Under the lining crack condition in this study, the flow velocity at the fracture point of the lining was on the order of 0.01 m/s, the nonlinear resistance accounts for 2%, and the non-Darcy effect could be ignored. According to the report by Kim, when the crack width was between 0.2 mm and 0.5 mm, an increase of 0.1 mm in the crack width would increase the permeability by an order of magnitude. Here, the crack width of lining was artificially assumed, and the crack location and width under the action of long-term cycle temperature and pressure were the contents of the next research. If the width of the lining fracture expanded to 0.5 mm, the gas flow velocity would increase to the order of 1 m/s. At this point, the viscous resistance term of gas flow was of the same order of magnitude as the inertial additional resistance term, and the non-Darcy effect cannot be ignored.
When previous scholars studied the leakage of the air reservoir, most of them took the outflow of gas from the inner surface of the lining as the leakage amount and believed that the outflow of gas from the lining indicated that the gas was unavailable. However, through the analysis of the seepage of the fractured surrounding rock of the air reservoir, this study found that the compressed air pressure in the air reservoir changed rapidly during the circulating operation. When the air reservoir discharged to the low-pressure section, the gas leaked into the lining and fractures would infiltrate back into the cavern, resulting in a small range of gas replenishment effect. Therefore, it was not reliable to regard the gas flowing out of the inner surface of the lining seal as leakage, and it was necessary to consider the gas backflow effect caused by the change of the inner and outer pressures of the lining.
6. Summary and Conclusions
The non-Darcy constitutive model could be calculated by using the non-Darcy correction factor without the flow velocity term to modify the fracture Darcy model of COMSOL, and the modified model is verified. The results showed that the modified method can be used to simulate the non-Darcy seepage well. The unsteady boundary of compressed air temperature and pressure coupled with a seepage heat transfer term in the air reservoir was defined by using a custom global ODE method and applied to the inner boundary of the tunnel. The seepage and heat transfer characteristics of non-Darcy gas in cracks surrounding rock were calculated when the lining cracks. (1)After cracking of the concrete lining of the CAES gas storage in the example of this study, under the condition of 3 fractures which is 0.3 mm wide and 50 m long on the lining, 0.7% of the compressed air in the cavern was lost in each operation cycle, which still met the engineering needs. After 30 operation cycles, the compressed air pressure was reduced by 1.5 MPa. If the air reservoir was operating under this condition, the charging rate should be increased by 1.2 kg/s during each cycle of charging stage. After the lining cracks, the gas seepage rate in the lining around the cracking point decreased, and in the discharging stage, the gas in the lining would flow back into the cavern for supplement(2)During the seepage and heat transfer in the fractured surrounding rock of the air reservoir, the gas seepage velocity in the fractured medium was on the order of 0.01 m/s, and the gas seepage velocity in the porous media of surrounding rock was on the order of 10-6 m/s. The influence range of heat transfer was 6 m, and the influence range of gas seepage was 50 m in the surrounding rock porous medium. The influence range of gas seepage in the surrounding rock fractures was related to the fracture connectivity, and the pressure distribution formed a high-pressure extension zone along with the fractures in the surrounding rock that were connected with the gap of the lining and rock. When the crack width of the lining was 0.3 mm, the gas seepage in the fracture is a Darcy flow, and the non-Darcy effect can be ignored. When the fracture width reached to 0.5 mm, the non-Darcy effect could not be ignored(3)The relative position of the surrounding rock fracture and lining fracture had great influence on the distribution of seepage pressure and gas seepage velocity field in the surrounding rock. When the surrounding rock fracture and lining fracture were directly opposite, the gas would flow through the fracture at a high speed, resulting in a high-pressure extension zone, and the gas velocity in the proximal surrounding rock affected by the fracture was higher, but the gas velocity in the distal surrounding rock was lower. When the surrounding rock fracture was dislocated with the lining fracture, the high-pressure extension zone would also be caused, but the gas velocity in the proximal and distal surrounding rock was low. When the surrounding rock fracture was far away from the lining, the pressure in the porous medium between the fracture and the lining was higher and the gas flow rate was faster
Appendix
The non-Darcy correction factor without the flow rate term is derived as follows: the following expression can be obtained from the transformation of Equation (1):
After multiplying around Equation (A.1) together, the above equation is transformed into
Multiply the right side of Equation (A.2) and then transfer and simplify it:
Substituting Equation (A.3) into Equation (4), another expression of can be obtained:
Nomenclature
Symbols| : | Pressure gradient of gas seepage (Pa/m) |
| : | Aerodynamic dynamic viscosity (Pa·s) |
| : | Permeability of porosity (m2) |
| : | Non-Darcy coefficient |
| : | Density (kg/m3) |
| : | Gas seepage velocity (m/s) |
| : | Non-Darcy correction factor |
| : | Temperature (K) |
| : | Specific heat capacity (J/(kg·K)) |
| : | Volume of gas cavern (m3) |
| : | Convective heat transfer coefficient (W/(m2·K)) |
| : | Convective heat transfer area (m2) |
| : | Thermal conductivity of surrounding rock (W/(m·K)) |
| : | Quality (kg) |
| : | The gas constant (J/(kg·K)) |
| : | length of the cave (m) |
| : | Cavity radius (m) |
| : | Lining thickness (m) |
| : | Width of lining fractures |
| : | Porosity |
| f: | Fracture |
| r: | Rock |
| l: | Lining |
| g: | Gap between lining and rock |
| w: | Cave wall |
| in: | Charging |
| out: | Discharging |
| e: | Leakage |
| N: | Non-Darcy |
| 0: | Initial. |
Data Availability
In our manuscript, the data of the validation model and the numerical model used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare no conflict of interest.
Acknowledgments
This project was financially funded by the National Natural Science Foundation of China (Grant No. 51778070) and Hunan Provincial Innovation Foundation for Postgraduate (Grant No. CX20190662).