# A Coupled Thermo-Hydro-Mechanical Model of Jointed Hard Rock for Compressed Air Energy Storage

Authors: Xiaoying Zhuang; Runqiu Huang; Chao Liang; Timon Rabczuk
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

Renewable energy resources such as wind and solar are intermittent, which causes instability when being connected to utility grid of electricity. Compressed air energy storage (CAES) provides an economic and technical viable solution to this problem by utilizing subsurface rock cavern to store the electricity generated by renewable energy in the form of compressed air. Though CAES has been used for over three decades, it is only restricted to salt rock or aquifers for air tightness reason. In this paper, the technical feasibility of utilizing hard rock for CAES is investigated by using a coupled thermo-hydro-mechanical (THM) modelling of nonisothermal gas flow. Governing equations are derived from the rules of energy balance, mass balance, and static equilibrium. Cyclic volumetric mass source and heat source models are applied to simulate the gas injection and production. Evaluation is carried out for intact rock and rock with discrete crack, respectively. In both cases, the heat and pressure losses using air mass control and supplementary air injection are compared.


---

## Body

## 1. Introduction

Renewable energy such as wind, solar, tidal, and wave only produces electricity intermittently and with low power and energy density, thus, nondispatchable and difficult to use at large scales as the modern society requires [1]. That is why many renewable energy technologies are lacking the economies of scale, which reduces their competitiveness and delays the transition to a low carbon economy. Therefore, economic solutions to bulk energy storage are urgently needed in order for renewable energy to take a significant share in the total energy mix. A critical issue for renewable energy to be integrated into grids with satisfactory stability is appropriate energy storage to defer electricity demand from peak to off peak times.Most energy storage systems are expensive, either in terms of Capex and Opex or in terms of energy losses incurred in storing and retrieving the energy. For example, batteries are costly, fly wheels are suitable for short-duration storage only. The CAES, besides pumped-hydro, is the only conceivable technology able to provide the very large scale energy storage deliverability above 100 MW in single unit sizes while free from adverse environmental effects of pumped-hydro. Hence, CAES has recently received lots of attention [2, 3] and it has been recently proposed that large-scale solar-CAES and wind-CAES deployment can enable renewable energy to compete against coal-fired electricity generation [4, 5]. In CAES, a source energy is stored in the form of highly pressurized air in underground rock caverns and the compressed air is released through turbines to generate electricity when needed [6] as shown in Figure 1. Underground rock caverns are the mostly chosen type of reservoirs in the CAES which provide initial ground stress against air pressure, strengths to withstand cyclic loadings, and large volumes of storage space and therefore are more economical than other types of storage such as steel tanks [7].Figure 1
Compressed air energy storage plant using hard rock caverns.Existing and planning CAES facilities for utility purpose are mainly in salt rock formation where vacuums are formed by pumping pressurized water to dissolve the rock medium. The only existing CAES built in hard rock, however, for experimental purpose, is a 1600 m3 cave constructed in deplete mine in Hokkaido, Japan. The air pressure during cyclic charging and discharging is between 4 and 8 MPa and rubber lining is used for air tightness. Another small scale CAES using water sealing in hard rock is now under plan in Switzerland and air pressure will be around 3.3 MPa. Though the CAES technology has come into existence for over 30 years [8], the exploration of utilizing hard rock is important for wider application of this technique since the requirement of salt rock is too restrictive geological constraint in many areas. In this paper, we carried out the investigation of utilizing hard rock cavern for CAES through a coupled thermal-hydro-mechanical modelling for physical quantities.


## 2. Geological Conditions for CAES

