**1. Introduction**

Large scale energy storage systems are required to facilitate the penetration of variable renewable energies in the electricity grids [1–4]. Underground space from abandoned mines can be used as underground reservoirs for underground pumped storage hydropower (UPSH) and compressed air energy storage (CAES) systems [5–11]. Pumped storage hydropower (PSH) is the most mature large-scale energy storage technology, and the round trip efficiency is typically in the range of 70–80% [12,13]. Diabatic CAES (D-CAES) is an alternative to PSH that requires lower capital costs and the round trip energy efficiency is around 40–50% [14]. D-CAES systems use natural gas to heat the compressed air in the decompression period. However, Adiabatic CAES (A-CAES) allows the storage of the thermal energy during the air compression period, avoiding the consumption of natural gas, therefore increasing the round trip energy efficiency up to 70–75% [15–17].

**Citation:** Prado, L.Á.d.; Menéndez, J.; Bernardo-Sánchez, A.; Galdo, M.; Loredo, J.; Fernández-Oro, J.M. Thermodynamic Analysis of Compressed Air Energy Storage (CAES) Reservoirs in Abandoned Mines Using Different Sealing Layers. *Appl. Sci.* **2021**, *11*, 2573. https://doi.org/10.3390/app11062573

Academic Editor: Luisa F. Cabeza

Received: 16 February 2021 Accepted: 11 March 2021 Published: 13 March 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

During the operation of A-CAES plants, the air is stored in the underground reservoir at high pressures during the charge period. Then, in the discharge period, the pressurized air is released, heated (stored heat), and expanded to generate electricity as the pressure within the reservoir is reduced. During the storage period, between charging and discharging, the air pressure and temperature vary depending on the storage time. The amount of stored energy depends on the reservoir volume and the thermodynamic conditions. The air temperature variations depend on the thermal conductivity of the sealing layer, concrete lining and rock mass. Therefore, the air temperature and pressure fluctuations are essential to design the reservoir volume, to select the appropriate compressor and turbines and to ensure safe operating conditions. Kushnir et al. [18] developed a study to determine the temperature and pressure variations within compressed air energy storage caverns and showed that the heat transfer reduces the temperature and pressure variations during the compression and decompression periods, leading to a higher storage capacity. Thermodynamic models were performed to determine temperature and pressure variations within adiabatic caverns of compressed air energy storage plants [19]. The results showed that the storage volume is highly dependent on the air maximum to minimum pressure ratio. Safaei and Aziz [20] carried out a thermodynamic analysis of three compressed air energy storage systems. They concluded that A-CAES with physical heat storage is the most efficient option with an exergy efficiency of 69.5%.

A pilot-scale demonstration of A-CAES was built in Switzerland [21,22]. The cavern was 120 m long and 4.9 m in diameter with a usable volume of 1942 m<sup>3</sup> . A packed-bed of rocks for thermal energy storage was located inside the cavern. Charging temperatures of 550 ◦C and air pressures between 0 and 0.7 MPa were employed in the test and constant heat transfer coefficients of 5 and 10 W m−<sup>2</sup> K <sup>−</sup><sup>1</sup> were considered in the simulations. The estimated round trip efficiency of the pilot plant reached 63–74%. Schmidt et al. [23] analyzed the effect of the cyclic loading on the geomechanical performance of CAES systems for lined and unlined tunnels in closed mines. The cycling loading operation was simulated for 10,000 cycles considering an operating pressure within the reservoir from 4.5 to 7.5 MPa. They observed moderate deformations and small thickness of plastic zones, while an increase of the initial volume of less than 0.5% was reached.

