Simulation of a fully coupled thermo–hydro–mechanical system in freezing and thawing rock

International Journal of Rock Mechanics and Mining Sciences 36 (1999) 563±580

www.elsevier.com/locate/ijrmms

Simulation of a fully coupled thermo±hydro±mechanical system in freezing and thawing rock K.M. Neaupane*, T. Yamabe, R. Yoshinaka Department of Civil and Environmental Engineering, Saitama University, 255 Shimo-Ohkubo, Urawa Saitama, Japan Accepted 3 March 1999

Abstract This paper describes the application of a numerical model of thermal±mechanical±¯uid ¯ow coupling system to the simulation of laboratory freezing and thawing experiments on rocks. A theoretical formulation that accommodates a linear stress±strain constitutive relationship is presented and a two-dimensional (plane stress) numerical modeling is performed based on the ®nite element method applied to thermo±poro±elasticity. As the primary objective of this research is to simulate the freezing and thawing process, the developed code takes into account the phase change of pore water during freezing. It is found from the numerical simulation that a relatively good prediction can be made of temperature transfer and deformation behavior within the elastic deformation limit. In some cases, however, large deformation results from the freezing processes, leading to material failure. Future research should therefore take into account irreversible and plastic e€ects. # 1999 Elsevier Science Ltd. All rights reserved.

1. Introduction 1.1. General Increased consumption of lique®ed natural gas (LNG) in Japan, Europe and America has focused attention on the problem of its storage. Semi-ground storage tanks have generally been used for the storage of LNG at a temperature of ÿ1628C. This has involved use of large areas of reclaimed land in Japan, where availability of land is scarce. An alternative solution to this problem would be to store LNG underground [1]. This, in turn, requires a comprehensive understanding of frozen ground phenomena such as induced deformation, temperature distribution, interaction of frozen ground and structure (in-ground storage tank), etc. for the design and safety of the storage system. When the ground freezes, thermal stresses and displacement are induced, along with volume change * Corresponding author. Tel.: +81-48-8583544; fax: +81-488587374.

due to expansion of water upon freezing and pore water moves into unfrozen regions towards the freezing fronts. The response of rock masses to this engineering problem is a coupled phenomenon involving thermal (T), hydrological (H) and mechanical (M) processes. The natural processes of freezing and thawing in snowy mountainous areas in Japan and elsewhere are responsible for slope failures, slope deterioration and weathering. Floor heaving in the plain areas of lower temperature zones is also attributed to these freezing and thawing processes. As the ground undergoes cycles of freezing and thawing, the pore water volume is altered and deformation is induced. The ground response to these phenomena may also require a THM feedback coupling system to describe them fully. Moving-boundary problems that are also known as Stefan problems are associated mainly but not exclusively with ¯uid ¯ow through porous media [2]. Thermo±hydro±mechanical coupling in conditions below freezing clearly demands to be treated as a moving-boundary two-phase case in which heat ¯ow occurs

0148-9062/99/$ - see front matter # 1999 Elsevier Science Ltd. All rights reserved. PII: S 0 1 4 8 - 9 0 6 2 ( 9 9 ) 0 0 0 2 6 - 1

564

K.M. Neaupane et al. / International Journal of Rock Mechanics and Mining Sciences 36 (1999) 563±580

Nomenclature Cijkl C (v) C (w) K K (T) S (r) T Ui V f g h h…x† i kij n P qi s t u v (r) i v (f) i v (s) i

elastic constant tensor speci®c heat at constant volume speci®c water content kinetic energy coecient of heat conduction degree of saturation temperature displacement vector volume Helmholtz free energy acceleration due to gravity total head radiant heat ¯ux vector permeability tensor porosity hydrostatic pressure conductive heat ¯ux vector entropy, density time speci®c internal energy relative velocity vector ¯uid velocity vector solid matrix velocity vector

Greek symbols a linear thermal coecient thermal coecient tensor bij compressibility b (P) expansivity b (T) Kronecker's Delta dij w Bishop's parameter

both in water and ice. As the primary objective of this paper is to solve the problem involving freezing and thawing of porous media, the developed numerical code takes into account the phase change of ¯uid. The modeling of thermo±hydro±mechanical (THM) processes is at an early stage of development; considerable work is yet to be done for simulating real ®eld problems. We still lack applicable test cases for the veri®cation and validation or evaluation of the developed code to apply them in a realistic geological system. This paper aims at providing freeze±thaw experiments as a means to test the validity of the fully coupled THM numerical model based code. A two-dimensional computer code HEAT has been developed for plane stress poro±elastic assumption in which deformation, ¯uid ¯ow and heat ¯ow are fully coupled. The numerical scheme employed in the code accommodates the phase change of water from liquid to solid and vice versa. Veri®cation of the program

eij j g l, m y r sij s(C) ij s(D) ij x

strain tensor pressure head unit weight Lame's elastic constants volumetric water content mass density stress tensor conservative stress tensor dissipative stress tensor elevation head

Superscripts f (properties associated with) water, ¯uid s (properties associated with) solid constituent Subscripts i,j,k,l

@i (
)=(
),i @i (
)i=(
)i,i Relations h=x+j f= g…Pf †

indicial notation; summation convention has been used to represent vectors and tensors by indexed base written in indicial notation gradient of (
) divergence of (
)i total head in terms of pressure head and elevation head pressure head expression in terms of hydrostatic pressure P and density of ¯uid g… f †

has been carried out for a simple one-dimensional consolidation case, the case of a thermo±elastic problem associated with a thick-walled circular cylinder and a heat transfer problem. The primary aim of this paper is to compare the prediction obtained from the numerical code and that obtained from the freezing and thawing experiment. 1.2. Review of feed-back coupling system Geo-engineering problems related to radioactive waste storage, geothermal energy extraction, reservoir engineering, heat extraction and heat storage, etc. involve thermal (T), hydrological (H) and mechanical (M) processes [3,4]. The one-way coupling mechanism is a simple case that re¯ects the continuing e€ect of one or more processes on the others, whereas the twoway coupling process reveals continuing reciprocal (feed back) interaction among di€erent processes in a

K.M. Neaupane et al. / International Journal of Rock Mechanics and Mining Sciences 36 (1999) 563±580