To date, there have been only two CAES plants in operation at utility scales. The first is the 290 MW Huntorf plant in Germany with a total storage volume of 310,000 m3, operational since 1978 [9, 10]. And the second of 110 MW-scale with a storage space of 500,000 m3 was built in McIntosh, AL, in 1991 [6]. While energy conversion equipment including compressors, recuperators, and gas turbines is well established technology and procurable from the market, present challenges of CAES lie in finding geo-storage spaces such as rock caverns that are competent for cost effective and reliable performance during long-term plant operation. Though using subsurface rock as oil/gas storage started from the 1900s, new issues arise for CAES including the cyclic actions of high air pressure and thermalfluidics. Other types of storage such as above ground gas tanks, though technically available, are too expensive, about 5 times the cost of caverns [11]. Moreover, the security of the storage is also a major concern for highly pressurized tanks above ground. Hence, the implementation of CAES technology waits for the exploring and characterizing sites suitable for storage.To the best of the authors’ knowledge, the CAES projects declared over the world are listed in Table1. From the successful experiences of the existing facilities, the caverns are ideally constructed in salt rock formations due to their good air tightness. However, the requirements of salt rocks are too restrictive and not available in many areas. Recently, research interests in exploring the potential of using bedrocks other than salt rock have been motivated such as crystalline rock [12], sedimentary porous rock [13], aquifer, depleted oil/gas wells [14], and limestone [15]. More CAES plants built on various types of rock are under plan around the world such as the 300 MW plant by Pacific Gas & Electric in CA, USA [15], 300 MW plant in China [11], and other pilot plants in Korea [12] and Israel [13]. In Australia, sandstones are typical hard rock from depleted copper, gold, and silver mines. Sandstones are normally porous enough to store large quantities and can be used as valuable aquifers for reservoirs. Such CAES can provide sufficient strengths to sustain the stress induced by the rise and drop of air pressure and temperatures. The depleted reservoir formation is economically attractive because it allows the reuse, with suitable modification, which reduces the start-up costs.Table 1
CAES projects over the world (data collected from various sources).
Project name
Location
Power capacity (MW)
Geological formation
Type
Cavern volume (m3)
Embedded depth
Max. pressure (MPa)
Status
Huntorf
Breman, Germany
290
Salt Rock
Diabetic
310,000
600
4.3–7.0
Operation
Norton
Norton,USA
2700
Hard rock (limestone)
Diabatic
9,600,000
670
5.5–11
Construction
Iowa Energy Park
Iowa,USA
270
Porous sand stone
Diabatic
—
914
—
Construction
ADELE
Germany
300
Salt Rock
Adiabatic
—
—
~10
Planning
Matagorga
Texsas,USA
540
Salt Rock
Diabatic
—
—
—
Planning
Seneca
New York,USA
150–270
Salt Rock
Diabatic
150,000
760
8–11
Planning
PG&E
California,USA
300
Porous aquifer
Diabatic
—
—
—
Planning
Datang CAES
Inner Mongolia,China
300
Sand stone
Diabatic
900,000
500 m
5–8
PlanningThe most significant challenge in CAES is to ensure the air tightness and prevent the air leak-off that may occur due to crack propagation under long-term cyclic actions resulting in pressure and heat losses. On the other hand, rock is natural mineral aggregate with distributed discontinuities such as joints, cracks, or inclusions. The CAES operation exerts complex and cyclic mechanical, thermal, and wet-dry cycles. These excitations may cause monotonic failure and require close examination useful for engineering design of cavern. The stagnation of the 2700 MW CAES plant in Huston, the US, is an example [15]. The construction of the plant was launched in 2001 and stalled later due to the air leakage problem in the limestone caverns occurred [11, 14]. The understanding of the air leak-off requires both understanding on the mechanisms, geological conditions, the development of physical model, and advanced numerical simulation for engineering design [16]. In the past two decades, there have been significant advancements in fracture simulation methods such as the Extended Finite Element Method (XFEM) or extended meshfree methods [17–36] to allow for the application of discrete crack modelling techniques to THM simulations [37, 38] as well as fluid structure interaction-induced cracks [39]. In this paper, we aim to find out the relation between cavern characteristics, energy efficiency, and heat and pressure losses and stresses of rock induced by hydraulic and thermal changes and the existence of crack.


## 3. Governing Equations of the THM Modelling