Jiang et al. [24] carried out experimental and numerical investigations of a lined rock cavern to analyze the thermodynamic process within the underground cavern. The volume of the cavern was 28.8 m<sup>3</sup> with 5 m length and 2.9 m in diameter. A 2 cm thick RFP was employed as sealing layer around the compressed air. The air pressure reached 6 MPa at constant mass flow rate (0.12 kg s−<sup>1</sup> ). Zhou et al. [25] developed a numerical simulation for the coupled thermo-mechanical performance of a lined rock cavern for CAES systems. An air pressure range from 4.5 to 6.9 MPa was considered for 100 cycles with a thickness of concrete lining of 0 and 0.1 m. They concluded that significant air temperature fluctuation was observed in sealing layer and concrete lining and no fluctuation occurred in the rock mass. An analytical study was proposed to estimate the geomechanical performance caused by air temperature and pressure variations in CAES systems [26]. Khaledy et al. [27] studied the thermo-mechanical responses of CAES systems in rock salt caverns. They concluded that the stability of the salt cavern is affected by the internal operating conditions with the air temperature increasing the volume convergence. Li et al. [28] carried out a failure study for gas storage by thermo-mechanical modelling in salt caverns. The results showed an affected area by the air and temperature operating conditions up to 10 m in the rock mass from the cavern walls. Round trip efficiencies between 70.5 and 75.1% were estimated by Barbour et al. [29] in A-CAES systems for continuous operation with packed bed thermal energy storage. The thermodynamic and geomechanical performance was studied by Rutqvist et al. [30] in lined rock caverns. They found that 97% of the thermal energy injected in the charge period could be recovered in the discharge period. Deng et al. [31] designed a new CAES systems with constant gas pressure and temperature in closed coal mines. A power output of 18 MW and a generating time of 1.76 h were obtained for a mine tunnel with 10,000 m<sup>3</sup> of volume at 500 m depth. The ADELE project was studied

in Germany to install an A-CAES plant with a storage capacity of 360 MWh and output power of 90 MW [2]. was studied in Germany to install an A-CAES plant with a storage capacity of 360 MWh and output power of 90 MW [2]. In this paper, abandoned mines are proposed as underground reservoirs for large

closed coal mines. A power output of 18 MW and a generating time of 1.76 h were obtained for a mine tunnel with 10,000 m3 of volume at 500 m depth. The ADELE project

*Appl. Sci.* **2021**, *11*, x FOR PEER REVIEW 3 of 19

In this paper, abandoned mines are proposed as underground reservoirs for large scale energy storage systems. A 200 m<sup>3</sup> tunnel in an abandoned coal mine was investigated as compressed air reservoir for A-CAES plants, where the ambient air is stored at high pressure. The thermodynamic response of A-CAES reservoirs was analyzed considering three solids around the pressurized air: a 20 mm thick sealing layer, a 35 cm thick concrete lining, and a 2.5 m thick sandstone rock mass. One-dimensional analytical and threedimensional CFD numerical models were conducted considering fiber-reinforced plastic (FRP) and steel as sealing layers and a typical operational pressure from 5 to 8 MPa. In the analytical model, air mass flow rates of 0.22 and <sup>−</sup>0.45 kg s−<sup>1</sup> have been considered for 30 cycles of compression and decompression, respectively. To reduce the computational time, the numerical simulation has been conducted for five cycles with air mass flow rates of 50 and <sup>−</sup>75 kg s−<sup>1</sup> . The air temperature and pressure fluctuations were estimated in the simulations. The temperature and the heat flux were also analyzed on the contact surfaces of the solids considering both FRP and steel sealing layers and the heat convection from the air to the sealing layers was calculated. Finally, the energy balance through the sealing layer was obtained during air charging and discharging. scale energy storage systems. A 200 m3 tunnel in an abandoned coal mine was investigated as compressed air reservoir for A-CAES plants, where the ambient air is stored at high pressure. The thermodynamic response of A-CAES reservoirs was analyzed considering three solids around the pressurized air: a 20 mm thick sealing layer, a 35 cm thick concrete lining, and a 2.5 m thick sandstone rock mass. One-dimensional analytical and three-dimensional CFD numerical models were conducted considering fiber-reinforced plastic (FRP) and steel as sealing layers and a typical operational pressure from 5 to 8 MPa. In the analytical model, air mass flow rates of 0.22 and −0.45 kg s−1 have been considered for 30 cycles of compression and decompression, respectively. To reduce the computational time, the numerical simulation has been conducted for five cycles with air mass flow rates of 50 and −75 kg s−1. The air temperature and pressure fluctuations were estimated in the simulations. The temperature and the heat flux were also analyzed on the contact surfaces of the solids considering both FRP and steel sealing layers and the heat convection from the air to the sealing layers was calculated. Finally, the energy balance through the sealing layer was obtained during air charging and discharging.