complex way. The two-way interaction mechanism as explained by Hart and John is illustrated in Fig. 1 [5]. Thermo±mechanical coupling results from thermal stresses imposed on the rock by heat, whereas thermal¯uid ¯ow coupling results from expulsion of pore ¯uid from the rock matrix due to the thermal expansion of water [6]. It is to be noted that thermal currents generated in saturated media may be entirely attributed to the density di€erences of ¯uid but in unsaturated zone it is more complex and relative permeability and capillary e€ects are to be considered [7]. 2. Research background The e€ects of heat on mechanical and ¯uid ¯ow processes have long been the subject of research. Terzaghi's general theory of hydrodynamic consolidation itself implies the basic in¯uence of temperature on the rate of pore pressure dissipation and resultant strain. Later, Paaswell carried out thermal consolidation experiments and studied the e€ect of temperature on strain [8]. Since then, there has been substantial research done on this ®eld. In general, most of the research related to thermal e€ect on ¯uid ¯ow and mechanical processes falls into two categories: (1) consolidation of soil assuming ®xed soil water temperature and (2) heat ¯ow through a static soil water system. Thus the processes occurring within the system are assumed to be uncoupled. Geological materials typically exhibit non-linear behavior. Therefore, for many geo-engineering pro-

565

blems, it may not be sucient to consider that the thermal, mechanical and ¯uid ¯ow processes occurring within a geological system are essentially uncoupled. Rock mass response to these engineering problems may not be predicted by considering each process independently and hence a coupling of two or more processes may be required depending on the engineering problems involved. However, an analytical solution for a fully coupled problem is not possible without making several simplifying assumptions. Booker and Savvidou developed a semi-coupled solution for consolidation around a point heat source buried in a poroelastic medium using the Laplace transform techniques, though this solution neglects the contribution from mechanical terms to the heat equation, which may be insigni®cant for most practical cases [9]. Seneviratne et al. presented a numerical solution based on the ®nite element method to describe the fully coupled consolidation and heat ¯ow around a rigid cylindrical heat source buried in clay with the aim of understanding the major mechanism of soil behavior close to a buried canister of hot radioactive waste [10]. Giraud and Rouset also concentrated their work on the numerical solution of a one-dimensional non-isothermal consolidation problem for fully saturated, porous, thermoelastic media. They concluded that coupling e€ects are most pronounced for low permeability and high porosity media and applied this concept for the prediction of temperature transfer and displacement ®elds in deep compressible clays with permeability as low as 10ÿ9 msÿ1 [11].

Fig. 1. Interaction mechanisms in a fully coupled thermal±mechanical±¯uid ¯ow system [12].

566

K.M. Neaupane et al. / International Journal of Rock Mechanics and Mining Sciences 36 (1999) 563±580

Hart and John presented a computational model that can simulate fully coupled thermal mechanical ¯uid ¯ow behavior for non-linear geological systems and analyzed the model by an explicit ®nite di€erence method [12]. Ohnishi and Kobayashi developed a ®nite element numerical method for coupled thermo±hydro± mechanical problems and implemented it in a twodimensional computer code on the assumption that the media is poro±elastic and in plain strain condition, that water in the ground does not change its phase [13]. As seen above, most research related to coupling problems has, so far, been in the higher temperature range. The little literature available on ground response behavior in sub-zero temperatures is concerned with the uncoupled analysis of stress, ¯uid and heat and is mostly associated with LNG storage problems.

3. Theoretical formulation In this paper, theoretical formulation of the fully coupled numerical model is based on the following assumptions . the medium is an anisotropic poro±elastic continuum and is fully saturated, . Darcy's law is valid for ¯uid ¯ow, . the rate of heat transfer between the solid and liquid phases is very slow and hence is neglected, . Fourier's law holds for heat ¯ux, . water density is a function of temperature and hydrostatic pressure, . solid/water mixture properties are described by the theory of mixture. The theory of mixture views material as a superposition of two or more constituents and at any time, each place in the mixture is simultaneously occupied by a particle from each constituent. For any property (
), volumetric fraction has been assumed to describe its mean value in the mixture as described by …
† ˆ n…
†…f† ‡ …1 ÿ n†…
†…s† ,

…1†

where `n' is porosity and superscript `f' and `s' are used to describe ¯uid and solid constituents respectively. For the mathematical formulation, Biot's theory of consolidation has been utilized [14]. It is assumed that the rock formation in consideration can be treated by the continuum approach. Fundamental laws of continuum mechanics are applied to describe responses associated with heat transfer, ¯uid ¯ow and mechanical stress problems in porous media. Details of these fundamental laws can be found in any text of solid mechanics or continuum mechanics [15]. An indicial

notation system has been used largely to describe the mathematical formulation. For details on the notation system refer to Chen and Saleeb [16]. 3.1. Continuity ¯ow equation The true velocity of a ¯uid experiences no matrix interference to ¯ow, but ¯ow through a porous medium is always resisted by the porous matrix in its way. Therefore the measurable velocity of the ¯uid is the relative velocity. Relative velocity is expressed in terms of actual ¯uid velocity and porous matrix velocity as …f† …s† v…r† i ˆ vi ÿ vi ,

…2†

(f) (s) where v(r) i is the relative velocity vector and vi and vi are actual ¯uid and porous matrix velocity vectors respectively. Mass is the property associated with every material continuum. In the continuity equation, mass is de®ned by a continuous function known as mass density, r(x, t ). The law of conservation of mass requires that the mass of a speci®c portion of the continuum remains constant and hence the material derivative of it is zero. The mass balance equation for a solid constituent takes the form

@ …1 ÿ n†r…s† ‡ @ i …1 ÿ n†r…s† v…s† i ˆ 0, @t

…3†

where r (s) is the mass density of the solid constituent and n is the porosity. As the derivative of mass density (r (s)) is zero, the solid mass balance equation after some manipulation becomes @n ˆ …1 ÿ n†@ i v…s† i : @t

…4†

Similarly, the ¯uid mass balance equation takes the form @ …f† nr ‡ @ i nr…f† v…f† i ˆ 0, @t

…5†

where r (f) is the ¯uid mass density. The equation of continuity for groundwater in the saturated±unsaturated zone is derived from Eq. (5) as @ ‰nS …r† r…f† Š ‡ @ i ‰nS …r† r…f† v…f† i Š ˆ 0, @t

…6†

where S (r) is the degree of saturation. Substituting Eqs. (2)±(5) into the second term of the left hand side in Eq. (6), the continuity equation can be rewritten as @ …s† …f† …r† …f† …r† @ n ‰nS …r† r…f† Š ‡ r…f† y@ i v…r† ˆ 0, i ‡ r S @ i vi ÿ r S @t @t