In this paper, three fields, namely, mechanical, thermal, and hydro fields, are considered and coupled. The hydro or seepage field represents the gas flow behaviour, which changes the temperature distribution by the thermalfluidic flow through the porous media. The thermal field in turn will change the density and pressure of the fluid and hence the velocity of the fluid. The thermal strain will result in the internal stress of the mechanical field, which influences the seepage field through volumetric strain. The seepage will change the pore pressure of the fluid with the porous media. The complete coupling of the three fields can involve a large quantity of field variables and material properties associated with each other in the three fields, which is computationally expensive and unnecessary from an engineering perspective. For example, the viscosity of the fluid is temperature dependent and can add another relation between seepage and thermal flied. Therefore, we introduce some assumptions in the present THM modelling framework as follows.(1)
The rock is treated as isotropic elastic porous media and the mechanical parameters of the rock are not changing with respect to the temperature.
(2)
The seepage of the air through the rock follows Darcy’s law where the inertia and gravity of the air are not considered.
(3)
The physical state of compressed air conforms with the ideal gas law.
(4)
There is no ground water considered in the model and it is a single phase fluid flow problem, that is, air. The rock pore is saturated by air.
(5)
The heat mass is transferred between the solid and fluid through both convection and conduction. The heat capacity and conductivity coefficients of all the media are not changing with respect to the temperature and pressure. The volumetric strain energy of the solid and the work produced by the fluid from the air compression are taken into account for the energy balance.
(6)
The volume change of the porous media is mainly due to the deformation of the pores.The above assumptions are used both in the derivation of the governing equations in the following, and the model settings of the numerical studies in Section4.


### 3.1. Fluid Field: Mass Conservation

The seepage of the air in rock satisfies the continuity equation described by(1)
∂
(
ρ
gas
ϕ
)
∂
t
+
∇
·
(
ρ
gas
q
)
=
Q
m
,
where ρ
gas is the density of the air, ϕ is the pore ratio of rock, q is the seepage velocity vector of the fluid, Q
m is the source of mass (positive for inputting fluid and negative for outputting). The first term on the LHS of (1) can be further calculated by
(2)
∂
(
ρ
gas
ϕ
)
∂
t
=
ϕ
∂
ρ
gas
∂
t
+
ρ
gas
∂
ϕ
∂
t
=
ϕ
∂
ρ
gas
∂
t
+
b
·
ρ
gas
∂
ε
V
∂
t
,
where ε
V is the volumetric strain and b is Biot’s consolidation coefficient. The last term on the RHS of (2) shows the direction relation between the volumetric change of the porous media and the porosity ratio. The density of the air is changing with the temperature T and pressure P which is described by the ideal gas law
(3)
ρ
gas
=
P
M
R
T
,
where P is the pore pressure, M is Molar mass of air that M
=
29
g
/
mol, R is ideal gas constant that R
=
8.314
J
/
(
mol
·
K
), T is absolute (thermodynamic) temperature and P is the pore pressure. The momentum of air seepage conforms to Darcy’s law. If the inertia, gravity, and adhesion of air is not considered, the seepage equation of air is expressed by
(4)
q
=
-
k
η
∇
P
,
where k is the permeability coefficient of rock and η is the coefficient of viscosity.By substituting (3) and (4) into (2), the governing equations of seepage field can be obtained by
(5)
b
·
ρ
gas
∂
ε
V
∂
t
+
ϕ
M
R
T
∂
P
∂
t
-
ϕ
M
P
R
T
2
∂
T
∂
t
-
∇
(
P
M
R
T
k
η
∇
P
)
=
Q
m
.The changing density is the main difference between the governing equations of air seepage and impressible fluid seepage. The density of air changes with the temperature, and therefore the thermal field is associated with the fluid field. It can be seen from (5) that the temperature T appears in the coefficient term ϕ
M
/
R
T for the pressure temporal derivatives ∂
P
/
∂
t and also adjoins with the spatial gradient of pressure term ∇
P. This results in the high nonlinearity of the governing equation which brings numerical convergence and stability problem. It can also be seen from the first term on LHS that volumetric strain of the solid also couples the stress field with fluid field.


### 3.2. Thermal Field: Energy Conservation