### **2. Materials and Methods 2. Materials and Methods**

### *2.1. Problem Statement 2.1. Problem Statement*

Underground space in abandoned mines may be used as compressed air storage systems for CAES plants. The simplified schematic diagram of the CAES system is shown in Figure 1. The compressor and turbine facilities are installed above the ground, while the compressed air reservoir is underground. The ambient air is compressed during off-peak periods and stored at high pressure in the underground reservoir (charge period). Then, when the electricity is required, the compressed air is released, heated, and expanded during peak hours in conventional gas turbines, driving a generator for power production (discharge period). The charge and discharge processes are carried out with a typical operating pressure range from 5 to 8 MPa. To reduce air leakage, two different sealing layers, FRP, and steel, have been employed in the present study. The air temperature and pressure fluctuations are estimated for both FRP and steel sealing layers. Underground space in abandoned mines may be used as compressed air storage systems for CAES plants. The simplified schematic diagram of the CAES system is shown in Figure 1. The compressor and turbine facilities are installed above the ground, while the compressed air reservoir is underground. The ambient air is compressed during off-peak periods and stored at high pressure in the underground reservoir (charge period). Then, when the electricity is required, the compressed air is released, heated, and expanded during peak hours in conventional gas turbines, driving a generator for power production (discharge period). The charge and discharge processes are carried out with a typical operating pressure range from 5 to 8 MPa. To reduce air leakage, two different sealing layers, FRP, and steel, have been employed in the present study. The air temperature and pressure fluctuations are estimated for both FRP and steel sealing layers.

**Figure 1.** Schematic diagram of compressed air energy storage (CAES) system in abandoned underground mines. Compressor and turbine facilities installed on the surface and underground compressed air reservoir with an operating pressure range from 5 to 8 MPa. **Figure 1.** Schematic diagram of compressed air energy storage (CAES) system in abandoned underground mines. Compressor and turbine facilities installed on the surface and underground compressed air reservoir with an operating pressure range from 5 to 8 MPa.

### *2.2. Analytical Model 2.2. Analytical Model*

A one-dimensional analytical model has been developed in MATLAB to study the thermodynamic performance during the operation time of the CAES system in a closed mine. Figure 2 shows the scheme of the 50 m long model. A lined tunnel with a usable volume of 200 m<sup>3</sup> and a cross section of 8 m<sup>2</sup> has been considered. A 20 mm thick sealing layer, a 35 cm thick concrete lining and a 2.5 m thick sandstone rock mass have also been considered in the model around the fluid. The evolution of the temperature (*Ta*), density (*ρa*) and air pressure (*Pa*) over time has been analyzed considering different air mass flow rates (*m˙* ). The heat flux through the contact surfaces, i.e., air-sealing layer ( . *Q*<sup>1</sup> ), sealing layer-concrete ( . *Q*<sup>2</sup> ), and concrete-sandstone ( . *Q*<sup>3</sup> ), has been estimated. Finally, the temperature on the sealing layer wall (*T*1), concrete wall (*T*2), and sandstone wall (*T*3) have also been estimated during the operation time considering an external temperature (*T*4) of 300 K. A one-dimensional analytical model has been developed in MATLAB to study the thermodynamic performance during the operation time of the CAES system in a closed mine. Figure 2 shows the scheme of the 50 m long model. A lined tunnel with a usable volume of 200 m3 and a cross section of 8 m2 has been considered. A 20 mm thick sealing layer, a 35 cm thick concrete lining and a 2.5 m thick sandstone rock mass have also been considered in the model around the fluid. The evolution of the temperature (*Ta*), density (*ρa*) and air pressure (*Pa*) over time has been analyzed considering different air mass flow rates (*ṁ*). The heat flux through the contact surfaces, i.e., air-sealing layer (ܳሶ ଵ), sealing layer-concrete (ܳሶ ଶ), and concrete-sandstone (ܳሶ ଷ), has been estimated. Finally, the temperature on the sealing layer wall (*T*1), concrete wall (*T*2), and sandstone wall (*T*3) have also been estimated during the operation time considering an external temperature (*T*4) of 300 K.