…7†

where volumetric water content (y ) is a function of degree of saturation and porosity as described by

K.M. Neaupane et al. / International Journal of Rock Mechanics and Mining Sciences 36 (1999) 563±580

y ˆ nS …r† :

…8†

The complete ¯uid mass balance equation after some manipulation, becomes nS …r†

@ r…f† @ S …r† …s† …f† …r† ‡ r…f† n ‡ r…f† y@ i v…r† i ‡ r S @ i vi ˆ 0: …9† @t @t

Of the four terms in Eq. (9), the ®rst derivative gives the rate of change of pore water mass density, the second gives the rate of change of degree of saturation, the third gives the transformation rate of solid constituent and the last term in the expression is for kinematics ¯ow. Before the ®nal expression is given for continuity ¯ow, the following points are to be considered: . the equation for the density of water at given temperature (T) and hydrostatic pressure (P) is …T† …P† r…f† ˆ r…f† 0 ‰1 ÿ b …T ÿ T0 † ‡ b …P ÿ P0 †Š

…10†

where b (T) and b (P) are expansivity and compressibility of water respectively and P0 and T0 describe the reference state of pressure and temperature respectively, . the relationship between changes in degree of saturation with respect to head change is r…f† n

@S …r† @S …r† @ j @j ˆ r…f† n  r…f† C …w† …j† @t @ j @t @t

…11†

where j is the pressure head and C (w) is de®ned as the speci®c water content, . the relative transformation of solid matrix during ¯ow is de®ned in terms of displacement components, Ui, as v…s† i ˆ

@ Ui , @t

…12†

. and the behavior of ¯uid movement through porous media is assumed to be de®ned by generalized Darcy's law ˆ ÿk…y†ij h,j , v…r† i

…13†

where k is the relative permeability and h is the total head. From Eqs. (10)±(13), expressing pressure head in terms of the pore pressure i.e. j=(P/r (f)g ), the ®nal form of the equation for continuity ¯ow is: ÿfnS

…r†

…T† @ T gr…f† 0 b

‡ r…f† C…j†

@t

‡ fnS

…r†

…P† …f† @ j gr…f† 0 b …r g†

@t

@j @Ui,i ÿ r…f† k…y†ij h,ji ˆ 0: ‡ r…f† S …r† @t @t

…14†

567

3.2. Equilibrium equation It is assumed that the medium is a continuous anisotropic elastic body and a volume of the elastic porous media is ®lled with water. The total applied stress distribution within this volume satis®es the equilibrium equation, @ j sij ‡ rbi ˆ 0,

…15†

where the component from acceleration has been omitted as it vanishes for quasi-static equilibrium. For ¯uid di€using through porous media, total stress may be decomposed into e€ective and hydrostatic components i.e. sij=sij 'ÿPdij where the hydrostatic component `P' is given by skk/3. Expressing hydrostatic pressure in terms of pressure head, the constitutive equation of equilibrium becomes …s 0 ij ÿ wg…f† jdij †,j ‡ rbi ˆ 0,

…16†

where w is Bishop's parameter which is a function of the degree of saturation, j and g (f) are pressure head and unit weight of ¯uid respectively and (wg (f)jdij),j is the term which means that the change of pressure head in¯uences the equilibrium equations. Deformations are also accompanied by a change in the temperature which can occur either from the deformation process itself or from external causes. For isotropic linear elastic materials, the Duhamel±Neumann relationship can be used and the following constitutive equation may be obtained [17] s 0 ij ˆ Cijkl ekl ÿ bij …T ÿ T0 †,

…17†

where Cijkl is the elastic constant tensor and bij is the thermal coecient de®ned by bij=(3l+2m)adij, in which l and m are two independent constants known as LameÁ constants and a is the coecient of linear thermal expansion. The last term of the equation means that heat di€erence in¯uences the equilibrium equation. In linear elastic theory, linear strain tensor eij is expressed in terms of the displacement vectors (ui ) such that eij=(ui,j+uj,i)/2. The ®nal form of the equilibrium equation is derived by substituting Eq. (17) into Eq. (16) and expressing the linear elastic tensor with equivalent displacement vector as   1 …f † Cijkl …uk,l ‡ ul,k † ÿ bij …T ÿ T0 † ÿ wg jdij ,,j ‡ rbi ˆ 0: …18† 2 3.3. Energy conservation equation To derive the energy conservation equation for a porous medium, the energy conservation equations for solid and ¯uid phases are derived independently and then these two derivations are summed up to yield the

568

K.M. Neaupane et al. / International Journal of Rock Mechanics and Mining Sciences 36 (1999) 563±580

  …C† ds C …v† dT 1 @sij ˆ ÿ T dt dt r @T

!

®nal form of the energy conservation equation. At ®rst, the energy conservation equation for the solid phase is derived. From the energy conservation law, the rate of change of internal energy may be expressed as the sum of stress power plus the heat added to the continuum [12]. When the stress power is not fully recoverable, it is sometimes helpful to split stress tensors into two components i.e. conservative or recoverable and dissipative stress components. It is not at all clear that such a separation can always be made. In ¯uid, for example, it is assumed that the recoverable stress component is equal to hydrostatic stress and the remainder of Cauchy stress is the dissipative stress component as in

is the rate of heat ¯ux per unit area by conwhere h…x† i duction and q (T) is called the radiant heat constant. But from Eq. (17), i.e. the Duhamel±Neumann equation of thermo±elasticity and Eq. (24), the energy conservation law without the e€ect of viscous dissipation becomes

…C† s…D† ij ˆ sij ÿ sij ,

rC …v†

…19†

where s(D) and s(C) are dissipative and conservative ij ij stress components respectively. The energy equation may be written as du 1 1 dq ˆ s…C† e_ ij ‡ s…D† e_ ij ‡ , dt r ij r ij dt

…20†

where u is called the speci®c internal energy, eij is the strain tensor and q is the heat in¯ux per unit mass into the continuum. It is supposed that only the dissipative stress component, s(D) ij , contributes to the internal entropy production as given by Clausius±Duhem inequality principle ds 1 …D† 1 dq ˆ s e_ ij ‡ , dt rT ij T dt

…21†