According to the first law of thermodynamics, the heat added to a system should be equal to the change of the internal energy of the system and the work done in the system. By considering the energy conservation of the skeleton of the rock as porous media, we can get(6)
ρ
s
C
p
,
s
∂
T
∂
t
+
T
K
s
β
∂
ε
V
∂
t
=
∇
·
(
α
s
∇
T
)
+
Q
s
,
where ρ
s, C
p
,
s, α
s, Q
s, K
s, and β are the density, specific heat, thermal conductivity, heat source, volumetric modulus, and thermal expansion coefficient of the porous media. It can be seen that the energy transportation in the solid is through the thermal conduction as the first term on RHS. The volumetric work induced by the temperature as the second term on LHS also results in the dissipation of energy. From the energy conservation of fluid, we can get
(7)
ρ
gas
C
p
,
gas
∂
T
∂
t
+
ρ
gas
C
p
,
gas
q
·
∇
T
+
P
(
∇
·
q
)
=
∇
·
(
α
gas
∇
T
)
+
Q
gas
,
where ρ
gas, C
p
,
gas, α
gas, and Q
gas are the density of air, heat specific of air, thermal conductivity, and heat source of the fluid. The energy transportation of air includes the heat conduction as the first term on RHS, convection as the second term on LHS, and external work of air as the third term on LHS. For a single phase flow, we assume that the solid and air are in thermal steady state, that is, the temperature of the fluid and solid are equal at the interface. The heat exchange along the interface of two phases is not considered. By combining (6) with (7), we can get the following:
(8)
(
ρ
C
p
)
eq
∂
T
∂
t
+
ρ
gas
C
p
,
gas
q
·
∇
T
+
(
1
-
ϕ
)
T
K
s
β
∂
ε
V
∂
t
+
P
(
∇
·
q
)
=
∇
·
(
α
eq
∇
T
)
+
Q
,
(9)
(
ρ
C
p
)
eq
=
ϕ
·
ρ
gas
C
p
,
gas
+
(
1
-
ϕ
)
·
ρ
s
C
p
,
s
,
(10)
α
eq
=
ϕ
·
α
gas
+
(
1
-
ϕ
)
·
α
s
,
(11)
Q
=
ϕ
·
Q
gas
+
(
1
-
ϕ
)
·
Q
s
,
where in (8) (
ρ
C
p
)
eq is the equivalent specific heat capacity of rock as a porous media saturated with air, α
eq is the equivalent heat conduction coefficient, and Q is the equivalent heat source of saturated porous media. Equations (8)–(11) describe the governing equations of porous media and the fluid as the thermal field governing equations.


### 3.3. Mechanical Field: Equilibrium Condition

Since the rock is assumed to be saturated with air, the constitutive model of thermal porous media should consider the contribution of the pore pressure and temperature of the air as follows:(12)
σ
=
(
λ
+
μ
)
(
∇
·
u
)
+
μ
∇
u
-
b
P
I
-
(
3
λ
+
2
μ
)
β
(
T
-
T
0
)
I
,
where σ is the Cauchy stress tensor, λ and μ are Lamé coefficients of rock, u is the displacement vector, I is a second-order identity tensor, β is linear thermal expansion coefficient of rock, and b is Biot’s consolidation coefficient. The relation between E and υ with respect to λ and μ is well known as
(13)
λ
=
E
υ
(
1
+
υ
)
(
1
-
2
υ
)
,
μ
=
E
2
(
1
+
υ
)
.The equilibrium of porous media is expressed as(14)
∇
σ
+
f
=
0
,
where f is the body force vector. Substituting (12) into (14) leads to the governing equation of the stress field as
(15)
(
λ
+
μ
)
∇
(
∇
·
u
)
+
μ
∇
2
u
-
b
·
∇
P
-
(
3
λ
+
2
μ
)
β
·
∇
T
+
f
=
0
.


### 3.4. The Governing Equation of Fluid Flow along Single Discrete Crack (Joint)

The existence of crack will have influence on all the three fields, especially the mechanics field and stress field. The thermal field will also change due to the crack; however, it will be relatively minor compared to stress and fluid field.Due to the existence of the crack in rock, the distribution of the air is changing in the rock and consequently results in the change of effective stress and thermal field. In the analysis, the joint is modeled as a channel, that is, internal boundaries in rock. We have not considered the influence of crack over the distribution of thermal and mechanical fields. Therefore, the stress field and thermal field over the crack are continuous. The coupling of the three fields for discontinuities is an interesting and challenging topic for the future. It is generally recognized that a continuum based equivalent elastic porous media, such as Barentblatt and Gil’man [40] nonequilibrium two-phase immiscible flow mode, is not capable of capturing the discontinuity or jump in the thermalfluidic and mechanical field. In this type of approach, the mass and energy is transferred with higher gradient compared to base rock by modifying the material parameters. The crack is modeled as equivalent porous media and there exists numerical instability of extra refined mesh near the crack. In the present paper, the discrete crack model is adopted for a single joint. The upper and lower surfaces of the crack are modeled as internal boundaries inside the domain and the fluid is allowed to flow between the boundaries. Darcy’s law for tangential flower along crack is defined as
(16)
q
f
=
-
k
f
η
d
f
(
∇
T
P
)
,
where q
f is the flow rate along the unit length of a crack, k
f is the permeability coefficient, d
f is the width of the crack opening, and ∇
T
P is the pressure gradient along the length of the crack. The continuity of the fluid, that is, the mass conservation, between the two surfaces of a single joint is described by
(17)
d
f
∂
(
ρ
gas
ϕ
f
)
∂
t
+
∇
T
·
(
ρ
gas
q
f
)
=
d
f
Q
m
.
It can be seen that the same variable of air pressure, namely, P, appears in the governing variable in both the governing equation of the fluid flow along the fracture (joint) and the rock as the porous media.