**Figure 2.** Scheme of the deep tunnel with a sealing layer, concrete lining and sandstone rock mass. (**a**) 3D model (not to scale); (**b**) Cross section. **Figure 2.** Scheme of the deep tunnel with a sealing layer, concrete lining and sandstone rock mass. (**a**) 3D model (not to scale); (**b**) Cross section.

The energy equation has been applied following Equation (1), where ܳሶ is the net heat transfer, ܹ ሶ is the net work done, *e* is the specific energy, *ρ* is the density in kg m−3, *V* is the reservoir volume, *v* is the air velocity and *S* is the cross section of the tunnel. In the air domain, the first term of Equation (2) represents the heat convection from the air to the sealing layer, which depends on the film coefficient of heat transfer (*h*), in Wm−2 K−1, the sealing layer surface in contact with the fluid (*A*1) and the temperature difference between the air and the concrete wall (*Ta-T*1). The second term of Equation (2) corresponds to the energy variation within the reservoir where the air pressure increases, which depends on the air mass inside the reservoir (*m*), the specific heat at constant volume (*Cv*) and the air temperature (*Ta*). Finally, the third term represents the input and output of energy in the control volume due to the energy provided by the air jet in the form of heat (*T*0) and kinetic energy. The air mass inside the reservoir at an instant of time, *t*, depends on the initial air mass (*m*0) and the air mass flow rate (*ṁ*), and is obtained by applying Equation (3). The energy equation has been applied following Equation (1), where *Q* is the net heat transfer, . *W* is the net work done, *e* is the specific energy, *ρ* is the density in kg m−<sup>3</sup> , *V* is the reservoir volume, *v* is the air velocity and *S* is the cross section of the tunnel. In the air domain, the first term of Equation (2) represents the heat convection from the air to the sealing layer, which depends on the film coefficient of heat transfer (*h*), in Wm−<sup>2</sup> K−<sup>1</sup> , the sealing layer surface in contact with the fluid (*A*1) and the temperature difference between the air and the concrete wall (*Ta-T*1). The second term of Equation (2) corresponds to the energy variation within the reservoir where the air pressure increases, which depends on the air mass inside the reservoir (*m*), the specific heat at constant volume (*Cv*) and the air temperature (*Ta*). Finally, the third term represents the input and output of energy in the control volume due to the energy provided by the air jet in the form of heat (*T*0) and kinetic energy. The air mass inside the reservoir at an instant of time, *t*, depends on the initial air mass (*m*0) and the air mass flow rate (*m˙* ), and is obtained by applying Equation (3).

$$
\dot{Q} - \dot{W} = \frac{\partial}{\partial t} \int \rho e \, dV + \oint \rho \mathcal{E} \, (\mathfrak{d}dS) \tag{1}
$$

.

$$-hA\_1(T\_a - T\_1) = \frac{d}{dt}(m\mathbb{C}\_{\mathbb{P}}T\_a) - \dot{m}\left(\mathbb{C}\_{\mathbb{P}}T\_0 + RT\_0 + \frac{v^2}{2}\right) \tag{2}$$

*m* = *m*<sup>0</sup> + . *m t* (3)

After some algebra, the air temperature (*Ta*) is calculated by applying Equation (4). The air temperature and density increase as the air pressure increases in the reservoir, and therefore the wall temperature increases by heat convection.

$$\frac{dT\_a}{dt} = \frac{\dot{m}}{m} \left[ \frac{\mathbb{C}\_p}{\mathbb{C}\_v} T\_0 + \frac{1}{2\mathbb{C}\_v} \left( \frac{\dot{m}}{m} L \right)^2 - T\_a \right] - \frac{h2\pi r\_1 L}{\mathbb{C}\_v m} (T\_a - T\_1) \tag{4}$$

The temperatures on the sealing layer wall (*T*1), the concrete wall (*T*2) and sandstone wall (*T*3) are obtained applying the equations of non-stationary heat transfer in sealing layer, concrete and sandstone domains (solid regions) considering an external temperature (*T*4) of 300 K and a temperature of the air mass flow inlet (*T*0) of 310 K. The heat flux reaches the sealing layer by convection. A part of this heat is transmitted to the concrete by conduction while another part increases the temperature of the sealing layer. A similar mechanism takes place for the concrete lining. Equations (5)–(7) represent the heat transmission on the sealing layer, concrete and sandstone formation, respectively.

