7 Thermal-Hydraulic-Mechanical Coupling Analysis of Rock Mass Y. OHNISHI Kyoto University, Japan and AKIRA KOBAYASHI Hazama Corporation, Ibaraki, Japan
7.1
INTRODUCTION
7.2
THEORIES
7.2.1 7.2.2 7.2.3 7.2.4 7.2.5
191 193
General Remarks Equilibrium Equation Continuity Equation for Ground Water Energy Conservation Law Governing Equations and Numerical Technique
193 193 194 196 197
7.3
VERIFICATION O F THE C O D E
197
7.4
NUMERICAL EXAMPLES
198
7.4.1 7.4.2 7.4.3
Macro-permeability Test in the Stripa Project Buffer Mass Test in the Stripa Project An Hypothetical Deep Nuclear Waste Repository
198 201 203
7.5
SUMMARY AND CONCLUSIONS
206
7.6
NOMENCLATURE
207
7.7
REFERENCES
207
7.1
INTRODUCTION
In this chapter, methods for the analysis of coupled thermal, hydraulic and mechanical phe nomena are described. Coupled behavior in a rock mass is one of the most complicated phenomena to be studied in the field of rock mechanics. This subject has to be examined not only with experiments and measurements but also with powerful numerical analyses. M a n y coupled phe nomena have been observed and pairs of hydraulic, mechanical and thermal phenomena have been studied. However, for the projects planned in recent years, e.g. high-level radioactive waste disposal in deep geology, compressed air energy storage systems and super magnetic electric energy storage systems, the fully coupled hydraulic, mechanical and thermal behavior of the rock mass induced from the structure should be assessed. The prediction of the phenomena is most important for the plan and design of a structure constructed in a rock mass. F o r example, the behavior of the foundations of a d a m or the behavior of the rock mass around ground openings and tunnels is related to coupled hydraulic and mechanical phenomena. The discontinuities which exist in the rock mass, i.e. fractures or faults, have a signific ant effect on this behavior. F o r such a case, a model to explain the coupled hydraulic and mechanical behavior in a fractured rock mass is needed for prediction. 191
192
Rock Mechanics Continuum
Modeling
F o r a long time, solution techniques for these coupled problems lagged behind analytical solutions, and were restricted to simple geometries. The recent development of computers has greatly advanced numerical techniques in the field of coupled analysis. F o r coupled hydraulic and mechanical problems of saturated porous media, the consolidation problem of clay has been studied numerically for quite some time, and notable among these are studies by Christian and Boehmer [1] and Sandhu and Wilson [ 2 ] . Recently, a model considering the elasto-plasticity or visco-plasticity of soils has been developed. However, many problems related to hydraulic and mechanical analyses of fractured rock still exist. One of them is the modeling of a fractured medium in which the deformation and the ground water flow are strongly influenced by the geometry of the fractures. Moreover, it is not well understood whether or not the principle of effective stress can be applied for a fractured rock mass. The state of the art studies on the various coupled processes for rock joints have been summarized by Tsang [3]. In such a situation, it seems important to examine if the model used for the analysis is suitable for the object being evaluated and to understand the limitations and assumptions present in the model. Based on a thorough under standing of the model, we have to consider the results of the analyses. Currently, there exist a few approaches for fractured porous media. One is the discrete fracture approach, which regards the fractured medium as a continuum in which discrete discontinuities exist [4, 5]. In addition, the distinct element method (DEM), which considers the fractured medium as an assembly of blocks, has been used for many projects [6, 7]. The other approach is the equivalent continuum model, which models the anisotropic properties produced by existing fractures [8]. F o r coupled thermal and mechanical problems, there are many examples in which classical continuum mechanics are applied with numerical calculation techniques. In such cases, the medium is often considered as elastic. The influence of the mechanical behavior on energy conservation is often ignored in the analyses. In some cases, thermal stress is imposed by using an initial stress and is obtained separately before solving the mechanical problems [9]. F o r fractured rock masses, D E M has also been used for assessing the coupled thermal and mechanical behavior of the natural barrier of a radioactive waste repository [10]. When plastic deformation due to thermal stress has to be considered, the yield function is needed as a function of temperature in addition to stress. In many cases in the geotechnical field, plastic deformation is conventionally considered apart from the thermal stress [ 1 1 , 1 2 ] . There are many engineering problems related to coupled hydraulic and thermal behavior. Geothermal and freezing soil problems are typical problems in the geotechnical field. Vaporization is considered as a geothermal problem because of the high temperature and high pressure in the ground [ 1 3 , 1 4 ] . Frost action is intensively studied as a freezing soil problem. Such phase change behavior is very complicated and involves many unknown mechanisms. In order to analyze a complex phenomenon such as vaporization or freezing, the limitations and assumptions of the numerical model must be well understood before examining the numerical results. Coupled hydraulic and thermal problems related to high-level radioactive waste disposal in deep geology has also been intensively studied recently. The effect of increasing temperatures on the ground water flow around a repository, due to radioactive decay of the waste, has been studied by many researchers [ 1 5 , 1 6 ] . The phenomena examined in their analyses are mainly dependent on the ground conditions in which the repository is constructed. In cases where the repository may be constructed in an unsaturated region (e.g. Apache Leap tuff site in the U.S.A.), vaporization of the ground water is the most important phenomenon to be analyzed. O n the contrary, when a reposi tory is constructed in saturated deep crystalline rock, convection of the ground water induced by increasing temperatures is the most important matter for assessment. In more recent studies, the problems of coupled hydraulic, mechanical and thermal phenomena have been explored. Bear and Corapcioglu [17] derived basic equations which describe the thermoelastic behavior of the ground due to hot water injection into confined and leaky aquifers. Hart [18] presented a model which fully describes coupled thermal, mechanical and hydraulic behavior in nonlinear porous geological systems and which calculates the model by an explicit finite difference method. Noorishad et al. [19] applied a finite element method with joint elements to fully coupled phenomena for a saturated fractured rock mass. Ohnishi et al. [20] developed a finite element code to handle the problems of coupled hydraulic, thermal and mechanical behavior of a saturated-unsaturated geological medium. In this chapter, T H A M E S (2D and 3D Thermal, Hydraulic And MEchanical System analysis computer code), which has been developed to examine coupled hydraulic, thermal and mechanical phenomena of rock mass, is introduced. Firstly, the governing equations used in the code are explained. Secondly, the verification of the code is carried out by comparison with the available analytical and experimental results. Lastly, a few examples are solved with the code. Macro-
Thermal-Hydraulic-Mechanical
Coupling Analysis of Rock Mass
193
permeability tests and Buffer mass tests conducted in the Stripa project are examined and the numerical results compared with observations, and an imaginary deep nuclear waste depository is assessed and the results between the two-dimensional and three-dimensional analysis results are compared.
7.2 7.2.1
THEORIES General Remarks
The coupled phenomena considered in the T H A M E S code are shown in Figure 1. Heat transport, water flow and deformation behavior are related to each other in a fully coupled manner. The governing equations are constructed in such a way as to satisfy these relations. The governing equations are derived under the following assumptions, (i) The medium is porous and elastic, (ii) Darcy's law is valid for flow of water in a saturated-unsaturated medium, (hi) Heat transport occurs only in a solid and liquid phase. Heat transferred among three phases (solid, liquid and gas) is neglected, (iv) Phase change of water between liquid and vapor is not considered, (v) Fourier's law holds for heat flux, (vi) Water density varies depending upon the temperature and pressure of water. The way to introduce these assumptions into the governing equations is described in the following sections.
7.2.2
Equilibrium Equation
In general, the rock mass has many discontinuities and presents anisotropic and heterogeneous properties. Although methods to express the discontinuities in the ground, e.g. fractures and cracks, have been studied by many researchers as mentioned above, an authoritative model has not been established. The continuity model is one of the models, in which the anisotropic properties induced by the existing fractures are introduced into the constitutive law in various ways. In this chapter, the ground is considered as a continuous isotropic elastic body because of the simplicity and conveni ence for the model of complex coupled phenomena. Consider a volume of a n elastic porous medium filled with a homogeneous fluid. The distribution of total stresses within this volume can be shown to satisfy the equilibrium equation (1)
where is the total stress, p is the density of a soil-water mixing medium and b is a body force. Terzaghi defined the effective stress principle. Bishop and Blight [21] extended his definition and proposed the following equation for a saturated-unsaturated medium t
where o'ij is an effective stress, <5 - is Kronecker's delta, p is a unit weight of water and ^ is a pressure head. The parameter % is defined as follows 0
f
saturated zone unsaturated zone
HEAT FLOW
WATER
FLOW
DEFORMATION Storativity change
Figure 1
Coupled phenomena considered in T H A M E S
(3)
194
Rock Mechanics Continuum
Modeling
X is a nonlinear function of S (degree of saturation). It is usually determined by laboratory experiments. The validity of equation (2) is not definite and is still under debate. Here, we assume that x approximately equals S . Substituting equation (1) into equation (2), we obtain the equilib rium equation for an effective stress in a saturated-unsaturated geological medium r
r
Wj + Xtijpf^j
+ pbt = 0
(4)
where (x^tjPf^) is the term which means that a change of pressure head influences the equilibrium equation. Temperature effects can be implemented in a constitutive law for a solid medium. For an isotropic linear elastic material, the Duhamel-Neuman relationship can be used and the following constitu tive law is obtained °ij = C e - Pdij(T - T ) (5) ijkl
kl
0
where /? = (3A + 2/z)a. C is an elastic matrix, s is a strain tensor, T i s a temperature, k and \x are Lame's constants and a is the thermal expansion coefficient. Subscript 0 means that the parameter is in the reference state. The infinitesimal strain-deformation relationship is ijkl
kl
«« = i ( w j + Mi,*)
(6)
k
where u is a deformation vector. Substituting equations (5) and (6) into equation (4), we get the stress equilibrium equation, which also takes into account the effects of heat transfer and pore pressure change x
LiC (u j ijkl
[—P$ij(T 7.23
k
+ u,. )-/M<,(r-7o) + xSvPtift.j
+ P^ = 0
fc
(7)
— r ) ] , j is the term which means that heat transfer influences the equilibrium equation. 0
Continuity Equation for Ground Water
Since the ground is considered as a continuous body, the water movement is described by Darcy's law, which is valid for laminar flow in saturated porous ground. The behavior of water movement in an unsaturated zone is a phenomenon in a medium consisting of multiphases, i.e. liquid, gaseous and solid phases. However, water movement in both saturated and unsaturated media can be treated as a liquid water flow by assuming that the permeability is a function of the volumetric water content and that the storativity in an unsaturated medium changes as a relation of the volumetric water content and suction. This method has been commonly used for flow analyses in a saturatedunsaturated medium. An equation of continuity for ground water in a saturated-unsaturated zone is derived from the mass conservation equation of the fluid in a porous medium as follows —
=
~ (pf i).t
()
v
8
where 9 is the volumetric water content, t is time and v is a velocity vector. An equation of motion for ground water can be expressed by Darcy's law. That is t
v
t
= -kiO^hj
(9)
in which fc(0) - is a permeability tensor which is a function of 9 and h is a total head. The total head is the sum of pressure head ij/ and elevation head Z as follows o
h =
\jf + Z
(10)
Volumetric water content is a function of degree of saturation S and porosity n and is expressed as follows r
9 = nS
(11)
T
Substituting equations (9), (10) and (11) into equation (8), we obtain the following equation {Pfk(9)ij(^ + Z)j}
ti
=
7-pfnS
r
(12)
Thermal-Hydraulic-Mechanical
Coupling Analysis of Rock Mass
195
The right-hand side of equation (12) is expanded to (13) The first term of the right-hand side represents a density change of pore water. The second one represents a change of the skeleton of a porous medium. The third one represents a change of water storage capacity. Considering the compressibility and thermal expansivity of water, the density of water can be expressed as follows Pt = PfoLl - P A T - T ) + PP(P - P ) ] 0
(14)
0
where P is the pore water pressure, p is the reference density at P = P and T = T . /? and /? are the thermal expansivity and compressibility of water, respectively. They are defined as f 0
0
0
. . (P = constant)
h
. (T= constant)
=
T
P
(15)
(16)
Eaton [22] assumed that buoyancy could be neglected in an unsaturated zone. Adopting this assumption, we set /? = 0 in a unsaturated zone. Combination of the first term of equation (13) and equation (14) gives T
(17) The pressure head ij/ is related to the pore water pressure as follows (18) Taking into account equations (10) and (18), equation (17) can be modified t o give nS
r
p nS f0
r
(19)
where the first and second terms of the right-hand side of equation (19) indicate a density change of pore water due to temperature change and a total head change, respectively. Assuming that the strain is infinitesimal, the second and third terms of the right-hand side of equation (13) are expressed as follows
(20) Equation (12) is modified by using equations (19) and (20) (21) where C(^) is the specific water content and is defined as (22) Equation (21) is an equation of continuity for ground water which takes into account the compress ibility of the ground and density change by heat transfer.
196 7.2.4
Rock Mechanics Continuum
Modeling
Energy Conservation Law
In general, the ground consists of materials with three phases: solid, liquid and gas. It is n o t easy to analyze heat transfer behavior since the thermal process is different in each phase. W e assume that a pore in a porous medium contains only material in the liquid phase. This means that the ground water does not change its phase from liquid to gas or vice versa. With the above assumptions, we derive an energy conservation law without the effect of viscous dissipation for ground water. This process is based upon the one proposed by Bear a n d Corapcioglu [17]. Taking into account the existence of an unsaturated zone, an equation of energy conservation is written as follows (23)
where C is the specific heat and / is the heat flux by conduction. Subscript f means 'fluid'. In equation (23), the first term of the left-hand side shows the time dependency of the energy, and the second term shows the energy change due to heat convection. T h e first term on the right-hand side expresses the energy change by heat conduction and the second term shows the reversible energy change caused by compression. Similarly, the energy conservation law for a solid is written as v
-
V(l-*)/.-(1
-n)P7
(24)
where the subscript s means 'solid'. In equation (24), the second term on the right-hand side indicates a reversible energy change caused by deformation. Faust and Mercer [23] proposed that the movement of water through porous media is so slow and the surface areas of all phases are so large that it is reasonable to assume that local thermal equilibrium among phases is achieved instantaneously. If this assumption is permitted, the following equation holds T
=
T
s
=
(25)
T
f
Equations (23), (24) and (25) can then be combined to obtain the following energy conservation equation {nS p C T
=
{
y{
-V-{nS / r
f
{ n S p C t * + (1 - n ) p C t > } - V r
n)p C
+ (1 -
s
+ (1 -
vs
r
f
vf
Viv-(1
nS l
n)/ } -
r
B
s
v s
8
-n)Tp
(26)
When it is assumed that Fourier's law is valid for heat conduction, the following equations are given /
f
/
s
=
-
K VT Tf
= - K VT
(27)
Ts
where K is the coefficient of heat conduction. The term 3 P / 9 r i n equation (26) can be modified by using equations (15) and (16) as follows T
=
constant
(28)
Neglecting the velocity of the solid, equation (26) may be rewritten by means of equations (6), (9), (27) and (28) in the following form h nStPfCvfVuTi
-
(KjnijTijli
4
(29)
where ( p C ) v
m
and K
Tmij
are expressed as follows (pC ) v
Kjmij
m
=
nSrPfCvf + (1 -
=
nS K T
Tfij
+
(1 —
n)p Cys
(30)
n)K
(31)
s
Taij
Thermal-Hydraulic-Mechanical
Coupling Analysis of Rock
Mass
197
Equation (29) is an energy conservation law in which stress-deformation and ground water flow are considered. The first term on the left-hand side shows a dependency on energy and the second one shows an energy change due to heat convection. The first, second and third terms in the right-hand side express an energy change due to heat conduction, pore water pressure change and reversible energy change caused by solid deformation, respectively.
7.2.5
Governing Equations and Numerical Technique
Equations (7), (21) and (29) represent the governing equations for a coupled thermal, hydraulic and mechanical problem. We use these equations by means of the total head expression as follows + n,. ) - pd^T
liC (u ijkl
{pfk(0)tjhj}
ti
+
kJ
k
-
- T ) + zSijPth].j + PM = 0 0
(32b)
Pf nS p gP ^ 0
nSrPfCtfVfiTi
-
r
{
(32a)
P
(KjnijTijli
+
(32c) where p = (1 — n)(p — S p ) and p is the density of the solid phase. It is necessary to set the following initial and boundary conditions to solve equations (32). Initial conditions s
s
r
f
s
u (x,t) t
=
u (x,0)
h(x,t)
= h(x,0)
t
(33)
Boundary conditions displacement u (x,t)
(34a)
=u {x,t)
t
t
traction =
(Tij(x,t)nj(x)
(34b)
fi(x t) 9
total head h(x,t)
=
(34c)
h{x,t)
flow rate {*(0)iA;}"i =
temperature
T(x t) 9
=
(34d)
-C(*0
(34c)
f{x,t)
heat flow KjmijTjni
=
-
(34f)
Q (x,t) T
where x is a position vector, n is the unit normal vector, M, is a known^displacement, h is a known head, f is a known surface traction, Q is a prescribed flow rate and Q is a prescribed heat flow. The Galerkin-type finite element technique is employed to formulate a finite element discretiza tion. This method for deriving the element matrix equation is a general approach which can be used for problems for which the variational principles d o not exist. Moreover, in the case that an unsaturated medium is included in the problem, and iterative calculation is needed for the nonlinear permeability. As shown in the equations (32), the element matrix is not symmetric. This means that the coupled terms between each pair of three phenomena, i.e. water flow, heat transport and deformation, d o not have the same effect on the overall behavior. }
x
7.3
T
VERIFICATION O F THE C O D E
T H A M E S has to be verified for fully coupled hydraulic, mechanical and thermal problems. However, these problems can not be solved analytically at the present time. Therefore, the code has been verified with the available analytical and experimental results. The problems considered in the verification are as follows.
Rock Mechanics Continuum
198
Modeling
(i) A one-dimensional consolidation problem. An analytical solution to this problem can be found in Terzaghi [24]. (ii) The stress-strain problem of a thick-walled circular cylinder subjected to uniform internal and external pressures. An analytical solution to this problem can be found in Timoshenko and Goodier [25]. (iii) The thermal stress problem of a thick-walled circular cylinder subjected to a uniform temper ature gradient. An analytical solution to this problem can be found in Boley and Weiner [26]. (iv) A one-dimensional heat conduction problem. An analytical solution to this problem can be found in Smith [ 2 7 ] . (v) A coupled thermal and hydraulic problem. Experimental results can be found in Sato [28]. He experimented with hot water seepage in a saturated sand sample. (vi) A saturated-unsaturated seepage problem. Quasi one-dimensional, nonsteady seepage flow through a sand box was examined by Akai and U n o [ 2 9 ] . Excellent agreement was achieved between the finite element solutions by THAMES and the analytical and experimental solutions. The numerical results are presented in Figures 2 - 7 . From these calculations, we conclude that THAMES is verified for the fundamental functions contained in it. In the following section, we try to solve actual rock mechanics problems to show the capability of this numerical code.
7.4
N U M E R I C A L EXAMPLES
7.4.1 Macro-permeability Test in the Stripa Project At Stripa in Sweden, macroscopic permeability test was carried out in the Stripa project. A 33 m length of the ventilation drift was sealed off and equipped with a ventilation system whose temperature was controlled to evaporate all water seeping into the room. The water seepage was r
Degree of consolidation
O
0
0.2
0.4
0.6
0.8
Time factor Figure 2
Comparison with the analytical solution of the one-dimensional consolidation
10
R
•
Radial stress
problem
Numerical
Analytical
Radial distance Figure 3
Comparison with the analytical solution of the stress-strain problem of a thick-walled circular cylinder subjected to uniform internal and external pressures (E = 1000, v = 0.3)
199
Coupling Analysis of Rock Mass
Tangenfiol
stress
Thermal-Hydraulic-Mechanical
Figure 4
Comparison with the analytical solution of the thermal stress problem of a thick-walled circular cylinder subjected to a uniform temperature gradient: T = 5r (T is temperature, r is radial distance), a = 0.0001
Analytical Numerical
t = O
R =
1 (0
< x<
* = O , I R = O ( R > O )
1) I
tf
J C
T
V
=1 =1
t
•
0.001
V
0.002
•
0.005
•
0.010
A
0.050
O 0.100 A 0.200
x Figure 5
Comparison with the analytical solution of the one-dimensional heat conduction problem
Figure 6 Comparison with the experimental results of the coupled thermal and hydraulic problem (L is the length of the specimen, x is the distance from the top of the specimen, T is the temperature of the injected water and T is the temperature of the water in the specimen at the initial state) 0
x
200
Modeling
(cm)
Rock Mechanics Continuum
(cm) Figure 7
Comparison with the experimental results of the saturated-unsaturated seepage problem
determined from measurements of the mass flow rate and the difference in the humidities of the drift. The pressure gradients in the rock walls were measured in holes that radiated out from the sealed room. For a coupled finite element analysis, the rock cavern is simply modeled as shown in Figure 8. Excavation of the underground space is simulated with the following stress-dependent permeability proposed by Iwai [30] (35) where fe is the intrinsic permeability under the stress
c
od
co
Table 1 Parameter E r ^vs A
Value
Parameter
Value 4.9x10* MPa 8371kg- X 0.0276 1
Parameters at the Initial State
-
V 1
0.33 3xl0" ' k J m - ^ - ^ C 0.00217 {
^Ts
c
Parameter
Value
k a t'
lxKT m 8xl0" 0.728
0
- 1
8
30°C test, k
0
nonlinear and all nonlinear
20°C test, k
0
nonlinear and oil nonlinear
- 2 0 ° C test, / /
6
2
6
nonlinear
36°C test, ju nonlinear 20°C test, observed 30°C test, observed
The pressure head at the drift
is set to be 0
,18 m |
Figure 8
Schematic model of macro-permeability test and comparison of the water tables by numerical results and results observed in the field
Thermal-Hydraulic-Mechanical
Coupling Analysis of Rock
Mass
201
Computed distribution of the horizontal permeability at the middle height of the drift is shown in Figure 9, in comparison with in situ measurements. Reduction of permeability in the figure seems to indicate closure of the fracture due to an increase in the circumferential stress. When the room temperature goes up to 30 °C, the rock matrix block expands and the fractures close. As a result, the permeability of the rock mass reduces as shown in Figure 9. It also causes the higher water level as shown in Figure 8. The calculated flow rate seeping into the room at different temperatures is shown in Table 2. Since the hydraulic boundary condition is not well known, a quantitative examination is difficult. However, the general trend of the change in inflow rate is similar to the measured trend. In a case where only the thermal and mechanical dependencies of k are considered, the reduction rate of the water inflow at a room temperature of 20 to 30°C is 39.4%, while the reduction rate is 32.9% in the case where all the dependencies are considered. This is due to the effect of the increase in p. It is found from these results, therefore, that a change in the inflow rate is considerably influenced by the temperature dependency of the kinematic viscosity, p. Figure 10 shows the calculated vertical permeability distribution in a horizontal section from the wall of the vertical hole. It is seen that the high permeability zone, which corresponds to a loosened area near the wall, is developed out to a radius of 7 m from the wall surface. 0
7.4.2
Buffer Mass Test in the Stripa Project
1
Horizontal permeability (m s )
Canisters of high-level waste will be stored in shallow vertical boreholes in the floor of a horizon tal tunnel. The annulus between the canisters and the borehole will consist of a highly compacted clay (HCC). The buffer material, bentonite clay, has a very low permeability and high swelling capacity when absorbing water. The objective of the Buffer mass test is to verify the barrier function of highly compacted bentonite and sand under real conditions. Therefore, it is necessary to test the behavior of the integrated system of heat-producing canisters {i.e. electric heaters), buffer material, rock and ground water in a representative crystalline rock mass.
Kelsall et al. 20°C rest, k nonlinear and all nonlinear 20°C test, /u nonlinear 30°C Test, k nonlinear 30°C rest, all nonlinear 30°C Test, /x nonlinear 20°C Test, observed 30°C tesT, observed 0
0
5
10
15
R a d i a l distance ( m ) Figure 9
Comparison of horizontal hydraulic conductivities calculated by T H A M E S with those inferred from the macropermeability test.
Table 2
Variation of Fluid Flow into the Drift with Increase in Temperature
Nonlinear
property
k only (calculated) p. only (calculated) k and p (calculated) Observed 0
0
Fluid inflow (mLmin ) 20 °C 30 °C 1
31.0 58.5 31.0 50.0
18.8 63.3 21.0 42.0
202
Modeling
o
1
Vertical permeability (m s" )
Rock Mechanics Continuum
o • • •
test, test, test, test, test,
k nonlinear nonlineor k nonlinear all nonlinear fi nonlinear 0
0
q
o,
A
20°C 20°C 30°C 30°C 30°C
10
, c
5
10
15
Radial distance (m) Figure 10
Vertical hydraulic conductivities calculated by T H A M E S
Figure 11
Finite element mesh of Buffer mass test analysis
Numerical analyses have been performed with the axisymmetric finite element mesh shown in Figure 11. The electric heater intensity is assumed to be 600 W. A bentonite powder is placed in the slot between the HCC and the host rock. Table 3 shows the parameters of the HCC. Figure 12 shows the distributions of porosity n, degree of saturation S , heat conductivity K and permeability k in the borehole. The swelling HCC produces an increase in porosity. S goes up by 100%, after a reduction in its value, depending on the combination of the conditions of n and swelling. The heat conductivity gradually decreases with the increase in the porosity. r
T
T
Table 3 Property Density Porosity Young's modulus (E) Poisson's ratio (v)
Value
Property
2150 K g m " 0.42 9.8xl0 MPa 0.3 3
3
Data of H C C Used for Analysis
Thermal expansion coefficient Specific heat Thermal conductivity Permeability
Value
fcOxio-^cr
1
1220 J kg" ^ C "
1
11 0. 4" 6 xml < r k J n r s - C r 3
1 9
2
1
l0
1
Thermal-Hydraulic-Mechanical
Coupling Analysis of Rock 1 . 0 r-
(a) saturation
V
\
Porosity
\
months
1 day
/
o
Degree
p
of
\
203
(bj_ 4
4 months
Mass
weeks
/
1 day
n — 1
10
20
0
Radial distance (cm)
Figure 12
*
1
(m s~ )
\
(d)
\ \ ^4
months
o
\
\
\ \
-
O
Hydraulic conductivity
1
Thermal conductivity ( W m
10
20
\
1
°C" )
en
io'V
Radial distance
10
Radial distance (cm)
20
(cm)
10
Radial
20
distance (cm)
Calculated property profiles as a function of radial distance from the heater
The comparison between the numerical results and field test results [31] is shown in Figures 13(a-c) for temperature, water content and swelling pressure, respectively. The temperature distribu tion obtained in the field measurement is well predicted by the computation, as shown in Figure 13(a). The distribution of water content disagrees with the measurement. This is possibly because the evaporation of water is neglected in the analysis. Therefore, it may be concluded that the water content may be strongly influenced by water movement due to evaporation. Figure 13(c) indicates that the initial swelling pressure calculated is larger than the measured one. This may be caused by overestimation of suction at the initial state.
7.4J
An Hypothetical Deep Nuclear Waste Repository
High-level radioactive wastes are planned to be disposed of in deep geological formations after an interim storage. The behavior of ground around such a repository is analyzed here. In the model calculations, the region is considered as having an area of 2000 by 2000 m and a height of 1500 m, as shown in Figure 14. The rock mass is assumed to be granite. The repository is assumed to be located within a 500 m by 500 m square at a depth of 1000 m. Table 4 shows the parameters used in the analysis. This model is analyzed using the two-dimensional and the three-dimensional codes. Figure 15 is a schematic diagram of the two-dimensional model and Figure 14 is that of the three-dimensional one. The repository is simulated by the heat sources indicated by black circles in Figures 14 and 15. The heat source is kept at 100 °C during the analyses. The temperature in the initial state is a function of the depth (GL — 1000 m), and the temperature at the repository is set at 50 °C. The
204
Rock Mechanics Continuum
Oo)
Water contents distribution (%)
(a)
Modeling
b)
uoj4.nqu4.sip a_n4.D_9duj8i
Radiol distance (cm)
Radial distance (cm)
(DdlAI) 4.uawdo|8Asp ajnssajd 6UJII9MS
(c)
(d) Pressure head distribution (m)
• Numerical — BMT hole no. 1 — BMT hole no. 5 — BMT hole no. 4 S 6
/ -
0
10
Time after start Figure 13
• Numerical, o Numerical,
2.1 years 0.9 years
10
20
Radial distance
(months)
(cm)
Comparison with observed values in Buffer mass test
Fixed temperature and total head condition
z (m)
No heat and water flux condition
0
500 700 1000
2000
x (m) Figure 14
Three-dimensional finite element mesh of hypothetical repository
permeability is also assumed to be a function of depth, as defined by equation (35). Total water head is constant for all regions in the initial state. The boundary conditions are indicated in Figures 14 and 15. The distribution of displacement vectors at y = 0 (frontal profile) due to thermal stress after an elapsed time of 100 years is shown in Figures 16 and 17. The three-dimensional analysis gives similar
Thermal-Hydraulic-Mechanical Table 4 Parameter Young's modulus Poisson's ratio Intrinsic permeability x direction y direction z direction
Parameter
4.9 x 1 0 M P a 0.3 KT m 10" m KT m 1 2
1 2
Value 2700kgm" 0.02 417Jkg- C" 6.0xl0- kJm2.1xlO- kJm6.0 x l O " ^
Unit weight Initial void ratio Specific heat of soil Heat conductivity of fluid Heat conductivity of soil Thermal expansion coefficient
3
2
2
2
No heat and water flux condition
3
l 0
Fixed temperature and total head condition
g 1000 700
N
HI No water flux and fixed temperature condition |
500 H
r i
300 0 -I i 0
500
1 — i — i — 700 1000
x (m) Two-dimensional finite element mesh of hypothetical repository
/
Figure 16
////
Distribution of displacement vectors by three-dimensional analysis
1 / //// /
, i i
(////
10~ m 2
/ / A S
Figure 17
1
4
l o
3
1 0
6 0
_____
1500 -
Figure 15
205
Parameters Used in the Analysis of Hypothetical Repository
Value
1 2
Coupling Analysis of Rock Mass
Distribution of displacement vectors by two-dimensional analysis
1
C sC- s_ 1
1
1
1
206
Rock Mechanics Continuum
Velocity, m s
Figure 18
*
* * *
*
% * »
t
• * >
*
•
•
Modeling
1
>
Distribution of ground water velocity vectors
results as the two-dimensional one at the vicinity of the heat source. However, with increasing distance from the heat source the three-dimensional analysis yields smaller displacements. The reason for the divergence of the results is due to the assumption of plane strain condition in the two-dimensional case. The thermal stresses are isotropic and the repository has a finite length for the y direction in Figure 14. The ground water velocity is also influenced by the heat generated at the repository. The distribution of velocity vectors at an elapsed time of 10000 years is given in Figure 18. For a location close to the heat source, the ground water flows upward and, for a remote location, the effect of the heat source on the ground water velocity is small.
7.5
SUMMARY A N D CONCLUSIONS
In this chapter, the coupled thermal, hydraulic and mechanical problem has been discussed. There are many remaining subjects to be studied in coupled problems. For example, the validity of the principle of the effective stress for the fractured rocks is one such problem. If the deformation of the rock mass is strongly dependent on the fracture structure, the pore pressure in the fracture has a great effect on the deformation of the rock mass. However, it seems very difficult to observe the pressure in a fracture under the deformation process in an in situ test. Since the ground water does not homogeneously fill the volume of a fractured rock, the effect of the ground water on the deformation behavior may be anisotropic. Ground water flow in an unsaturated fractured medium is a difficult phenomenon to observe. There have not been many studies for the measurement of the unsaturated properties of a rock matrix. The unsaturated flow problem in a fracture also remains to be solved. The effect of the stress field on the permeability of the rock mass is not well understood. The exponential function of a void ratio has been used for the permeability of clay and sand strata in many cases [32]. For fractures, many studies have been carried out [30, 33-35]. However, a new concept of fluid flow in a fracture plane, such as a channel flow, has been introduced as a result of a recent field observation [36]. When fracture flow is considered as a phenomenon restricted to a plane, the effect of the deformation on the permeability has to be considered as a change in the area of the channels in the fracture plane. The effect of temperature change on the deformation and permeability of a rock mass has not been well investigated. Plastic deformation of the rock matrix due to temperature change was observed by Shimooka et al. [37]. However, a quantitative examination of the phenomenon has not been carried out, and it is difficult for the geomechanical problem to form a constitutive model considering the plastic deformation induced by the temperature change, because not much data have been obtained from the fundamental experiments.
Thermal-Hydraulic-Mechanical
Coupling Analysis of Rock
Mass
207
Given this state of the art, numerical analyses are helpful in understanding the phenomena by comparing the model results with observations. By examining the difference between observed and calculated results and the assumptions used in the model, the phenomena can be studied from various points of view. These complex phenomena will then be better understood by refined experiments based on the results of the analyses.
7.6
NOMENCLATURE
tit ft T i Q Q (pCv)m a P b PP P X C(i//) Ci C dij e h J k(9)ij K f^Tmij X and p n n P 9 p p alj Gij S T t u Vi ij/ t
T
t
T
ijk v
kl
Tij
t
f
T
t
known displacement components known head known surface traction components position vector prescribed flow rate prescribed heat flow defined in equation (30) thermal expansivity coefficient = (3A + 2p)a body force vector components compressibility of water thermal expansivity of water defined in equation (3) specific water content elastic matrix components specific heat Kronecker's delta strain tensor components total head heat flux by conduction permeability tensor components which are a function of 9 coefficient of heat conduction defined in equation (30) Lame's constants porosity unit normal vector components pore water pressure volumetric water content density of a soil-water mixing medium unit weight of water effective stress components total stress components degree of saturation temperature time deformation vector components velocity vector components pressure head
Subscript f means 'fluid', s means 'solid' and 0 means that the parameter is in reference state. 7.7 1. 2. 3. 4. 5.
REFERENCES Christian J. T. and Boehmer J. W. Plane strain consolidation by finite elements. J. Soil Mech. Found. Div., Am. Soc. Civ. Eng. 96, 1435-1457 (1970). Sandhu R. S. and Wilson E. L. Finite element analysis of seepage in elastic media. J. Eng. Mech. Div., Am. Soc. Civ. Eng. 95, 6 4 1 - 6 5 2 (1969). Tsang C. F. Coupled behavior of rock joints. In Rock Joints (Edited by N. Barton and O. Stephansson), pp. 5 0 5 - 5 1 8 (1990). Noorishad T. N., Witherspoon P. A. and Brekke T. L. A method for coupled stress and flow analysis of fractured rock masses. Geotechnical Engineering Publication, N o . 71-6, University of California, Berkeley (1971). Ohnishi Y. and Ohtsu H. Coupled stress flow analysis of discontinuous media by finite elements (in Japanese). Trans. JSCE 322, 111-120 (1982).
208 6. 7. 8. 9. 10. 11. 12. 13. 14. 15. 16. 17. 18. 19. 20.
21. 22. 23. 24. 25. 26. 27. 28. 29. 30. 31. 32. 33. 34. 35. 36. 37.
Rock Mechanics Continuum
Modeling
Barton N., Chryssanthakis P. and Monsen K. Validation of M U D E C against Colorado school of mines block test data. SKB Technical Report, 88-14 (1988). Fairhurst C. and Hart R. D . Verification and validation of coupled mechanical/water flow effects in rock m a s s e s - s o m e possibilities and limitations. Proc. GEOVAL-87 Symp., Stockholm, pp. 527-545 (1987). Oda M. An equivalent continuous model for coupled stress and fluid flow analysis in jointed rock masses. Water Resour. Res. 22, 1945-1956 (1986). Shimooka K., Ishizaki K., Okamoto M., Kumata M., Araki K. and Amano H. Thermal characteristics of rocks for high-level waste repository (in Japanese). JAERI-M 9247 (1980). Shen B. and Stephansson O. Rock mass response to glaciation and thermal loading from nuclear waste. Proc. GEOVAL-90 Symp., Stockholm, pp. 550-558 (1990). Utsugida Y. Coupled analysis of flow and heat around a high-level nuclear waste repository. In 5th Int. Conf. Numerical Methods in Geomechanics, Nagoya, pp. 711-716. Balkema, Rotterdam (1985). Hamajima R., Hunyu K., Kusabuka M., Watanabe M. and Koide H. Coupled analysis of stress and heat transfer for three-dimensional rock masses (in Japanese). Proc. 19th Symp. Rock Mech., Tokyo, pp. 221^225 (1987). Mercer J. W. and Faust C. R. Simulation of water- and vapor-dominated hypothermal reservoirs. In 50th Annual Fall Meeting of the Soc. Pet. Eng. AIME, Dallas, Paper SPE 5520 (1975). Witherspoon P. A., Neuman S. P., Sorey M. L. and Lipman M. J. Modeling geothermal systems. Technical Report N o . LBL-3263 (1975). Hodgkinson D . P. A mathematical model for hydrothermal convection around a radioactive waste depository in hard rock. Ann. Nucl. Energy 1, 3 1 3 - 3 3 4 (1980). Wang J. S. Y., Tsang C. F., Cook N . G. and Witherspoon P. A. A study of regional temperature and thermohydrologic effects of an underground repository for nuclear wastes in hard rock. J. Geophys. Res. 86, 3759-3770 (1981). Bear J. and Corapcioglu M. Y. A mathematical model for consolidation in a thermoelastic aquifer due to hot water injection or pumping. Water Resour. Res. 17, 7 2 3 - 7 3 6 (1981). Hart R. D . A fully coupled thermal-mechanical-fluid flow model for nonlinear geologic system, Ph.D. Thesis, University of Minnesota (1981). Noorishad J., Tsang C. F. and Witherspoon P. A. Coupled thermal-hydraulic-mechanical phenomena in saturated fractured porous rocks-numerical approach. JGR, J. Geophys. Res. 89, 1 0 3 6 5 - 1 0 3 7 3 (1984). Ohnishi Y., Shibata H. and Kobayashi A. Development of finite element code for the analysis of coupled thermohydro-mechanical behaviors of saturated-unsaturated medium. In Proc. Int. Symp. Coupled Processes Affecting the Performance of a Nuclear Waste Repository, Berkeley, pp. 263-268 (1985). Bishop A. W. and Blight G. E. Some aspects of effective stress in saturated and partly saturated soil. Geotechnique 13, 177-193 (1963). Eaton R. R. A numerical method for computing flow through partially saturated porous media. Numerical Methods in Thermal Problem, Vol. 3, pp. 9 1 1 - 9 2 0 (1983). Faust C. R. and Mercer, J. W. Geothermal reservoir simulation 1. Mathematical models for liquid- and vapordominated hydrothermal system. Water Resour. Res. 15, 2 3 - 3 0 (1979). Terzaghi K. Settlement and consolidation of clay. Eng. News Rec. Vol. 26 (1925). Timoshenko S. and Goodier J. N . Theory of Elasticity, p. 59. McGraw-Hill, N e w York (1951). Boley B. A. and Weiner J. H. Theory of Thermal Stresses, pp. 288-291. Wiley, London (1960). Smith G. D . Numerical Solution of Partial Differential Equations. Clarendon Press, Oxford (1965). Sato K. Experimental determination of transfer parameters of heat flow through porous media by means of a new designed apparatus in laboratory (in Japanese). Trans. JACE 320, 5 7 - 6 5 (1982). Akai K. and U n o T. Study on the quasi-one-dimensional, non-steady seepage flow through soil (in Japanese). Trans. JSCE 125, 14-22 (1966). Iwai K. Fundamental studies of fluid flow through a single fracture, Ph.D. Dissertation, University of California, Berkeley (1976). Push R. and Borgesson L. Final report of the Buffer mass t e s t - V o l u m e 11; test results. SKB Technical Report 85-12 (1985). Lambe T. W. and Whitman R. V. Soil Mechanics. Wiley, London (1969). Snow, D. T. Rock fracture spacings, openings and porosities. J. Soil Mech. Found Div., Am. Soc. Civ. Eng. 94, 7 3 - 9 1 (1968). Ohnishi Y. Laboratory measurement of induced water pressures in jointed rock, Ph.D. Thesis, University of California, Berkeley (1973). Louis C. Introduction a Thydraulique des roches, Ph.D. Thesis, Paris (1976). Abelin H., Neretnieks I., Tunbrant S. and Moreno L. Final report of the migration in a single fracture-experimental results and evaluation. Stripa Project 85-03, SKB (1985). Shimooka H. Numerical methods for analysis of temperature rises and thermal stresses around high level radioactive waste repository in granite (in Japanese). J. At. Energy Soc. Jpn. 24, 889-895 (1982).