## 3.1. Fluid Field: Mass Conservation

The seepage of the air in rock satisfies the continuity equation described by(1)
∂
(
ρ
gas
ϕ
)
∂
t
+
∇
·
(
ρ
gas
q
)
=
Q
m
,
where ρ
gas is the density of the air, ϕ is the pore ratio of rock, q is the seepage velocity vector of the fluid, Q
m is the source of mass (positive for inputting fluid and negative for outputting). The first term on the LHS of (1) can be further calculated by
(2)
∂
(
ρ
gas
ϕ
)
∂
t
=
ϕ
∂
ρ
gas
∂
t
+
ρ
gas
∂
ϕ
∂
t
=
ϕ
∂
ρ
gas
∂
t
+
b
·
ρ
gas
∂
ε
V
∂
t
,
where ε
V is the volumetric strain and b is Biot’s consolidation coefficient. The last term on the RHS of (2) shows the direction relation between the volumetric change of the porous media and the porosity ratio. The density of the air is changing with the temperature T and pressure P which is described by the ideal gas law
(3)
ρ
gas
=
P
M
R
T
,
where P is the pore pressure, M is Molar mass of air that M
=
29
g
/
mol, R is ideal gas constant that R
=
8.314
J
/
(
mol
·
K
), T is absolute (thermodynamic) temperature and P is the pore pressure. The momentum of air seepage conforms to Darcy’s law. If the inertia, gravity, and adhesion of air is not considered, the seepage equation of air is expressed by
(4)
q
=
-
k
η
∇
P
,
where k is the permeability coefficient of rock and η is the coefficient of viscosity.By substituting (3) and (4) into (2), the governing equations of seepage field can be obtained by
(5)
b
·
ρ
gas
∂
ε
V
∂
t
+
ϕ
M
R
T
∂
P
∂
t
-
ϕ
M
P
R
T
2
∂
T
∂
t
-
∇
(
P
M
R
T
k
η
∇
P
)
=
Q
m
.The changing density is the main difference between the governing equations of air seepage and impressible fluid seepage. The density of air changes with the temperature, and therefore the thermal field is associated with the fluid field. It can be seen from (5) that the temperature T appears in the coefficient term ϕ
M
/
R
T for the pressure temporal derivatives ∂
P
/
∂
t and also adjoins with the spatial gradient of pressure term ∇
P. This results in the high nonlinearity of the governing equation which brings numerical convergence and stability problem. It can also be seen from the first term on LHS that volumetric strain of the solid also couples the stress field with fluid field.


## 3.2. Thermal Field: Energy Conservation