$$\frac{dT\_1}{dt} = \frac{2\mathcal{U}\_\text{lt}}{m\_\text{s}\mathcal{C}\_{ps}}(T\_a - T\_1) - \frac{2\mathcal{U}\_\text{s}}{m\_\text{s}\mathcal{C}\_{ps}}(T\_1 - T\_2) - \frac{dT\_2}{dt} \tag{5}$$

$$\frac{dT\_2}{dt} = \frac{2\mathcal{U}\_\mathrm{s}}{m\_\mathrm{c}\mathcal{C}\_{pc}}(T\_1 - T\_2) - \frac{2\mathcal{U}\_\mathrm{c}}{m\_\mathrm{c}\mathcal{C}\_{pc}}(T\_2 - T\_3) \tag{6}$$

$$\frac{dT\_3}{dt} = \frac{2U\_\text{c}}{m\_\text{ss}\mathcal{C}\_{\text{pss}}}(T\_2 - T\_3) - \frac{2U\_\text{ss}}{m\_\text{ss}\mathcal{C}\_{\text{pss}}}(T\_3 - T\_4) \tag{7}$$

where *U<sup>h</sup>* , *Us*, *Uc*, and *Uss* are the convection transmittance, sealing layer transmittance, concrete transmittance, and sandstone transmittance, respectively. These thermal coefficients are evaluated according to

$$\mathbf{U}\_h = h 2 \pi r\_1 \mathbf{L} \tag{8}$$

$$\mathcal{U}\_s = \frac{2\pi \mathcal{K}\_s L}{Ln\left(\frac{r\_2}{r\_1}\right)}\tag{9}$$

$$\mathcal{U}\_{\mathbb{C}} = \frac{2\pi K\_{\mathbb{C}}L}{Ln\left(\frac{r\_2}{r\_2}\right)}\tag{10}$$

$$\mathcal{U}\_{\rm ss} = \frac{2\pi K\_{\rm ss} L}{L n\left(\frac{r\_4}{r\_3}\right)}\tag{11}$$

In these expressions, *h* is the film coefficient of heat transfer, in Wm−<sup>2</sup> K −1 ; *L* is the tunnel length, in m; *K<sup>s</sup> , K<sup>c</sup>* and *Kss* are the thermal conductivity of sealing layer, concrete and sandstone, in Wm−<sup>1</sup> K −1 ; *Cp<sup>s</sup> , Cpc*, and *Cpss* are the specific heat at constant pressure for the sealing layer, reinforced concrete and sandstone, in J kg−<sup>1</sup> K −1 ; and *r*1*, r*2*, r*3, and *r*<sup>4</sup> are the (equivalent) radius of sealing layer, concrete, sandstone and external walls, in m (Figure 2b). Note that due to the geometrical characteristics of the tunnel, the heat transfer formulation in cylindrical coordinates has been used for Equations (8)–(11).

For the evaluation of the convection film coefficient, the correlation proposed by Woodfield et al. [32] and Heath et al. [33] for the measurement of heat transfer coefficients in high-pressure vessels during gas charging was considered. In particular, a mixed (natural and forced) convection is modelled as a combination of both unsteady Reynolds and Rayleigh numbers. For the present investigation, a slight correction has been introduced for the exponent of the Rayleigh number, resulting the Nusselt number in the following expression:

$$Nu = 0.56 \, Re^{0.67} + 0.104 \, Ra^{0.34} \tag{12}$$

where *Re* is defined as the Reynolds number of the incoming mass flow rate (*Re* = 4 . *m*/2*µπr*1) and *Ra* is computed from the instantaneous thermal properties of the air (*Ra* = *gβ*(*T<sup>a</sup>* − *T*1)*Cpρ* 2*L* <sup>3</sup>/(*µK*)), being *g* the gravity acceleration, *β* is the volumetric thermal expansion coefficient, *Cp* the specific heat coefficient at constant pressure, *ρ* the density, *L* the tunnel length, *µ* the dynamic viscosity and *K* the thermal conductivity. The film coefficient is thus computed as