where s is the entropy density and T is the absolute temperature. A potential `f' called the Helmholtz free energy is now introduced such that f ˆ u ÿ sT. The derivative of the Helmholz free energy and equation of equilibrium (15) gives rise to the following expression: df ˆ

1 …C† s deij ÿ sdT: r ij

…22†

For thermoelasticity the state variables would be the deformation gradient (eij) and absolute temperature (T ). Thus for an ideal thermoelastic behavior, both the Helmholtz potential ( f ) and speci®c internal energy (u) are a function of absolute temperature (T ) and deformation gradient (eij) such that f ˆ f…T,eij † and u ˆ u…T,eij † respectively. From the partial derivative of Eq. (22), expressing speci®c internal energy in terms of speci®c entropy and thereby speci®c heat at constant volume, the equation for the rate change of entropy reduces to

@eij : eij dt

…23†

By equating Eq. (21) and Eq. (23), after some manipulation, we get ! …C† @s dT @eij ij …D† …T† …24† ˆ ÿh…x† ‡T rC …v† i,i ‡ sij eij ‡ rq @T dt dt

dT @ eij …T† : ˆ ÿh…x† ÿ bij T i,i ‡ rq dt dt

…25†

Finally, the energy conservation equation for a solid constituent takes the form   …s† …s† …s† …s† …v†…s† dT ‡ vi T,i …1 ÿ n†r C dt ÿ …1 ÿ n†bij T …s† ˆ ÿ…1 ÿ n†h…x†…s† i,i

@e…s† ij ‡ r…s† q…T†…s† , @t

…26†

where superscript `s' is axed to indicate that the properties are associated with a solid matrix. To arrive at the conservation law for ¯uid, we begin with a constitutive equation in which the recoverable stress component is equal to the hydrostatic stress components i.e. s(C) ij =ÿPdij. With simple manipulation the constitutive relation for ¯uid can be readily expressed as s…C† deij ˆ …ÿP †dij deij ˆ ÿPdekk ,

…27†

where ekk is the volumetric strain. We now consider speci®c entropy as a function of absolute temperature `T' and volume `V' such that s=s …T, V †. The expression for rate of change of speci®c entropy is then deduced following the procedure that was used to derive Eq. (24),     ds C …v† dT 1 @ P dV ˆ ‡ , …28† T dt dt r @ T V dt  where (@ P @ T)V means that the partial derivative is taken at constant volume (V ). Equating Eq. (22) and Eq. (28), we arrive at ! …T† dT b …D† …T† ˆ ÿh…x† ‡ T …P† vi,i , …29† rC …v† i,i ‡ sij eij ‡ rq dt b where b (T) and b (P) are  expansivity  and compressibility of water such that @P @ T ˆ b…T† b…P† is constant. The

K.M. Neaupane et al. / International Journal of Rock Mechanics and Mining Sciences 36 (1999) 563±580