According to the first law of thermodynamics, the heat added to a system should be equal to the change of the internal energy of the system and the work done in the system. By considering the energy conservation of the skeleton of the rock as porous media, we can get(6)
ρ
s
C
p
,
s
∂
T
∂
t
+
T
K
s
β
∂
ε
V
∂
t
=
∇
·
(
α
s
∇
T
)
+
Q
s
,
where ρ
s, C
p
,
s, α
s, Q
s, K
s, and β are the density, specific heat, thermal conductivity, heat source, volumetric modulus, and thermal expansion coefficient of the porous media. It can be seen that the energy transportation in the solid is through the thermal conduction as the first term on RHS. The volumetric work induced by the temperature as the second term on LHS also results in the dissipation of energy. From the energy conservation of fluid, we can get
(7)
ρ
gas
C
p
,
gas
∂
T
∂
t
+
ρ
gas
C
p
,
gas
q
·
∇
T
+
P
(
∇
·
q
)
=
∇
·
(
α
gas
∇
T
)
+
Q
gas
,
where ρ
gas, C
p
,
gas, α
gas, and Q
gas are the density of air, heat specific of air, thermal conductivity, and heat source of the fluid. The energy transportation of air includes the heat conduction as the first term on RHS, convection as the second term on LHS, and external work of air as the third term on LHS. For a single phase flow, we assume that the solid and air are in thermal steady state, that is, the temperature of the fluid and solid are equal at the interface. The heat exchange along the interface of two phases is not considered. By combining (6) with (7), we can get the following:
(8)
(
ρ
C
p
)
eq
∂
T
∂
t
+
ρ
gas
C
p
,
gas
q
·
∇
T
+
(
1
-
ϕ
)
T
K
s
β
∂
ε
V
∂
t
+
P
(
∇
·
q
)
=
∇
·
(
α
eq
∇
T
)
+
Q
,
(9)
(
ρ
C
p
)
eq
=
ϕ
·
ρ
gas
C
p
,
gas
+
(
1
-
ϕ
)
·
ρ
s
C
p
,
s
,
(10)
α
eq
=
ϕ
·
α
gas
+
(
1
-
ϕ
)
·
α
s
,
(11)
Q
=
ϕ
·
Q
gas
+
(
1
-
ϕ
)
·
Q
s
,
where in (8) (
ρ
C
p
)
eq is the equivalent specific heat capacity of rock as a porous media saturated with air, α
eq is the equivalent heat conduction coefficient, and Q is the equivalent heat source of saturated porous media. Equations (8)–(11) describe the governing equations of porous media and the fluid as the thermal field governing equations.


## 3.3. Mechanical Field: Equilibrium Condition

Since the rock is assumed to be saturated with air, the constitutive model of thermal porous media should consider the contribution of the pore pressure and temperature of the air as follows:(12)
σ
=
(
λ
+
μ
)
(
∇
·
u
)
+
μ
∇
u
-
b
P
I
-
(
3
λ
+
2
μ
)
β
(
T
-
T
0
)
I
,
where σ is the Cauchy stress tensor, λ and μ are Lamé coefficients of rock, u is the displacement vector, I is a second-order identity tensor, β is linear thermal expansion coefficient of rock, and b is Biot’s consolidation coefficient. The relation between E and υ with respect to λ and μ is well known as
(13)
λ
=
E
υ
(
1
+
υ
)
(
1
-
2
υ
)
,
μ
=
E
2
(
1
+
υ
)
.The equilibrium of porous media is expressed as(14)
∇
σ
+
f
=
0
,
where f is the body force vector. Substituting (12) into (14) leads to the governing equation of the stress field as
(15)
(
λ
+
μ
)
∇
(
∇
·
u
)
+
μ
∇
2
u
-
b
·
∇
P
-
(
3
λ
+
2
μ
)
β
·
∇
T
+
f
=
0
.


## 3.4. The Governing Equation of Fluid Flow along Single Discrete Crack (Joint)

The existence of crack will have influence on all the three fields, especially the mechanics field and stress field. The thermal field will also change due to the crack; however, it will be relatively minor compared to stress and fluid field.Due to the existence of the crack in rock, the distribution of the air is changing in the rock and consequently results in the change of effective stress and thermal field. In the analysis, the joint is modeled as a channel, that is, internal boundaries in rock. We have not considered the influence of crack over the distribution of thermal and mechanical fields. Therefore, the stress field and thermal field over the crack are continuous. The coupling of the three fields for discontinuities is an interesting and challenging topic for the future. It is generally recognized that a continuum based equivalent elastic porous media, such as Barentblatt and Gil’man [40] nonequilibrium two-phase immiscible flow mode, is not capable of capturing the discontinuity or jump in the thermalfluidic and mechanical field. In this type of approach, the mass and energy is transferred with higher gradient compared to base rock by modifying the material parameters. The crack is modeled as equivalent porous media and there exists numerical instability of extra refined mesh near the crack. In the present paper, the discrete crack model is adopted for a single joint. The upper and lower surfaces of the crack are modeled as internal boundaries inside the domain and the fluid is allowed to flow between the boundaries. Darcy’s law for tangential flower along crack is defined as
(16)
q
f
=
-
k
f
η
d
f
(
∇
T
P
)
,
where q
f is the flow rate along the unit length of a crack, k
f is the permeability coefficient, d
f is the width of the crack opening, and ∇
T
P is the pressure gradient along the length of the crack. The continuity of the fluid, that is, the mass conservation, between the two surfaces of a single joint is described by
(17)
d
f
∂
(
ρ
gas
ϕ
f
)
∂
t
+
∇
T
·
(
ρ
gas
q
f
)
=
d
f
Q
m
.
It can be seen that the same variable of air pressure, namely, P, appears in the governing variable in both the governing equation of the fluid flow along the fracture (joint) and the rock as the porous media.