$$h = Nu \frac{K}{L} \tag{13}$$

Finally, a backward Euler explicit discretization has been employed over Equations (3)–(7) to resolve the coupled system of equations. A typical time-step of 0.01 s was applied to provide an accurate temporal description of the different temperatures. The model was resolved iteratively until a maximum prescribe pressure was reached inside the reservoir.

In addition, according to Kushnir et al. [18] and Zhou et al. [25], the air has been considered as a real gas through a compressibility factor (Z) that it is computed using the Berthelot gas state equation. Finally, considering a density value uniformly distributed in the volume, the pressure value of the CAES may be obtained at any instant by applying Equations (14) and (15).

$$P\_a = \left(\frac{R}{V}\right) Z m T\_a \tag{14}$$

$$Z = 1 - \frac{9}{128} \left(\frac{P\_a}{P\_c}\right) \left(\frac{T\_c}{T}\right) \left(\frac{6T\_c^2}{T^2} - 1\right) \tag{15}$$

where *T<sup>c</sup>* and *P<sup>c</sup>* the air temperature and the air pressure at critical conditions, assumed to be 132.65 K and 3.76 MPa. The contact surfaces between air-sealing layer (*A*<sup>1</sup> = 2*πr*1*L*), sealing layer-concrete (*A*2~2*πr*2*L*), concrete-sandstone (*A*3~2*πr*3*L*), and sandstone-exterior (*A*4~2*πr*4*L*), are 359, 365, 545, and 1610 m<sup>2</sup> , respectively. Note that, in order to model a complete cycle of compression (storing energy) and expansion (releasing energy) of the CAES system, the energy balance for the air domain—Equation (4)—needs to be rewritten as a function of the discharge mass flow rate (*m˙ out < 0*) and the air temperature within the reservoir, according to Equation (16):

$$\frac{dT\_a}{dt} = \frac{\dot{m}\_{out}}{m} \left[ \left( \frac{\mathbb{C}\_p}{\mathbb{C}\_v} - 1 \right) T\_a + \frac{1}{2\mathbb{C}\_v} \left( \frac{\dot{m}\_{out}}{m} L \right)^2 \right] - \frac{h2\pi r\_1 L}{\mathbb{C}\_v m} (T\_a - T\_1) \tag{16}$$

Air mass flow rates of 0.22 kg s−<sup>1</sup> and <sup>−</sup>0.45 kg s−<sup>1</sup> have been considered in the charge and discharge periods, respectively, for both RFP and steel sealing layers. The thermal properties and the volume of air, reinforced concrete, sandstone rock mass and sealing layers considered in the model are shown in Table 1 [24,25]. Note significant differences between the thermal conductivity of FRP and steel.

**Table 1.** Thermal properties and volume of air, concrete, sealing layers and sandstone rock mass.


### *2.3. CFD Numerical Model 2.3. CFD Numerical Model*

### 2.3.1. Model Geometry, Mesh and Boundary Conditions 2.3.1. Model Geometry, Mesh and Boundary Conditions

*Appl. Sci.* **2021**, *11*, x FOR PEER REVIEW 7 of 19