energy conservation equation for ground water now becomes   …f† …f† …f† …r† …f† …v†…f† dT ‡ vi T,i nS r C dt ! …T† b …f† …T†…f† ÿ nS …r† T …f† …P† v…f† …30† ˆ ÿnS …r† h…x†…f† i,i i,i ‡ r q b in which superscript `f' is used to describe properties related to ¯uid constituent. If the movement of ¯uid through porous media is suciently slow and the surface areas of all phases are suciently large, it is reasonable to assume that thermal equilibrium among phases is achieved. Thus assuming T=T (s)=T (f), the energy conservation equations for solid and water, i.e. Eq. (26) and Eq. (30), can be united to yield the conservation equation at the given temperature. With the assumption that Fourier's law is valid for heat ¯ux, )(f) )(s) =ÿK(T)(f) T,j and h(x =ÿK(T)(s) T,j, the coni.e. h(x i ij i ij stitutive equation for energy conservation is ®nally reduced to dT …T†…m† ‡ nS …r† r…f† C …v†…f† v…f† T,ji i T,i ÿ Kij dt ! b…T† @ Ui,i …r† ÿ q …T† ˆ 0, …31† ‡nS T …P† k…y†ij h,ji ‡ …1 ÿ n†bT @t b

…rC …v† †…m†

where

¯ow problems. It should, however, be noted that pore pressure distribution is unknown at the moving hydraulic boundaries. In such cases, it is reasonable to use a total head expression in which the reference point for the elevation is the lowest in the domain. Using this concept, total head may be de®ned as the sum of elevation and pressure head and Eq. (18) can then be written as:   1 …f † Cijkl …uk,l ‡ ul,k † ÿ bij …T ÿ T0 † ÿ wg jdij 2 ,j …33† ‡ r …s† bi ˆ 0, where the last term in the equation is a gravity contribution (external body force). The continuity ¯ow equation is rewritten as …T† ÿfnS …r† gr…f† 0 b

‡ r…f† C…j†


ÿ @T …P† …f † @j ‡ fnS …r† gr…f† r g o b @t @t

@j @ Ui,i ÿ r…f† k…y †ij h,ji ˆ 0, ‡ r…f† S …r† @t @t

…34†

and the energy conservation equation is …rC …v† †…m†

dT ‡ nS …r† r…f† C …v†…f† v…f† i T,i dt !

ÿ K…T†…m† T,ji ij ‡ …1 ÿ n†bT

…rC …v† †…m† ˆ nS …r† r…f† C …v†…f† ‡ …1 ÿ n†r…s† C …v†…s† ,

569

‡ nS …r† T

b…T† b…P†

k…y †ij h,ji

@ Ui,i ÿ q …T† ˆ 0 @t

…35†

ˆ nS …r† K…T†…f† ‡ …1 ÿ n†K…T†…s† and K…T†…m† ij ij ij q …T† ˆ r…f† q…T†…f† ‡ r…s† q…T†…s† :

…32†

Parameter K refers to the coecient of heat conduction. Equation (31) is the energy conservation equation in which stress-deformation and ¯uid ¯ow are considered. The ®rst term on the left hand side of the equation shows a dependency on energy and the second term shows energy change due to heat convection. The next three terms express energy change due to heat conduction, pore water pressure and reversible energy change caused by solid deformation respectively. The last term q …T† is the total heat ¯ux added to the continuum.

4. Governing equations Eqs. (14), (18) and (31) represent the governing equations for fully coupled thermal±mechanical-¯uid

5. Numerical solution procedure 5.1. Numerical formulation The numerical solution procedure employs a numerical model based on the ®nite element method (FEM). It is based on Galerkin's formulation, a weighted residual method in which trial functions themselves are chosen as weighting functions. The unknowns are displacement, head and temperature. The procedure essentially consists of two steps: ®rst, an arbitrary volume integration of all three governing Eqs. (32)± (35), is taken and then trial functions are de®ned for the unknowns Ui (displacement), h (head) and T (temperature) leading to a weighted residual approximation based on the Galerkin method. The ®nite element discretization of the displacement (U ), head (h ) and temperature (T ) is de®ned by the piece wise approximation of the trial functions in the equilibrium equation as

570

K.M. Neaupane et al. / International Journal of Rock Mechanics and Mining Sciences 36 (1999) 563±580

9 m …t† …1† m …t† > > ˆ N N…1† U U > m i m i > > > mˆ1 > > > X = …2† m …t† …2† m …t† N m hi ˆ Nm hi hh  : > mˆ1 > > > > X > > m …3† m > N…3† Th  m Ti …t† ˆ Nm Ti …t† > ; X

Uh 

Fin ˆ

…

… …36†

‡

Sa

… Hnm ˆ

Omitting the intermediate details, we arrive at the following system of ordinary di€erential equations expressed in matrix form as 2 3 2 m 3 0 0 0 0 U1 …t† 7 60 7 @ 6 Um 0 0 0 6 1 7 6 2 …t† 7 4 Anm A2nm Enm Dnm 5 @t 4 hm …t† 5‡ J1nm J2nm 0 Mnm T m …t† 2 11 32 m 3 2 1 3 U1 …t† C1nm G1nm Knm K11 Fn nm 2 76 Um …t† 7 6 F2 7 6 K11 K11 C2 G 2 nm nm nm nm 7 ˆ 6 n 7, …37† 6 1 76 4A A2nm Hnm 0 54 hm …t† 5 4 Qn 5 nm T m …t† J1nm J2nm Lnm Vnm Pn where the symbolic matrix notations represent the following in integral forms: … i …f† …r† …1† N…2† …38† Anm ˆ n r S Nm,i dVa Va

Enm ˆ

Va

Jinm ˆ

… Va

…f† …P† …r† …f† N…2† ‡ r…f† C…j†gN…2† n fr0 nS r gb m dVa

Va

Kik nm

ˆ

Cinm ˆ

Ginm

Va

… Va

… Va

… ˆ

Va

…f† …r† …f† …T† …3† N…2† n r0 nS r b Nm dVa

…1† N…3† n …1 ÿ n†bTNm,i dVa

… Mnm ˆ

Qn ˆ ÿ

…v† …m† …3† N…3† Nm dVa n …rC †

Sa

Va

N…3† n,i

… … Va

 …f† dSa N…2† n q

… ‡

(

…

Va

…46†

…r†

ÿ nS T

Va

…47†

…f† N…2† n r cdVa

b…T† b…P†

Sa

…48†

) kij N…2† m,j dVa

…49†

  …r† …f† …v†…f† …f† N…3† r C v nS N…3† n,i i m,i dVa ‡

…T†…m† …3† N…3† Nm,j dVa n,i Kij

… Pn ˆ ÿ

… Dnm ˆ ÿ

…

Vnm ˆ

Va

 …s† bi dVa N…1† n r

…2† …f† N…2† n,i r kij Nm,j dVa

Va

Lnm ˆ

…

ÿ
…3† N…1† n,j ÿ bTo dij Nm dVa

Va

mˆ1

…

^ N…1† n ti dSa ‡

…T† N…3† n Q dSa

…50†

… ‡

Va

…T† N…3† n q dVa

…51†

…39†

…40†

…41†

…42†

…1† N …1† n,j Cijkl N m,l dVa

…43†

…f† …2† N…1† n,j …ÿwdij g †Nm dVa

…44†

…3† N…1† n,j …ÿbdij †Nm dVa

…45†

On the left hand side of Eq. (37), the displacement (Ui), ¯uid pressure (h ) and the temperature (T ) are nodal unknowns at the nodal points. The matrices represented by letters K, C and G are related to displacement (Ui), displacement±head coupling (Ui±h ) and displacement±temperature coupling (Ui±T ) respectively. Similarly A, (E+H) and D are associated with head±displacement coupling (h±Ui), pressure head (h ) and head±temperature coupling (h±T ), whereas J, L and (M+V) are associated with temperature±displacement coupling (T±Ui), temperature±head coupling (T± h ) and temperature respectively, all in ®nite element terms. The vectors on the right hand side of the equilibrium equations are known quantities. We now introduce the transient character of the coupling phenomena, adding an extra dimension to the problem±time e€ects. In actual ®eld simulations, as for example of an underground radioactive waste disposal or in-ground LNG storage system, time e€ects may be partly attributed to thermal and hydraulic diffusion and partly due to creep. But for freezing and thawing laboratory simulations, the creep phenomenon has been neglected. The numerical solution can use a time-step scheme for which a given quantity (U, h and

K.M. Neaupane et al. / International Journal of Rock Mechanics and Mining Sciences 36 (1999) 563±580

T ) at two di€erent time instants t and t+Dt are interpolated linearly f…t ‡ yDt † ˆ …1 ÿ y †f…t† ‡ yf…t ‡ Dt †

…52†

with y 2 ‰1,0Š

The continuity ¯ow equation in time domain becomes  Aknm um k …t ‡ Dt † ‡ Enm ‡ yDtHnm hm …t ‡ Dt † ‡ Dnm T m …t ‡ Dt † ˆ Aknm um k …t†  ‡ Enm ÿ yDtHnm hm …t† ‡ Dnm T m …t† ‡…1 ÿ y †DtQn …t† ‡ yDtQn …t ‡ Dt †

…53†

and the energy balance equation is written as

571

m Jknm um k …t ‡ Dt † ‡ yDtLnm h …t ‡ Dt †

 ‡ Mnm ‡ yDtVnm T m …t ‡ Dt † m … † ˆ Jknm um k …t† ÿ 1 ÿ y DtLnm h …t†

 ‡ Mnm ÿ …1 ÿ y †DtVnm T m …t† ‡ …1 ÿ y †DtPn …t† ‡ yDtPn …t ‡ Dt †:

…54†

5.2. Initial and boundary conditions The ¯ow domain is discretized in a graded size mesh of Lagrangian elements as shown in Fig. 2(a). All but one are quadrilateral elements. The types of element used for displacement, temperature and head are described in Fig. 2(b) and Fig. 2(c) for quadrilateral and triangular elements respectively.

Fig. 2. Details of ®nite element mesh.

572

K.M. Neaupane et al. / International Journal of Rock Mechanics and Mining Sciences 36 (1999) 563±580

Table 1 Input parameter values used in calculation for rocks (Temperature ÿ20±+208C) Parameters

Tu€ (Tage stone)

Sandstone (Sirahama )

Young's modulus (E ) in kgf/cm2 Poisson's ratio (n ) Density r (s) in g/cm3 Porosity (n)% Coecient of permeability (k ) in cm/s Degree of Saturation S(r) in % Thermal conductivity (K ) Cal/cm s.8C Speci®c heat C (v)(s) in Cal/g 8C Thermal expansivity (a ) in /8C

1.0  104 0.25 1.83 22.4 1.0  10ÿ5 100 2.5  10ÿ3 0.195 9  10ÿ6

4.02  104 0.37 2.41 13.0 1.42  10ÿ8 100 5.28  10ÿ4 0.195 8.8  10ÿ6

5.2.1. Initial condition

5.3. Material properties input

9 Ui …x,y,t † ˆ Ui …x,y,0† on R = h…x,y,t † ˆ h…x,y,0† on R ; T…x,y,t † ˆ T…x,y,0† on R

The input material properties of the rock used in the code have been presented in Table. 1 and the properties of water/ice used have been tabulated in Table 2. Young's modulus and Poisson's ratio of the rock type speci®ed have been obtained from uniaxial compression tests on dry specimen samples. Porosity has been derived from density measurements. Thermal conductivity and speci®c heat data were obtained by thermal needle experimental procedure. Thermal expansivity coecient data for the water/ice used are given in Fig. 3. Values for the coecient of permeability have been taken from the literature [18].

…55†

5.2.2. Boundary condition Displacement: Ui …x,y,t † ˆ U^ i …x,y,t † Traction: sij …x,y,t †nj ˆ T^ i …x,y,t † Hydraulic: h…x,y,t † ˆ h^…x,y,t † Flow rate: ÿ k…y †ij h,j ni ˆ Q…x,y,t † Temperature: ÿ K…Tm† T,j ni ˆ Q…T† …x,y,t † ij

9 on Bu > > > > > > on Bs > > = on Bh > > > on BQ > > > > > on B ;

6. Stability and accuracy of the numerical scheme

q

…56† where (x, y, z ) is position vector, ni is the unit normal vector and UÃ, hÃ, Tà are known displacement, head and surface traction respectively; Qà and QÃ(T) are prescribed ¯ow rate and heat ¯ow respectively. R stands for the continuum region enclosed by a boundary B and subscripts s, h, Q and q are used to describe stress, hydraulic, ¯ow rate and thermal boundary conditions respectively.

A two-dimensional computer code HEAT has been developed in plane stress poro±elastic assumption in which deformation, ¯uid ¯ow and heat ¯ow are fully coupled. As the main objective of this research is to solve problems involving freeze±thaw, it is necessary to address the moving boundary problem at which phase

Table 2 Input parameter values for ice/water used in calculation Parameters

Values

Young's modulus (E ) in kgf/cm2 Compressibility of water in cm2/kgf Thermal expansivity of water Latent heat of ice(water) in Cal/cm3 Poisson's ratio (n) Density of water r (f) in g/cm3 Density of water/ice at 08C in g/cm3 3 Density of ice r(f) i in g/cm Speci®c heat of water C (v)(f) in Cal/g 8C Speci®c heat of water/ice at 08C in Cal/g 8C Speci®c heat of ice in Cal/g 8C

9.0  104 4.9  10ÿ5 See (Fig. 3.) 80.00 0.350 1.000 0.9585 0.917 1.000 0.725 0.450

Fig. 3. Thermal expansivity of water as a function of temperature.

K.M. Neaupane et al. / International Journal of Rock Mechanics and Mining Sciences 36 (1999) 563±580

573

Fig. 5. A comparison of numerical and experimental simulation results with and without latent heat (rock: Tu€, distance from heat source=4 cm).

Fig. 4. Program ¯ow chart.

transition takes place. Further, there may be no sharp boundary between two phases and the phase transition may not be distinct. Established mathematical formu-

lation can solve simple moving boundary problems but fails to explain the nature of the mushy region [2]. In this paper, mathematical formulation of moving fronts has not been addressed as such, but for every element described in time±space domain, freezing index (IPH) is set at ÿ1, 0 or 1 depending on whether the average temperature of the element is below freezing, at zero8C or above zero8C. If the freezing index of an element is ÿ1, a latent heat term is added to the heat ¯ux term of the energy balance equation. Material properties that are a function of temperature are changed accordingly and equilibrium conditions are updated, as shown in the program ¯ow chart (see Fig. 4). The e€ect of the latent heat term addition to the energy balance equation can be clearly seen in Fig. 5. The simulation is more realistic with the addition of the

Fig. 6. One dimensional consolidation problem.

574

K.M. Neaupane et al. / International Journal of Rock Mechanics and Mining Sciences 36 (1999) 563±580

Fig. 7. Thermoelastic problem of a thick walled cylinder.

latent heat to energy balance equation and is closer to the experimental results (see Fig. 10 also). Details on this have been presented in the temperature transfer section elsewhere. One of the fundamental diculties associated with numerical solutions of THM coupling systems is that it must be possible to compare a numerical solution with a reliable benchmark solution. Benchmark solutions are, however, available only for very simpli®ed cases. Such closed-form solutions are viewed as general cases of speci®ed problems and hence are important to give a good qualitative view. The performance of a developed code can not be tested with an analytical benchmark solution as such a solution for coupled problems in which stress, pore pressure and temperature are mutually linked, is not available. Therefore the code has been veri®ed under the following conditions:

. a one dimensional elastic consolidation problem, . the thermoelastic problem of a thick walled circular cylinder subjected to a uniform temperature gradient and . a heat transfer problem. It has been found that the analytical solution and solution from developed numerical code are in good agreement (Fig. 6±8). This veri®es the fundamental functions contained in the numerical code HEAT. Stability of such a numerical scheme has always been a matter of concern. For high permeable media, the ¯uid velocity vector is highly oscillative but for small permeability values of the order of 10ÿ5 cm/s to 10ÿ8 cm/s, this numerical scheme yields satisfactory results.

Fig. 8. A problem of simple heat transfer only.

K.M. Neaupane et al. / International Journal of Rock Mechanics and Mining Sciences 36 (1999) 563±580

Fig. 9. (a) Layout of the experimental set-up. (b) Plan view of the specimen with strain gauges and temperature sensors.

575

576

K.M. Neaupane et al. / International Journal of Rock Mechanics and Mining Sciences 36 (1999) 563±580

Fig. 10. Temperature transfer: comparison of numerical and experimental results (rock: Tu€).

Fig. 9(b). The experiments were designed to cover room temperature as well. The tests were performed on fully saturated specimen. The specimen was con®ned within a thermostat box so as to closely simulate the adiabatic condition. A circulation type low-temperature thermostat bath (TRL-N30L) was used for source temperature variation from ÿ208C to +208C through brine circulation through the center of the circular hole (f=4.6 cm) of the specimen. A thin cylinder of copper covered with aluminum foil separated the brine solution from the rock specimen boundary. The temperature of the thermo-static bath was maintained at ÿ208C and specimen was subjected to freezing for 72 h to achieve as low and stable a temperature as possible. Measurements were taken every minute. The low temperature circulation was then cut o€ and the frozen specimen was allowed to thaw for another 72 h. An automatic data logger was used for data recording. Results obtained from the experiment are presented in the following chapter with subsequent explanation.

7. Freeze±thaw experiment 7.1. Specimen samples In order to examine the validity of this numerical model-based code, freezing and thawing experiments were carried out. Results of the laboratory experiments and of the thermo±mechanical analysis from the code are compared to demonstrate the function of the developed code. Measurement from the laboratory tests includes strain along three directions and temperature, both at pre-de®ned locations. Two major tests were performed; one each for Tu€ and Sirahama sandstone. The Sirahama sandstone name derives from the Sirahama district in Japan from where it was sampled. The Tu€ (Tage stone) was sampled from Tochigi prefecture in Japan.

8. Comparison of results and discussion 8.1. Temperature transfer A comparison of the simulated temperature and experimental temperature results for Tu€ during freezing is shown in Fig. 10. This comparison has been done for radial distances of 1 and 4 cm from the heat source. Numerical estimations are shown by dashed lines while experimental results are shown with ®rmed lines. For both numerical and experimental cases, the

7.2. Experimental procedure The experimental set up for the freezing and thawing tests performed consisted essentially of a number of temperature sensors and strain gauges placed along pre-de®ned directions at speci®ed distances from the center of the heat source. A prismatic specimen of 30  45  15 cm was prepared with a hole of 4.6 cm diameter at its center. Then the temperature sensors and strain gauges were ®xed on the smooth surface of the specimen at pre-de®ned distances and directions and covered with water-proof resin before it was submerged into water for 3 days to fully saturate it. The complete experimental layout is shown in Fig. 9(a) and the position and orientation of the temperature sensors and strain gauges on the rock specimen are shown in

Fig. 11. Temperature transfer: comparison of numerical and experimental results (rock: Sirahama sandstone, distance from heat source=1 cm).

K.M. Neaupane et al. / International Journal of Rock Mechanics and Mining Sciences 36 (1999) 563±580

577

Fig. 12. Mechanism of deformation induction from phase change and moving boundary.

temperature gradients are steep at ®rst but becomes gentler to approach the steady state condition gradually. It can be seen, in this case, that the agreement between simulated transient pro®les and the experimental results is fairly good. Referring to Fig. 5, it shows the simulated temperature pro®les at a distance of 4 cm from the heat source. The simulated pro®le

without the e€ect of the latent heat shows a steeper gradient at ®rst and reaches the steady state condition earlier; this is not realistic as far as its comparison with the experimental result is concerned. Fig. 11 shows the simulated temperature pro®les alongside the experimental data for Sirahama sandstone during freezing and thawing. The results are for

Fig. 13. Temperature versus radial strain: results from numerical simulation and laboratory experiment (rock: Tu€, distance from heat source=1 cm).

Fig. 14. Temperature versus radial strain: results from numerical simulation and laboratory experiment (rock: Tu€, distance from heat source=4 cm).

578

K.M. Neaupane et al. / International Journal of Rock Mechanics and Mining Sciences 36 (1999) 563±580

generally be agreed that the numerical scheme can satisfactorily predict temperature transfer from freeze± thaw experiments. 8.2. Induced deformation

Fig. 15. Numerical simulation of thermal loading±unloading (rock: Tu€, distance from heat source=1 cm).

a radial distance of 1 cm from the heat source. A very good match between the simulated and experimental results can be seen at the early stage of freezing but the temperature from the laboratory experiment reaches the steady state while the simulated temperature still shows a small gradient. During thawing, nonlinearity can be seen at around zero 8C due to the phase change. The simulated and experimental results match well during the later stage of thawing as both reach the steady state condition. Thus the numerical simulation and experimental results exhibit similar trends as far as temperature transfer during freezing and thawing is concerned. From these results, it may

Fig. 16. Temperature versus radial strain: results from numerical simulation and laboratory experiment (rock:Sirahama sandstone, distance from heat source=1 cm).

The complete mechanism of the freezing and thawing process and induced deformation from it, may be explained with a simple scheme (see Fig. 12). It is logical that as the temperature decreases, the specimen is subjected to compression at ®rst because it is at the freezing stage. But if the temperature decreases below freezing, it should be accompanied with greater dilation as a result of an increase in volume due to phase change. Initially, the freezing boundary is close to the heat source. Therefore, at the observation point, as shown in Fig. 12, the strain gauge is subjected to contraction. As the moving freezing boundary approaches the observation point, the strain gage at the observation point experiences expansion because of volume increase. As the boundary moves further, the observation point may be subjected to compression again due to freezing. The characteristic curves between temperature and radial strain obtained from laboratory experiments as well as numerical simulation may be explained by this concept. 8.2.1. Temperature versus strain As it is inconvenient to measure induced stresses directly from laboratory experiments, e€orts were made to compare radial strains readily measurable from the experiments. Fig. 13 and 14 illustrate induced radial strain resulting from simulation and the experiment as a function of temperature during freezing for the Tu€. Fig. 13 shows results for a point 1 cm from the heat source and Fig. 14 shows the same for a point 4 cm from the heat source. Non-linear elastic properties of ice make the temperature±strain relationship highly non-linear and complex. In both cases, when the temperature decreases from the initial condition, the negative radial strain increases as the material is under compression. As water in pore spaces changes phase at around 08C, the non linear material behavior can be seen to develop during that time. When the temperature drops below freezing, positive strain starts to develop because of the expansion of water on freezing. Though the experimentally measured radial strains and simulated radial strains do not match with the simulated radial strains at places, di€erences are small and trends are similar. Fig. 15 is the result of the numerical simulation of the freezing and thawing process for the Tu€, analogous to the loading±unloading mechanism in a nonlinear elastic continuum. It is evident that the charac-

K.M. Neaupane et al. / International Journal of Rock Mechanics and Mining Sciences 36 (1999) 563±580

teristic curves from freezing and thawing are in good agreement except at the freezing boundary. This deviation shows that the boundary between phases may not be sharp and may be attributed to the latent heat involved during phase change. This is also shown in Figs. 14 and 16 (encircled region). Fig. 16 displays induced radial strain as a function of temperature for the Sirahama sandstone during a cycle of freezing and thawing. The comparison between simulated and experimental results has been done for a point at a radial distance of 1 cm. The laboratory test results show that a large deformation induced during thermal loading (freezing) would not be recovered on thawing. Residual strain of the order of more than 0.3% may be large enough to induce plastic deformation. Prediction from the numerical code, however, shows considerable variation from experimental results as it is based on a poro±elastic assumption (see Fig. 16). Therefore for a strain value as large as those obtained from the freeze±thaw experiment with Sirahama sandstone, plasticity may be the factor to be reckoned with.

9. Conclusion The work presented in this paper describes application of a numerical model of thermo±hydro±mechanical coupling to the simulation of laboratory freezing and thawing experiments. A theoretical formulation that accommodates a simple linear stress±strain constitutive relationship has been presented and a twodimensional (plain stress condition) numerical modeling was performed based on the ®nite element method applied to thermo±poro±elasticity. The freezing and thawing experiments described here on rock specimens and their evaluation by model computations have shown that with numerical simulation, a relatively good prediction can be made of the temperature transfer and deformation behavior if the materials under test do not fail. This is because deformation thus computed is purely based on elastic considerations and thus the ®nite element code fails to predict large deformation resulting from material failure. The developed numerical model based code is thought to be applicable for a number of coupled engineering problems namely: LNG storage in frozen ground, sinking of shafts in frozen ground or the same by freezing method in saturated/unsaturated loose ground, natural slope deterioration, weathering and ¯oor heaving in lower temperature zones, etc. In addition, the code is also applicable for coupled problems associated with nuclear waste repository, reservoir engineering, energy extraction, etc. However, the code may not be applicable for problems where the material

579

under consideration fails, as mentioned above. As far as future research is concerned, irreversible and plastic e€ects should be taken into account. Hence an elasto± plastic formulation may be required to fully describe the freezing±thawing mechanism.

Acknowledgements This research was carried out in the Rock Mechanics Laboratory of Saitama University. The authors gratefully acknowledge the contribution of graduate students M. Kakizawa, S. Kurosawa and H. Watanabe in carrying out experimental investigations as well as computer programming.

References [1] Inada Y, Yokata K. Some studies of low temperature rock strength. Int J Rock Mech Min Sci & Geomech 1984;21(3): 145±153. [2] Crank J. Free and Moving Boundary Problems. UK: Oxford Science Publication, Oxford University Press, 1984. [3] Amy MO, Rouset G. Thermo±hydro±mechanical modeling of an underground radioactive wastes disposal. In: Smith IM, editor. Numerical Methods in Geotechnical Engineering. Rotterdam, Netherlands: A.A. Balkema, 1994. [4] Marsily G. An overview of coupled process with emphasis on geohydrology. In: Tsang C-F, editor. Coupled Processes Associated with Nuclear Waste Repositories. Orlando, FL, USA: Academic Press Inc, 1987. [5] Hart RD, St. John CM. Formulation of a fully coupled thermal±mechanical±¯uid ¯ow model for non-linear geologic systems. Int J Rock Mech Min Sci & Geomech 1986;23(3): 213±24. [6] Thunvik R, Braester C, Hydrothermal conditions around a radioactive waste repository. In: Lutze W, editor. Scienti®c Basis for Nuclear Waste Management V. Material Research Society, Symposium proceedings, volume II. New York, USA: Elsevier Science Publishing Co., Inc, 1982. [7] Aboustit BL, Advani SH, Lee JK. Variational principles and ®nite element simulations for thermo±elastic consolidation. Int J Num Anal Methods Geomech 1985;9:49±69. [8] Paaswell RE. Temperature e€ects on clay soil consolidation. J. of Soil Mechanics and Foundation Div., Proceedings of ASCE, May 1967; 93 (SM3). [9] Booker JR, Savvidou C. Consolidation around a point heat source. Int J Num Anal Methods Geomech 1985;9:173± 84. [10] Seneviratne HN, Carter JP, Booker JR. Analysis of fully coupled thermo±mechanical behavior around a rigid cylindrical heat source buried in clay. Int J Num Anal Methods Geomech 1994;18:177±203. [11] Giraud A, Rouset G. Thermoelastic and thermoplastic response of a porous space submitted to a decaying heat source. Int J Num Anal Methods Geomech 1995;19:475±495 . [12] Hart RD. A fully coupled thermal±mechanical±¯uid ¯ow model for non-linear geologic systems, Ph.D. thesis, Univ. of Minnesota, 1981. [13] Ohnishi Y, Kobayashi A, Thermal±hydraulic±mechanical coupling analysis of rock mass. In: Hudson JA, editor.

580

K.M. Neaupane et al. / International Journal of Rock Mechanics and Mining Sciences 36 (1999) 563±580

Comprehensive Rock Engineering. Analysis and Design Method, Vol. 2. UK: Pergamon Press, 1993. [14] Biot MA. General theory of three dimensional consolidation. J Appl Phys 1941;12:155±64. [15] Malvern LE. Introduction to the Mechanics of a Continuous Medium. NJ: Prentice-Hall, 1977. [16] Chen Wai-Fah, Saleeb Atef F. Constitutive equations for engin-

eering materials. In: Elasticity and Modeling, vol. 1. Netherlands: Elsevier Amsterdam, 1994 Second ed. [17] Gould PL. Introduction to Linear Elasticity. USA: SpringerVerlag Inc, 1983. [18] Xue Z. Permeability and microgeometry of sandstones under hydrostatic pressure, Ph.D. Thesis, Faculty of Engineering, Hokkaido University, Japan, 1992.