## 4. Numerical Modelling of the CAES System in Hard Rock

The governing equations described in Section3 for mass balance, energy conservation, and equilibrium are discretised by the finite element method and coded to study the hard rock cavern performance for CAES. In this section, a hard rock cavern under the designed operational conditions of CAES without considering the existence of ground water is modeled. The aim of the numerical modelling is to investigate the heat loss, pressure distribution, and seepage of air near the vicinity of cavern.A rock cavern with a diameter of 5 m and embedded with 100 m depth is modeled as shown in Figure2(a). Initial ground stress for the equivalent stress state is modeled and excavation step is modeled and the rock near the cavern is regarded as excavation disturbed zone. Standard atmospheric pressure and ground temperature of 286.15 K (13°C) are applied to the four edges of boundaries. Six points of interest marked P1–P6 are sampled in the simulation results as shown in Figure 2(b). The operational cycle of the cavern is 24 hours. In the initial condition of the cavern, the air is injected into the design value of 5.5 MPa after 16 hours injection and followed by 8 hours storage. After that, the air injection and release follows a cycle of 8 hours storage, 4 hours storage, 4 hours release, and another 8 hours storage. In air compression stage, the air is injected at the rate of 1.12
×
10
-
3 kg
/
(
s
·
m
3
) to the cavern at the temperate of 296.65 K (23.5°C). In the energy recuperation stage, the air is released at the rate of 2.24
×
10
-
3 kg
/
(
s
·
m
3
). The time length of the present example is 80 operational days. The material parameters used in this case study are listed in Table 2.Table 2
Material parameters and boundary conditions used in the example.
Parameters
Young's modulusE of rock (GPa)
35
Poisson's ration of rockυ
0.3
Density of rockρ
s (kg/m3)
2800
Pore ratioϕ
1.0
Permeability coefficientk (m2)
1
×
10
-
9
Permeability coefficient along crackk
f (m2)
—
Crack widthd
f (m)
—
Viscosity of airη (Pa·s)
1.86
×
10
-
5
Biot's consolidation coefficientb
0.95
Heat conduction of rockα
s (W/m·K)
3
Heat conduction of airα
s (W/m·K)
1000
Heat specific of air under constant pressureC
p
,
gas (J/kg·K)
1000
Heat specific of rockC
p
,
s (J/kg·K)
900
Expansion coefficientβ (1/K)
1.0
×
10
-
5
Initial pressureP
0 (atm)
1
Initial temperatureT
0 (K)
286.15
Air injection temperatureT
in (K)
296.65Model setting of the hard rock cavern. (a) Geometry and profile of cavern. (b) Points of interest near the cavern.
(a)
(b)The pressure drop along the cavern ceiling with respect to the time is shown in Figure3. Due to the seepage of air through the pores of rock, pressure drop is observed in the consecutive days of operation. With the mass conservation control, that is, same injection and production of air in terms of mass, a substantial drop of pressure can be observed after 60 days of operation. At day 60, the minimum pressure inside the cavern is about 5 MPa lower than the desired operational pressure. This is different from the observation of CAES cavern in aquifer by Kim et al. [41]. It means that for hard rock cavern without water concealing, supplementary air injection is necessary to maintain operational pressure due to air seepage.Figure 3
Pressure drop of air after 60 operational days.It was also found that the influence of pressure fluctuation over the rock diminishes at about 2.1 m distant to the internal boundary of cavern (at P3 in Figure2). It can also be seen that the air seepage through the rock gradually converges as shown in Figure 4. On day 1, the pressure drop along the distance to the cavern is steep. The maximum pressure occurs at the interface of air and rock. At day 20, the maximum pressure is observed at about 0.2 m distant to the cavern boundary. After 60 days of operation, the maximum pressure point is shifted to about 0.5 m to the cavern boundary. The indication of the result to cavern design is that the validation of rock strength should be conducted for certain range to the cavern boundary. The change of pressure and pore pressure at different distances to the cavern with respect to the time is plotted in Figure 5. It can be seen that there is slight lagging behind for the development of stress in rock with respect to the pore pressure. This conforms to common phenomenon of the fluid seepage and saturation in porous media, and shows the validity of the present THM models. It can also be seen that at P3 which is 2.1 m distant to the inner boundary of cavern, the pore pressure vanishes showing the vanishing influence from cavern.Figure 4
Pressure variation along the cavern ceiling from Day 1 to 60.Figure 5
The maximum principle stress and pore pressure with respect to time.For the thermal field, the temperature gradually converges to the ground temperature of 13°C when it is distant from the cavern as shown in Figure6. The temperature on Day 1 at steady state is about 3 degrees than the injection temperature due the thermal conduction. It can be seen that the temperature decreases since the expansion of air is a heat consumption process.Figure 6
The temperature along cavern ceiling with respect to the time from Day 1 to 60.To study the pressure control strategy of the system, supplementary air condition is also modeled. By comparing the results of air density at day 60 with respect to day 1, the average density loss can be calculated which is about 1.048 kg/m3 daily. Therefore, in the model, supplementary air mass of 3.252% is injected. The results of cavern pressure with and without supplementary injection are compared in Figure 7. It can be seen that with supplementary air injection, the air pressure can be maintained between 5.5 and 8.5 MPa, which is ideal for the operational requirement. The maximum principle stresses at different points in rock are compared between mass conservation and air supplement conditions in Figure 8. The air supplement on one hand results in more stable air pressure inside cavern but also induces the increased stress in rock.Figure 7
The pressure difference with and without considering the volumetric strain induced by thermal field.Figure 8
The comparison of the maximum principal stress in rock between air supplement control and mass control during CAES operation.The difference of pressure between intact rock cavern and a carven with a single joint at the rock ceiling is compared as shown in Figure9. It can be seen that a single fracture with settings in Table 2 results in the pressure loss of about 4 MPa after 600 hours operation. For the energy loss, the cracked rock cavern has lost about 5 degrees of temperature as shown in Figure 10. It can also be seen that, cracked rock cavern will have larger fluctuation compared to intact rock in each operational cycle. The pores pressure distribution shows the seepage of air through the joint after 25 days of operation where the concentration takes place near the joint as shown in Figure 11. The internal pressure of air through dissipation into rock induces the pore pressure where the maximum value takes place near the joint of 2.5 MPa. In the present study, the discontinuous stress field and seepage field are not considered here.Figure 9
Comparison of pressure loss between intact and jointed rock cavern.Figure 10
Comparison of heat loss between intact and jointed rock cavern.Comparison of the maximum principle stress in the rock cavern with and without joint on the ceiling. (a) Intact rock and (b) with single joint on the ceiling.
(a)
(b)


## 5. Discussions

In this paper, we investigate the feasibility of utilizing hard rock for compressed air energy storage by a coupled THM model. The energy loss, stress distribution, and pore pressure distribution during cyclic operational condition are studied. It was found that mass control based CAES operation results in energy loss and is not feasible in substantial hard rock cavern. Supplementary air injection is needed to maintain the required pressure level. For rock cavern with single joint, higher fluctuation of pressure and temperatures are found which is unfavorable for CAES system stability. The energy loss is significantly increased due to the crack and the principle. In the air injection modelling, the injection rate is based on the density loss of mass control condition after 60 days of operational condition. A constant air supplement of 3.25% is adopted. However, it is only an average means of injection while in the real seepage field and the supplementary injection should be dynamical to maintain designed pressure value. Only a single static crack is studied in the present work. The crack propagation due to cyclic THM actions is not considered. These two problems can be the topics for further study.



---

*Source: 179169-2014-01-19.xml*