A 3D numerical model of a horseshoe-shaped tunnel located inside a closed mine was conducted to simulate the compression and decompression processes (air charging and discharging) that occur in the tunnel during the operation of the CAES plant. The computational domain is a tunnel with 8 m<sup>2</sup> of cross section and 50 m in length and includes both the fluid area and the solid areas around the fluid. The configuration of the simulated tunnel is illustrated in Figure 3. To improve the geomechanical performance, a 4 m<sup>2</sup> circular cross-section has been designed for the compressed air (blue area in Figure 3). The tunnel is reinforced with a 35 cm thick concrete lining and is finished with a dead-end. To avoid thermal leakage, a 20 mm thick sealing layer is considered between the air and concrete lining. Two types of sealing materials are considered: FRP and steel. The external part of the model corresponds to the 2.5 m thick sandstone rock mass existing in the mine. The total cross section of the model is 57.13 m<sup>2</sup> (8 <sup>×</sup> 8 m) and the model volume is 2587 m<sup>3</sup> . In addition, the useful volume of air reaches 200.53 m<sup>3</sup> . The entire geometry is meshed with 2,690,038 hexahedral and tetrahedral cells. Finer mesh was defined in the sealing layer, concrete and air zones, with higher grid density in those regions where the gradients of the flow characteristics are extremely important. A 3D numerical model of a horseshoe-shaped tunnel located inside a closed mine was conducted to simulate the compression and decompression processes (air charging and discharging) that occur in the tunnel during the operation of the CAES plant. The computational domain is a tunnel with 8 m2 of cross section and 50 m in length and includes both the fluid area and the solid areas around the fluid. The configuration of the simulated tunnel is illustrated in Figure 3. To improve the geomechanical performance, a 4 m2 circular cross-section has been designed for the compressed air (blue area in Figure 3). The tunnel is reinforced with a 35 cm thick concrete lining and is finished with a deadend. To avoid thermal leakage, a 20 mm thick sealing layer is considered between the air and concrete lining. Two types of sealing materials are considered: FRP and steel. The external part of the model corresponds to the 2.5 m thick sandstone rock mass existing in the mine. The total cross section of the model is 57.13 m2 (8 × 8 m) and the model volume is 2587 m3. In addition, the useful volume of air reaches 200.53 m3. The entire geometry is meshed with 2,690,038 hexahedral and tetrahedral cells. Finer mesh was defined in the sealing layer, concrete and air zones, with higher grid density in those regions where the gradients of the flow characteristics are extremely important.

Sealing layer (Steel) 500 45 7.24 76.5

**Figure 3.** Computational domain. Mesh and boundary conditions used in CFD simulations. **Figure 3.** Computational domain. Mesh and boundary conditions used in CFD simulations.

Four different type of materials (air, sealing layer, concrete, and sandstone) are used to simulate the heat transfer process between the air inside the mine and the surrounding media. In addition, as in the analytical model, FRP and steel are considered as sealing layers. The air is defined as real-gas to allow the simulation of the compression process. To reproduce the concrete layer existing between the sealing layer and the sandstone, a solid material zone is defined. The sandstone and sealing layer zones are also defined as Four different type of materials (air, sealing layer, concrete, and sandstone) are used to simulate the heat transfer process between the air inside the mine and the surrounding media. In addition, as in the analytical model, FRP and steel are considered as sealing layers. The air is defined as real-gas to allow the simulation of the compression process. To reproduce the concrete layer existing between the sealing layer and the sandstone, a solid material zone is defined. The sandstone and sealing layer zones are also defined as solid materials. In addition, a series of thermal properties are imposed on the sandstone, concrete and sealing zones to simulate the conduction heat transfer.

The numerical simulation of the compression/expansion process was carried out with the commercial software ANSYS Fluent® 16.2. This code was used to solve the full 3-D Unsteady Reynolds-Averaged Navier–Stokes (URANS) equations for compressible flow. The coupling between velocity and pressure was resolved by the SIMPLE algorithm. The spatial and temporal derivatives of the governing equations for the fluid flow were calculated by means of a second-order discretization, and the pressure interpolation PREssure STaggering Option (PRESTO) scheme has been used. Turbulent closure was established by the RNG k–ε model together with standard wall functions. In addition, the resolution of the energy equation for both fluid and solid volumes was activated. A constant mass flow rate was imposed as an inlet boundary condition for each model and an interface condition was established for solid/solid and solid/fluid surfaces. To reduce the computational time, air mass flow rates of 50 kg s−<sup>1</sup> and <sup>−</sup>75 kg s−<sup>1</sup> were considered in the numerical simulations. In Section 3.3, these mass flow rates are also employed in the analytical model to carry out the comparative analysis. In addition, a constant temperature of 300 K was considered in the sandstone surface. Moreover, a no-slip shear condition was imposed on the walls. The material properties used in the numerical model have been indicated in Table 1. A fixed time step of 5 <sup>×</sup> <sup>10</sup>−<sup>2</sup> s was selected for the simulations to ensure their stable convergence. The convergence criteria is based on the residual values of the solution solved in the numerical domain. A typical threshold of 10−<sup>5</sup> was set for continuity, momentum and energy equations.
