A non-linear analysis of non-isothermal wave propagation in linear-elastic fluid-saturated porous media

International Journal of Plasticity 18 (2002) 313±344 www.elsevier.com/locate/ijplas

A non-linear analysis of non-isothermal wave propagation in linear-elastic ¯uid-saturated porous media A. Gajo Dipartimento di Ingegneria Meccanica e Strutturale, UniversitaÁ di Trento - Via Mesiano, 77, 38050 Trento, Italy Received in ®nal revised form 9 December 1999

Abstract The non-isothermal dynamic behaviour of saturated porous media is analysed numerically employing the ®nite element method and taking energy convection due to large pore ¯uid displacements into account. A di€erent pore ¯uid reference temperature is introduced in order to allow properly for heat convection: this concept is usually neglected in the literature and is discussed and analysed herein. The numerical procedure is validated in a simple problem of hot ¯uid injection in a steady seepage ¯ow and by comparing the numerical results, neglecting energy convection, with those obtained with a novel solution of the linearised equations, presented herein, which is based on the transfer functions and Fourier transforms method. Finally, the e€ects of energy convection in wave propagation are analysed: in a pervious porous medium the ¯ux of energy due to energy convection is much greater than the one due to heat conduction; in any case, wave propagation can be considered completely adiabatic even when energy convection is taken into account. Thus the validity of the results presented in the literature and based on the linearised theory is demonstrated. # 2002 Elsevier Science Ltd. All rights reserved. Keywords: Saturated porous material; Non-isothermal conditions; Wave propagation; Finite elements; Analytical solution

1. Introduction The theoretical framework for the analysis of isothermal wave propagation in ¯uid-saturated elastic porous media was ®rst proposed by Biot (1941, 1956, 1962a,b). He examined both the high and the low frequency limits and showed the existence of two longitudinal waves and one shear wave, which are dispersive and dissipative. Viscous coupling and frequency content have important e€ects on the 0749-6419/02/$ - see front matter # 2002 Elsevier Science Ltd. All rights reserved. PII: S0749-6419(00)00100-5

314

A. Gajo / International Journal of Plasticity 18 (2002) 313±344

resulting velocities of these waves: when viscous coupling is small and the frequency range is high, relative movements between the two phases may occur; these movements are in phase for the fastest wave and out-of-phase for the slowest one. Conversely, when the permeability is low and low frequencies are involved, relative movements are prevented and, for wave propagation, the saturated porous medium behaves as a one-phase medium. In particular, the slowest longitudinal wave is only a real wave when viscous coupling between the two phases is minimal and the frequency range is high, otherwise it loses its wave character and is reduced to a di€usion phenomenon with a very high attenuation (Biot, 1956; Garg et al., 1974). Pecker and Dereziewicz (1973) and Bowen and Chen (1975) later performed a linearised analysis of non-isothermal wave propagation. In particular, Pecker and Dereziewicz (1973) showed that under non-isothermal conditions there are four longitudinal waves and one shear wave, and that these waves are dissipative and dispersive. Moreover, propagation of the ®rst two longitudinal waves is essentially adiabatic in the high-frequency regime and is dominated by interphase heat transfer in the low-frequency regime, with high dispersion and dissipation. Finally, in the framework of the linearised theory, Bowen and Reinicke (1978) analysed which waves have a non-dispersive low-frequency approximation. These theoretical ®ndings are very useful for interpreting dynamic measurements which, from the conclusions of Pecker and Dereziewicz (1973), provide the adiabatic elastic properties of the solid skeleton. The theoretical results have been veri®ed experimentally by many authors (Berryman, 1980; Plona, 1980; Hamdi and Taylor Smith, 1982; Johnson et al., 1982; Van Der Grinten et al., 1985, 1987; Gajo et al., 1997). Although further experimental evidence would be needed, some inconsistencies between the linearised theory and the experimental ®ndings appear to exist when the viscous coupling is small, i.e. when relative movements between the two phases are possible; conversely, the consistency between the theory and the experiments is excellent when viscous coupling is high and relative movements are prevented (Gajo et al., 1997). On the other hand, it has been shown that the interpretation of dynamic tests is extremely sensitive to even the smallest approximations, especially if the sti€ness of the solid skeleton is low (Gajo and MongiovõÁ, 1994); the use of a linearised theory in which the energy convection due to ¯uid ¯ow is neglected might consequently a€ect the interpretation of experimental ®ndings and explain the apparent inconsistencies between the theory and the experiments. The main purpose of this work is to evaluate the role of energy convection due to pore ¯uid ¯ow during wave propagation in ¯uid-saturated porous media and to ascertain any limits of the usual linearised analyses. Due to the non-linearity related to energy convection, the problem is solved numerically. To this purpose, a ®nite element method able to deal with the dynamic-thermo-hydro-mechanical coupling of ¯uid-saturated porous media, in any range of frequency and on the assumption of large pore ¯uid displacements, is developed and validated. The present work is composed of ®ve parts: ®rst the ®eld equations are brie¯y reiterated, discussing the need to introduce a new reference temperature for the pore ¯uid in order to take the e€ects of heat convection suitably into account. The second

A. Gajo / International Journal of Plasticity 18 (2002) 313±344

315

part presents a novel solution in time domain for the linearised problem based on the transfer functions and Fourier transforms method. The third part describes the non-linear ®nite element method for solving the dynamic-thermo-hydro-mechanical coupling in ¯uid-saturated porous media, developed for large pore ¯uid displacements. In the fourth part, the numerical procedures are validated against the solution in time domain described in the second part and against simple problems of heat convection in a steady seepage ¯ow. Finally, the e€ects of energy convection due to ¯uid ¯ow in wave propagation are analysed and discussed. Although this study is mainly devoted to the analysis of non-isothermal wave propagation, the topics dealt with and the numerical procedures developed and validated are of interest in a number of engineering and environmental problems, where both high and low frequencies could be involved, e.g. nuclear waste disposal, extraction of geothermal energy, mining and petroleum engineering, intrusion of hot magma in wet rocks, or response to frictional heating in fault zones. The interest in this kind of problem is con®rmed by the number of recently-published studies on the thermomechanical constitutive modelling of porous media (e.g. Saada et al., 1996; Hueckel et al., 1998), which could be employed in non-isothermal ®nite element calculations. 2. Field equations 2.1. Momentum balance equations The momentum balance equations of a saturated isotropic porous medium are given by Biot (1941, 1956, 1962a,b). A Lagrangian formulation is used for the solid skeleton, for which the displacements are assumed to be small, whereas a Eulerian formulation is adopted for the pore ¯uid, so attention is focused on the ¯uid motion through a control volume in movement with solid skeleton displacement. On the hypothesis of body forces being null, the equations can be written as follows: : : 11 u i ‡ 12 Vi ‡ b…ui

Vi †

ji; j ˆ 0

…1†

: 12 u i ‡ 22 Vi ‡ b…Vi

: ui †

;i ˆ 0

…2†

in which  ij,  are the partial stress and pressure in the solid and ¯uid phases per unit area of bulk (traction is positive); : : ui ; ui ; u i are the absolute displacement, velocity and acceleration of the solid skeleton; : Vi, Vi are the absolute velocity and acceleration of the pore ¯uid; n is the porosity; 11, 12, 22 are the following densities 11=(1-n) s+a, 12= a and 22 ˆ nf ‡ a , where s , f are the densities of the solid grains and pore ¯uid and a is the added mass (Biot, 1956; Gajo, 1996);

316

A. Gajo / International Journal of Plasticity 18 (2002) 313±344

b is related to the permeability kD (with the units of L2) and to the dynamic viscosity D (with the units of M/(LT)) of the pore ¯uid, through b=n2D/kD. The partial pressure  and the stress ij can be expressed in terms of the usual e€ective stress ij00 and pore pressure p through the relationships ij ˆ ij00 ‡ …

n†pij ;

 ˆ np

…3†

where ij is the Kronecker delta and  is a parameter depending on the bulk moduli of the solid skeleton and solid grains (Nur and Byerlee, 1971), which will be de®ned in Section 2.3. In Eqs. (1) and (2) the contribution of convective terms of the pore ¯uid in the transfer of momentum between the constituents has been disregarded, as in Biot (1956) and all the subsequent literature, because this contribution is negligible in view of the inaccuracy with which the permeability kD is evaluated (Zienkiewicz and Shiomi, 1984). In general, pore ¯uid viscosity is temperature-dependent, but this dependence will be neglected here for the sake of simplicity. Moreover, viscous coupling becomes frequency-dependent (Biot, 1956; Van Der Grinten et al., 1987; Gajo, 1995) for very high frequencies, but this e€ect will also be disregarded here. 2.2. Energy balance equations The energy balance equations for small temperature changes were stated by Pecker and Dereziewicz (1973) and by Bowen and Chen (1975). It is assumed that the solid and ¯uid temperatures are not locally equilibrated. In order to take energy convection due to ¯uid ¯ow into account, the ¯uid enthalpy gradient was introduced in the balance equations according to McTigue (1986), obtaining : T0 s ‡ k…Ts

Tf †

…1

: T0 f ‡ k…Tf

Ts †

nkf Tf;ii ‡ nf …Vi

n†ks Ts;ii ˆ 0

…4† : ui †hf;i ˆ 0

…5†

where: s, f are the entropies per unit volume of the solid and ¯uid phases; Ts, Tf are the absolute temperatures of the solid and ¯uid phases; T0 is the global reference temperature; ks, kf are the thermal conductivities of the solid grains and pore ¯uid; k is the coecient of interphase heat transfer; hf is the ¯uid enthalpy per unit of mass. In Eqs. (1) and (2) and (4) and (5), the phenomena analogous to the Dufour and Soret e€ects in solutions have been disregarded (Bowen and Chen, 1975; McTigue, 1986), but they could easily be introduced. As in Eqs. (1) and (2), the superposed dot in Eqs. (4) and (5) and in the ones that follow denotes the partial derivative with respect to time …@… †=@t†.

A. Gajo / International Journal of Plasticity 18 (2002) 313±344

317

According to Pecker and Dereziewicz (1973), the entropy per unit of volume of the two phases is assumed to be a linear function of s ˆ ss ii ‡ fs  ‡ …1

n†s Cps

f ˆ sf ii ‡ ff  ‡ nf Cpf

Ts T0

…6†

Tf T0

…7†

where ss and € are the isobaric coecients of linear and volumetric thermal expansion of the two phases [see also Eq. (11)], sf and fs are thermo-elastic coupling terms and Cps and Cpf are the speci®c heats of the two phases at constant pressure. The ¯uid enthalpy per unit of mass hf of the pore ¯uid can be expressed by hf ˆ Cpf Tf ‡ …T0 ff

1†

p : f

…8†

Eq. (8): can be derived from the ®rst law of thermodynamics which leads to : : T0 ff ˆ hf ‡ p=f (where p is assumed to be positive in traction and T  T0) and from the expression of ¯uid
entropy per unit of mass which, when linearised, leads : : : to ff ˆ …ff =f †p ‡ Cpf =T0 Tf . 2.3. Constitutive equations and the reference temperature for pore ¯uid For a linear elastic isotropic porous medium, the constitutive equations in terms of the e€ective stress ij00 and pore pressure p are as follows
…9† ij00 ˆ 2ij ‡ leij a011 …Ts T0 †ij a012 Tf T0f ij : p ˆ …

: : n†Qe ‡ nQE

: a021 Ts

T0



: a022 Tf

:  T0f

…10†

where
ij ˆ ui;: j ‡ uj;i =2 is the small strain tensor; e=ii, E=Vi,i are the solid skeleton dilation and the rate of pore ¯uid dilation; l,  are the isothermal LameÁ constant for a linear elastic isotropic solid skeleton; T0f is the reference temperature for the pore ¯uid;  is a parameter depending on the isothermal bulk modulus of the solid skeleton K ˆ …3l ‡ 2†=3 and solid grains (Ks) through  ˆ 1 K=Ks ; Q is a parameter depending on the isothermal bulk moduli of the solid skeleton (K), solid grains (Ks) and pore ¯uid (Kf) through 1=Q ˆ n=Kf ‡ … n†=Ks ; a011 , a012 , a021 , a022 are parameters related to the elastic moduli of the constituents and to the coecients of thermal expansion through the following relationships, derived from Pecker and Dereziewicz (1973): a021

ˆ 3ss …

a011 ˆ 3ss K; n†Q ‡ fs nQ;

a012 ˆ 3sf K a022 ˆ 3sf …

n†Q ‡ ff nQ

…11†

318

A. Gajo / International Journal of Plasticity 18 (2002) 313±344

The reference temperature for the pore ¯uid T0f is introduced in Eqs. (9) and (10) to take the heat convection due to ¯uid ¯ow properly into account. In fact, the state of stress ij00 and p is a€ected only by the variation in ¯uid temperature occurring inside the volume element, whereas the variation in ¯uid temperature occurring outside the volume element is unin¯uential. To clarify this concept, let us consider an ideal experiment: the injection of a hot ¯uid with Tf=310K into a volume element of a saturated porous medium with an initial temperature of Ts=Tf=300K (i.e. T0=300K); let us then approximate this injection to a sudden replacement of the cold pore ¯uid with the hot one. After the pore ¯uid replacement, due to the heat transfer between the hot ¯uid and the cold solid phase, the temperature of the hot ¯uid will decrease while the temperature of the solid phase will increase, leading respectively to a volumetric contraction of the ¯uid and a dilation of the solid phase. It is fairly obvious that this phenomenon can be appropriately taken into account by considering T0f=310K as the reference temperature for the pore ¯uid, and T0=300K as the reference temperature for the solid skeleton. Conversely, considering the reference temperature for the pore ¯uid as T0f=T0=300 K would imply an erroneous dilation of pore ¯uid inside the volume element (as if the ¯uid were heated inside the volume element, not outside, before being injected, as we have assumed here): the result would be an erroneous evaluation of the state of stress ij00 and p. That is why a di€erent pore ¯uid reference temperature has been introduced. Section 5 will examine the e€ects induced by T0f which, to the best of the author's knowledge, are usually disregarded in the literature. Obviously, when heat convection is null, T0f=T0. The need for a di€erent reference temperature for the pore ¯uid when heat convection is considered stems from the nature of the coupling between the two constituents: considering, for instance, the convective transport of a contaminant that does not interact mechanically with the solid skeleton (i.e. that does not a€ect the constitutive law), the contaminant concentration (which would correspond in this case to Tf would not require the introduction of anything analogous to the reference T0f. Since, as previously explained, T0f is related to the temperature of the pore ¯uid entered in the volume element, then, disregarding the in®nitesimal terms of higher order, the rate of change of T0f can be established from a simple volumetric balance of the temperatures of the ¯uid leaving and entering the volume element, as shown in Fig. 1, thus obtaining : T0f ˆ

…Vi

: ui †Tf;i

…12†

It is worth observing that both terms of Eq. (12) (multiplied by nfCpf) appear in the pore ¯uid energy balance [Eq. (5)]. 2.4. Resulting ®eld equations The following primary variables are considered ui s ˆ …Ts

T0 †

p f ˆ …Tf

Vi T0f † 0f ˆ …T0f

T0 †

A. Gajo / International Journal of Plasticity 18 (2002) 313±344

319

It is worth emphasising that f represents the variation in ¯uid temperature occurring inside the volume element, which a€ects the state of stress, whereas 0f represents the variation in ¯uid temperature due to the temperature di€erence of the ¯uid leaving and entering the volume element, which does not a€ect the state of stress. The resulting ®eld equations are . the momentum balance equation for the solid skeleton [from Eqs. (1), (3) and (9)]

: : 11 u i ‡ 12 Vi ‡ b…ui

Vi †

2ji; j

le;i

…

n†p;i ‡ a011 s;i ‡ a012 f;i ˆ 0

…13†

. the constitutive equation for the pore ¯uid [Eq. (10)]

: p

…

: n†Qe

: : : nQE ‡ a021  s ‡ a022  f ˆ 0

…14†

. the momentum balance equation for the pore ¯uid [from Eqs. (2) and (3)]

: 12 u i ‡ 22 Vi ‡ b…Vi

: ui †

np;i ˆ 0

…15†

. the energy balance equation for the solid skeleton [from Eqs. (3), (4), (6), (9) and (10)]

: : c11  s ‡ c12  f ‡ k…s

f

: : 0f † ‡ a11 T0 e ‡ a21 T0 E

…1

n†ks s;ii ˆ 0

…16†

Fig. 1. Schematic illustration of the volumetric balance of the temperatures of ¯uid leaving and entering the volume element for the evaluation of the rate of change of T0f.

320

A. Gajo / International Journal of Plasticity 18 (2002) 313±344

. the energy balance equation for the pore ¯uid [from Eqs. (3), (5), (7), (8), (9) and (10) minus Eq. (12) multiplied by nfCpf]

: : : : c21  s ‡ c22  f ‡ k…f ‡ 0f s † ‡ a12 T0 e ‡ a22 T0 E : ‡n…ff T0 1†…Vi ui †p;i ˆ 0


nkf f;ii ‡ 0f;ii ‡

…17†

. the rate of change of the pore ¯uid's reference temperature [Eq. (12)]

:  0f ‡ …Vi


: ui † f;i ‡ 0f;i ˆ 0

…18†

where cij are the speci®c heat at constant volume, given by a generalisation of the corresponding relationships for a single-phase medium (Pecker and Deresiewicz, 1973): c11 ˆ …1 c21 ˆ

n†s Cps T0 ‰3ss a11 ‡ fs a21 Š; c12 ˆ T0 ‰3ss a12 ‡ fs a22 Š T0 ‰3sf a11 ‡ ff a21 Š; c22 ˆ nf Cpf T0 ‰3sf a12 ‡ ff a22 Š

…19†

with a11 ˆ a011 ‡ … n†a021 ; a21 ˆ na021 ;

a12 ˆ a012 ‡ … a22 ˆ na022

n†a022

…20†

Obviously, the sum of Eqs. (17) and (18) multiplied by nfCpf gives the same equation as would be obtained from the pore ¯uid energy balance equation [Eq. (5)], after substituting Eqs. (3), (7), (8), (9) and (10). Finally, it is worth noting that for most of the low-frequency problems involving heat convection due to ¯uid ¯ow, the last term of Eq. (17), which is related to the pore pressure gradient p,i is usually negligible (Shao, 1997) and can be disregarded. Conversely, in the high-frequency regime generally involved in wave propagation, the contribution of the pore pressure gradient p,i may not be negligible. 3. Solution of the linearised problem When pore ¯uid displacements are small and consequently energy convection due to ¯uid ¯ow is negligible, the ®eld equations presented in the previous section are brought down to four linear equations (for a one-dimensional problem) and a Lagrangian formulation may be used also for the pore ¯uid (i.e. the absolute : displacement of the pore ¯uid U may be taken as the unknown ®eld, thus V ˆ U i i and i : Vi ˆ U i ). As a result, for a one-dimensional problem along the x-direction the four linear equations are  : : d11 u;xx d12 U;xx ‡ a11 s;x ‡ a12 f;x ˆ 0 …21† 11 u ‡ 12 U ‡ b u U

A. Gajo / International Journal of Plasticity 18 (2002) 313±344

: 21 u ‡ 22 U ‡ b U

: u



: : c11  s ‡ c12  f ‡ k…s

: : f † ‡ a11 T0 u;x ‡ a21 T0 U;x

…1

: : c21  s ‡ c22  f ‡ k…f

: : s † ‡ a12 T0 u;x ‡ a22 T0 U;x

nkf f;xx ˆ 0

d21 u;xx

d22 U;xx ‡ a21 s;x ‡ a22 f;x ˆ 0 n†ks s;xx ˆ 0

321

…22† …23† …24†

where the index notation was dropped for the sake of simplicity and dijs are d11 ˆ l ‡ 2 ‡ …a n†2 Q; d12 ˆ n… d21 ˆ n… n†Q; d22 ˆ n2 Q

n†Q

…25†

This system of equations can be solved in time domain by means of the transfer functions and Fourier transforms method. With this approach, all relevant quantities are assumed to depend on the harmonic function exp‰i!…t x=v†Š, where v is the phase velocity. The dispersion relation can be obtained from the null determinant condition of the matrix of the coecients deduced from Eqs. (21)±(24). 2 3 s11 ib=! s12 ‡ ib=! ia11 =…!v† ia12 =…!v† 6 s21 ‡ ib=! s22 ib=! ia21 =…!v† ia22 =…!v† 7 7ˆ0 det6 …26† 4 a11 T0 =v 5 a21 T0 =v t11 ‡ …1 n†ks =v2 t12 2 a12 T0 =v a22 T0 =v t21 t22 ‡ nkf =v where sij ˆ ij dij =v2 and tij ˆ icij =! ‡ k=!2 . It can be demonstrated that the resulting four longitudinal waves are dispersive and dissipative, as discussed by Pecker and Dereziewictz (1973). For each longitudinal wave and for each frequency component !, the initial amplitudes of ¯uid displacements [U0j (!)] and solid and ¯uid temperatures [s0j(!) and f0j (!)] can be expressed in terms of the respective initial amplitudes of solid displacements [u0j (!)], using the ratios ji (!) de®ned as follows: j1 …!† ˆ U0j …!†=u0j …!†;

j2 …!† ˆ s0j …!†=u0j …!†;

j3 …!† ˆ f0j …!†=u0j …!†

…27†

where the index j varies between 1 and 4 and refers to each longitudinal wave. The ratios ji (!) can be evaluated by solving the following linear system, obtained from just three equations (the ®rst three in this case) of the initial system [Eqs. (21)±(23)]. 2

3 2 3 12 d12 =v2j ‡ ib=! ia11 = !vj ia12 = !vj j1 …!†

7 6 2 ia21 = !vj ia22 = !vj 5 4 j2 …!† 5 4 22 d22 =vj ib=! 2 2 j3 …!† ic11 =! ‡ k=! ‡ …1 n†ks =vj ic12 =! k=!2 a21 T0 =vj 2 3 2 11 ‡ d11 =vj ‡ ib=! ˆ 4  ‡ d =v2 ib=! 5 21

21

j

a11 T0 =vj

…28†

322

A. Gajo / International Journal of Plasticity 18 (2002) 313±344

This system must be solved for each longitudinal wave (vj). Finally, the amplitudes of the longitudinal waves in terms of solid displacements (u0j) can be evaluated from the boundary condition. For instance, if we consider a boundary condition consisting of an applied ¯uid pressure p0 (!) with constant solid and ¯uid temperatures, the amplitudes of the longitudinal waves in terms of absolute solid displacements u0j (!) can be obtained by solving the following linear system: 4  X jˆ1 4  X jˆ1 4  X

! id11 vj

! id12 j1 …!† vj

! id21 vi

! id22 j1 …!† vj

a11 j2 …!†

 a12 j3 …!† u0j …!† ˆ …1

a21 j2 …!†

 a22 j3 …!† u0j …!† ˆ np0 …!†

n†p0 …!†

…29†

…30†

 j2 u0j …!† ˆ 0

…31†

 j3 u0j …!† ˆ 0

…32†

jˆ1 4  X jˆ1

The transfer functions for the solid and ¯uid temperature variations are consequently s …x; !† ˆ

4  X j2 …!†u0j …!†exp

i!x=vj




…33†

jˆ1

f …x; !† ˆ

4  X

j3 …!†u0j …!†exp

i!x=vj




…34†

jˆ1

Figs. 2 and 3 show the solid and ¯uid temperature variations obtained with the proposed method and corresponding to the transit of the ®rst two longitudinal waves generated by an applied ¯uid pressure pulse having the shape of a single step and an amplitude of 100 kPa. Low, medium and high viscous coupling values were considered (kD=2.010 9, kD=2.010 12 and kD=2.010 18 m2). The material properties are given in Section 5 and correspond to a sand at a low con®ning pressure saturated with oil. The e€ects of dispersion and damping on wave propagation induce a change in the shape of the single step pulse, which is consistent with the one evaluated analytically under isothermal conditions (Gajo, 1995; Gajo and MongiovõÁ, 1995).

A. Gajo / International Journal of Plasticity 18 (2002) 313±344

323

4. Finite element solution of the non-linear problem The di€erential solving system consisting of Eqs. (13)±(18) can be discretised spatially using the standard interpolation schemes (Hughes, 1987; Zienkiewicz and Taylor, 1989). u ˆ Nu u

p ˆ NT p

V ˆ Nu V

s ˆ NT s

f ˆ N T  f

0f ˆ NT 0f

where, in general, uT= [u1 , u2 , u3 ], VT=[V1, V2, V3]; u; p; V ; s ; f and 0f are the nodal solid displacements, pore pressure, ¯uid velocities and temperature values; Nu is the shape function for the displacement and velocity ®elds and NT is the shape function for the pore pressure and temperature ®elds.

Fig. 2. History of solid temperature at a depth of 1.00 cm, due to transit of the ®rst two longitudinal waves generated by a stepping variation in ¯uid pressure (results obtained using the transfer functions method).

324

A. Gajo / International Journal of Plasticity 18 (2002) 313±344

In the following paragraphs, we consider only a one-dimensional problem for the sake of simplicity. In this case, the one-dimensional discretisation was obtained with three-node quadratic ®nite elements for u and V, and two-node linear ®nite elements for p, s, f and 0f. In this way, the stresses, strains and temperatures in Eqs. (13)± (17) are approximated by interpolating functions having the same order. The ®rst ®ve di€erential equations [Eqs. (13)±(17)] are semidiscretised using Galerkin weighting functions. For the sixth equation [Eq. (18)], due to the presence of heat convection, a weighting function W of the discontinuous Petrov-Galerkin type is used (Brooks and Hughes, 1982). The components of W are connected to those of the linear shape function NT through Wi ˆ NTi ‡

h @NTi ‰sign…V 2 @x

: u†Š

…35†

Fig. 3. History of ¯uid temperature at a depth of 1.00 cm, due to transit of the ®rst two longitudinal waves generated by a stepping variation in ¯uid pressure (results obtained using the transfer functions method).

A. Gajo / International Journal of Plasticity 18 (2002) 313±344

325

where x is the spatial coordinate, the index i refers to the two nodes of the linear ®nite element, h is the length of the linear ®nite element and is a factor which should depend on the Peclet parameter, but it proves here to be =1, since there are no di€usive terms in Eq. (18). The resulting non-symmetrical system is 3 2 2 : 3 3 u 3 u 2 2 B 0 M12 0 0 0 6 : 7 M11 0 0 0 0 0 6 7 6 p 7 6 T 6 p 7 7 6 0 7 6 G1 P 0 A21 A22 0 7 0 0 0 0 0 76 76 6 : 7 7 7 6 76 7  7 6 6 T T 7 6  V V 7 6 6 M12 0 0 0 0 0 76 B 0 M 0 0 0 22 7 6 6 76 : 7 76 6 6 ‡ 7 7  6 A 6 0 0 0 C11 C12 0 7 0 0 0 0 07  s 7  s 7 76 76 6 11 6 7 6 6 7 : 7 7 6 7 6  7 4 0 0 0 C21 C21 0 56 0 0 0 0 0 56 6  f 7 4 A21 6  f 7 5 4 4 : 5 0 0 0 0 0 C22 0 0 0 0 0 0     0f  0f 32 3 2 3 2 K G1 fs B A11 A12 0 u 76 p 7 6 0 7 60 0 GT2 0 0 0 76 7 6 7 6 76  7 6 7 6 76 V 7 6 ff 7 6 0 G2 B 0 0 0 76 7ˆ6 7 ‡6 …36† 7 6 7 60 6 0 A12 L11 ‡ E E E 7 76  s 7 6 qs 7 6 76 7 6 7 6 40 D2 A22 ET L22 ‡ E L22 ‡ E 54  f 5 4 qf 5 0 0 0 0 D1 D1 0  0f The detailed expression of these matrices is given in the Appendix. Time integration of the di€erential system is resolved in a coupled way using Newmark's (1959) algorithm, in the form given by Bathe (1982). This algorithm is modi®ed to admit a speci®ed velocity boundary condition (Gajo et al., 1994). In all the examples presented in the next sections, the parameters of the Newmark scheme will be such as to ensure unconditional stability (which is related to the implicit nature of the computations and to the choice of the algorithm's coecients 50.5 and 50.25 (+0.5)2, when damping is neglected) and give a slight numerical damping in order to avoid spurious high-frequency oscillations in the numerical results (=0.70 and =0.36). The resulting system is obviously non-linear due to the presence of the terms D1 and D2, which take heat convection into account and which depend on the unknown values of the relative velocities of the two components … h : i
 …37† D1 ˆ WT nf Cpf Nu V u BT d



… D2 ˆ



h NTT n…1

 ff T0 †Nu V

: i u BT d

…38†

It is worth noting that the introduction of the new reference temperature for the pore ¯uid 0f has the advantage of concentrating the use of the discontinuous Petrov±Galerkin weighting functions in the last equation alone, thus involving very few terms. Thus, all the other equations are spatially discretised according to the standard procedures.

326

A. Gajo / International Journal of Plasticity 18 (2002) 313±344

The solution of the non-linear system is obtained by the direct iterative method, :
i.e. carrying out successive substitutions of the calculated relative velocities V u until the consecutive values are consistent to within a speci®ed tolerance. 5. Validation of the numerical procedures Two numerical examples are developed to validate the numerical scheme for static problems dominated by heat convection and for dynamic problems in which energy convection is disregarded. The numerical results concern the transient one-dimensional problems of hot ¯uid injection in a steady seepage ¯ow and non-isothermal wave propagation analysed using the linearised theory. The ®rst problem concerns the injection of hot ¯uid and is solved to validate the proposed theoretical approach based on a di€erent pore ¯uid reference temperature T0f and to ascertain the reliability of the numerical scheme in a practical engineering problem where heat convection is dominant. For this purpose, a high permeability value was considered, i.e. kD=2.010 9 m2. The second problem concerns wave propagation and is solved to study the accuracy of the ®nite element procedures in representing a dispersive, non-isothermal wave propagation in a saturated soil in which energy convection is disregarded. This is done by comparing the numerical results with those obtained using the transfer functions and Fourier transforms method presented in Section 3. Three permeability values were considered, i.e. kD=2.010 9, kD=2.010 12 and kD=2.010 18 m2, corresponding to low, medium and high viscous coupling, respectively. This test problem con®rms the reliability of the ®nite element method in solving practical problems with any range of viscous coupling and frequency content. In all the following examples, the material properties were deduced from the measurement by Pecker and Deresiewicz (1973), considering a small solid skeleton sti€ness value: these properties approximately correspond to a loose sand at a low e€ective con®ning pressure (about 100 kPa), saturated with oil. A solid skeleton with such a relatively low sti€ness was considered in order to enhance the thermal coupling e€ects in wave propagation. l ˆ 115:4 MPa Ks ˆ 36:0 GPa f ˆ 1000 kg=m3 ss ˆ 0:9  10 5 …K† 1 ff ˆ 1:152  10 3 …K† 1 Cpf ˆ 2093 J=…kg  K†

 ˆ 76:9 MPa Kf ˆ 1:9 GPa s ˆ 2650 kg=m3 sf ˆ 0:108  10 5 …K† 1
k ˆ 4:309  105 J= m3  s  K ks ˆ 3:057 J=…m  s  K†

n ˆ 0:40 D ˆ 2:0  10 3 Pa  s a ˆ 200 kg=m3 fs ˆ 2:88  10 4 …K† 1 Cps ˆ 1172 J=…kg  K† kf ˆ 0:149 J=…m  s  K†

These properties are such that in a pervious porous medium the longitudinal waves of the ®rst and second kind propagate at a velocity of 1547.0 and 264.4 m/s under non-isothermal conditions, and at 1363.8 and 262.8 m/s under isothermal conditions, i.e. the thermal coupling leads to an increase of about 20% in the velocity of the longitudinal wave of the ®rst kind and leaves the velocity of the longitudinal wave of the second kind almost unchanged.

A. Gajo / International Journal of Plasticity 18 (2002) 313±344

327

Since spatial and temporal re®nement are linked in some way, a spatial re®nement needs to be accompanied by a suitable temporal re®nement (Hiremath et al., 1988). In all the following examples, therefore, we have chosen time steps suitably smaller than the time interval needed for the faster wave to propagate between two adjacent nodes; the time steps were also small enough to ensure the convergence of the direct substitution method within a few iterations for the non-linear analysis. Due to the essential theoretical nature of the present study, at each iteration the linear system of discretised equations was solved simultaneously and no staggered or partitioned solution procedure was employed. Consistent mass sub-matrices were used. 5.1. Test problem 1: hot ¯uid injection We present the solution of the problem of hot ¯uid injection in a ®nite soil column 2.40 cm in length, shown in Fig. 4. In the ®rst 80 ms (i.e. from t=0 to t=80 ms), the ¯uid pressure at the top surface is increased from 0 to 0.25 kPa with respect to the (null) ¯uid pressure at the bottom surface. At the end of that period, the steady state condition in the seepage ¯ow is

Fig. 4. Representative soil column.

328

A. Gajo / International Journal of Plasticity 18 (2002) 313±344

brie¯y established due to the high permeability and the ®nal relative velocity proves : equal to …V u†  2.60 cm/s and is directed downwards. After the ®rst 150 ms, the external ¯uid temperature at the top surface increases steadily by 2 K in 150 ms (i.e. from t=150 to t=300 ms) with respect to the initial temperature of the saturated porous medium. Spatial discretisation is uniform and involves 20 one-dimensional quadratic ®nite elements for u and V, and 20 linear ®nite elements for p, s, f and 0f. The temporal integration involves 12 000 steps of 0.05 ms, for a total time of 0.6 s. Figs. 5 and 6 show the distribution of solid (s) and ¯uid (f and 0f) temperatures in t=0.55 sec and the history of the temperatures at a depth x=0.36 cm: the propagation of a hot ¯uid wave due to ¯uid motion and a simultaneous increase in solid temperature are observable. While the solid temperature (s) increases due to interphase heat transfer, the ¯uid temperature inside the solid skeleton (f) decreases and the reference temperature 0f consequently increases up to values greater than the increase in the external ¯uid temperature (2 K). This is consistent with the fact that, as long as the solid phase is cooler than the ¯uid one, there is a continuous delivery of thermal energy in the volume element through ¯uid ¯ow and thus 0f increases continuously. The ®nal value of f can be estimated by disregarding the work done by the pore pressure gradient p,i and by equating the thermal energy needed to heat the solid phase to the thermal energy lost by the ¯uid due to cooling, i.e. f ˆ

…1

n†Cps s nCpf

Fig. 5. Distribution of ¯uid and solid temperatures at time t=0.55 s, due to hot ¯uid injection (test problem 1).

A. Gajo / International Journal of Plasticity 18 (2002) 313±344

329

In the same assumptions, the ®nal values of s and f+0f will equate to the ®nal increase of the external ¯uid temperature at the top surface (i.e. s=f+0f=2 K). From the previous relationship, the ®nal value of 0f will thus be approximately   …1 n†Cps 0f ˆ 1 ‡ s : nCpf It is worth noting that, thanks to the use of a weighting function of the discontinuous Petrov±Galerkin type [Eq. (35)], the solution is very smooth and all spurious oscillations, which are known to occur when the usual Galerkin weighting functions are employed in a convection-dominated problem, are avoided.

Fig. 6. History of ¯uid and solid temperatures at a depth of 0.36 cm, due to hot ¯uid injection (test problem 1).

330

A. Gajo / International Journal of Plasticity 18 (2002) 313±344

Due to the small variation in the pore ¯uid temperature, the volumetric contraction of the ¯uid due to cooling can be disregarded (see Fig. 6b) and thus the expected transit time of the hot ¯uid wave can easily be evaluated from the velocity of the steady seepage ¯ow (t=363.2 ms): Fig. 6a illustrates the satisfactory comparison between the expected and computed hot thermal wave transit times. Let us consider what happens when the new pore ¯uid reference temperature T0f is disregarded (as is usually the case), assuming T0f=T0 in the constitutive equation. This is done by establishing that the damping and the sti€ness matrix of Eq. (36) are given respectively by 3 2 B 0 M12 0 0 0 6 GT1 P 0 A21 A22 A22 7 7 6 6 BT 0 M 0 0 0 7 22 7 6  6 A11 0 0 C11 C12 0 7 7 6  4 A 0 0 C21 C22 0 5 21 0 0 0 0 0 C22 …39† 2 3 K G1 B A11 A12 A12 T 60 7 0 G 0 0 0 2 6 7 6 0 G2 7 B 0 0 0 7 and6  60 0 A12 L11 ‡ E E E 7 6 7 40 D2 A22 ET L22 ‡ E L22 ‡ E 5 0 0 0 0 D1 D1 The results shown in Fig. 7 are thus obtained. The computed hot ¯uid wave transit time is clearly larger than was estimated: this is due to the erroneous reference

Fig. 7. The e€ects of T0f on the history of ¯uid temperature at a depth of 0.36 cm in the problem of hot ¯uid injection (test problem 1).

A. Gajo / International Journal of Plasticity 18 (2002) 313±344

331

temperature being assumed for the pore ¯uid. In fact, due to the wrong assumption that T0f=T0 in Eq. (36), the pore ¯uid dilates due to an erroneous temperature increase (Tf T0>0), instead of compressing due to the actual temperature decrease (Tf T0f<0). As a result, the hot thermal wave is slowed down by the wrong pore ¯uid dilation occurring during wave transit. This can be better appreciated in Fig. 8, which shows the history of the `change in ¯uid content' =n(e E) (obtained by numerically integrating the ®nite element results, which are in terms of rates). In the ®rst 80 ms, the ¯uid content decreases due to compression of the solid skeleton induced by seepage forces (see also Fig. 9); subsequently, if the new pore ¯uid

Fig. 8. The e€ects of T0f on the history of change in ¯uid content  at a depth of 0.36 cm in the problem of hot ¯uid injection (test problem 1).

Fig. 9. History of the vertical displacements of the solid skeleton evaluated at a depth of 0.36 cm in the problem of hot ¯uid injection (test problem 1).

332

A. Gajo / International Journal of Plasticity 18 (2002) 313±344

reference temperature T0f is taken into account as proposed here, the ¯uid content starts to increase due to thermal compression of the pore ¯uid and dilation of the solid skeleton; conversely, if the new pore ¯uid reference temperature T0f is disregarded (i.e. T0f=T0), the ¯uid content decreases due to thermal expansion resulting from the wrong reference temperature T0 being applied to the pore ¯uid. Fig. 9 shows the solid vertical displacements (u) evaluated at a depth of 0.36 cm: initially (from t=0 to t=80 ms) the soil layer is compressed due to the previouslymentioned seepage forces, then it expands due to the increase in solid skeleton temperature induced by the interphase heat transfer. For the sake of comparison, Fig. 9 also shows the vertical displacements of the solid skeleton, calculated on the assumption that the interphase heat transfer is null (k=0) and that the thermal conductivity of the ¯uid is negligible (kf  0). In this case, the two constituents are thermally uncoupled and the variations in the temperatures of the solid (s) and of the ¯uid inside the solid skeleton (f) are consequently due only to the adiabatic compression induced by seepage forces and are negligible; the resulting vertical solid displacements are practically null too, therefore, as shown in Fig. 9. It can be concluded that the proposed theoretical framework and numerical procedures are reliable for the analysis of problems dominated by heat convection. 5.2. Test problem 2: linearised analysis of non-isothermal wave propagation Its worth remembering that a validation of a similar mixed ®nite element formulation under isothermal conditions (u p U formulation) has been done elsewhere (Gajo et al., 1994) by comparing the numerical results with an analytical solution of Biot's equations (Gajo and MongiovõÁ, 1995). In order to validate the proposed numerical scheme for the linearised analysis of non-isothermal wave propagation, the energy convection was disregarded and the numerical results were compared with those obtained using the transfer functions method. Energy convection was disregarded by uncoupling the last equation from the others in the solving system (36): this is tantamount to assuming that the sti€ness matrix of Eq. (36) is as follows 2 3 K G1 B A11 A12 0 6 0 0 GT2 0 0 0 7 6 7 6 0 G2 B 7 0 0 0 6 7: …40† 6 0 0 A L11 ‡ E 7 E 0 12 6 7  T 40 0 A E L22 ‡ E 0 5 22 0 0 0 0 0 D1 In fact, the last equation of the solving system is already uncoupled in the mass and damping matrices. The following problem was considered: a pore pressure varying according to a stepping pulse was applied at the top free surface of a soil column 5.00 cm in length (Fig. 4). Spatial discretisation is uniform and involves 50 one-dimensional quadratic ®nite elements for u and V, and 50 linear ®nite elements for p, s, f and 0f. The temporal

A. Gajo / International Journal of Plasticity 18 (2002) 313±344

333

integration involves a uniform time step of 0.25 ms; this time step approximately equates to the time interval needed for the faster wave to propagate between two adjacent nodes. The calculated variations in solid and ¯uid temperatures (s and f) at a depth of 1.00 cm are shown in Figs. 10 and 11 and are compared with the results obtained using the transfer functions and Fourier transforms method for di€erent viscous coupling values. Figs. 10 and 11 show the temperature variations associated with the transit of longitudinal waves of the ®rst and second kind. Numerical results show a greater rise time than the analytical solution, but they could be improved by using a ®ner spatial and temporal discretisation. It is worth noting that the numerical results are able to capture the most important aspects of the results obtained using the

Fig. 10. History of solid temperature at a depth of 1.00 cm, due to wave propagation: comparison between the results obtained with the transfer functions method and the ®nite element method (test problem 2).

334

A. Gajo / International Journal of Plasticity 18 (2002) 313±344

transfer functions and Fourier transforms method, in any range of permeability and for both types of longitudinal wave. The ®nite element solution not only correctly reproduces the wave propagation forms for the two extreme permeability conditions (with two waves for kD=2.010 9 m2 and only 1 wave for kD=2.010 18 m2), it also correctly predicts behaviour for an intermediate permeability condition (kD=2.010 12 m2), in which dispersion is important. The overall response of the numerical procedures developed for non-isothermal dynamic conditions is therefore similar to the one presented by the isothermal ®nite element solutions for dynamic problems in saturated porous media (Gajo et al., 1994). From the good consistency between the ®nite element results and those obtained with the transfer functions method in the frequency range represented by the selected spatial and temporal discretisation, it can be concluded that the proposed

Fig. 11. History of ¯uid temperature at a depth of 1.00 cm, due to wave propagation: comparison between the results obtained with the transfer functions method and the ®nite element method (test problem 2).

A. Gajo / International Journal of Plasticity 18 (2002) 313±344

335

numerical scheme is a reliable tool for the analysis of dynamic non-isothermal problems in any range of viscous coupling and frequency content. 6. Non-linear analysis of non-isothermal wave propagation In the previous section, the ®nite element procedures were validated and proved to be a reliable tool for the analysis of dynamic non-isothermal problems in which energy convection is taken into account. It is then possible to analyse the e€ects induced by energy convection due to pore ¯uid ¯ow on a dispersive non-isothermal wave propagation in a saturated porous medium. The same problem as in Subsection 5.2 was analysed: a pore pressure varying according to a stepping pulse was applied at the top free surface of a soil column 5.00 cm in length (Fig. 4). Here again, three permeability values were considered, i.e. kD=2.010 9, kD=2.010 12 and kD=2.010 18 m2, corresponding to low, medium and high viscous coupling values, respectively. The same spatial and temporal re®nements as in Subsection 5.2 were employed. In particular, the selected temporal re®nement (0.25 ms) ensured the convergence of the direct substitution method within a few iterations. The propagation of waves generated by a maximum applied ¯uid pressure of 0.10 and 10.00 MPa was considered. Fig. 12 compares the results of the linearised theory (without energy convection) and non-linear analysis (with energy convection) obtained with the ®nite element method [energy convection was disregarded in the linearised analysis with the method explained in Subsection 5.2, Eq. (40)]. The results refer to the high permeability case (kD=2.010 9 m2), for which energy convection is more important. The results practically coincide even for the large wave generated by an applied ¯uid pressure pulse of 10 MPa; only the amplitude of the longitudinal wave of the second kind is slightly a€ected by the non-linear analysis. The e€ects induced by heat convection on non-isothermal wave propagation can thus be considered negligible, even for a wide-amplitude input pulse. Fig. 13 shows the history of 0f at a depth of 1.00 cm for two applied ¯uid pressure input pulses with amplitudes of p=10 MPa and p=0.10 MPa: the amplitude of 0f clearly varies with the square of the amplitude of the pulse applied. Let us analyse in Fig. 14 the variation in the new reference temperature 0f: in the case of high permeability (kD=2.010 9 m2), 0f is about ®ve orders of magnitude smaller than f for an applied ¯uid pressure pulse of 100 kPa (0f would be three orders of magnitude smaller than f for a pulse amplitude of 10 MPa). Such a large di€erence between 0f and f explains why the e€ects of heat convection are so negligible. This occurs despite the fact that the relative velocities between the solid and : ¯uid phases are as large as those of the steady seepage ¯ow …V u  2:30 cm=s†, as shown in Fig. 15. The fact that the ®nal value of 0f is small despite the large relative velocity between the two phases can be explained by the small spatial extension of the ¯uid temperature gradient (corresponding to the length of the stepping pulse in the ®nite element mesh) and by the small relative displacements between the two phases. From said dependence of 0f on the square of the applied pulse, it can be deduced that

336

A. Gajo / International Journal of Plasticity 18 (2002) 313±344

heat convection is expected to be relevant for a pressure pulse with a very large amplitude (larger than about 1000 MPa for the material properties considered here). It is worth noting that the computed 0f is consistent with a rough estimate that can be obtained from the linearised analysis. In fact, from the negligible contribution of heat convection observed (0f  0), the gradient f,i and the relative velocity : …V u† appearing in Eq. (18) can be estimated from the linearised analysis. Due to the higher frequency cut-o€ induced by spatial and temporal discretisation, the propagating pulse can be approximated to a pressure step with a limited rise time r (r  4 ms from Fig. 14). If the permeability is high and the distance travelled is small, as in this case (kD=2.010 9 m2), the pulse shape is not distorted :
(Gajo, 1995). Let us use ^ f to indicate the variation in ¯uid temperature and V^ u^ for the relative velocity after: the transit of the ®rst wave: from the linearised analysis ^ f =15.4610 3 K and V^ u^ =0.023 m/s. At a given position, during the rise time the gradient f,i is constant and equates to f;i ˆ ^ f =…v1 r† (where v1=1547 m/s is the velocity of the ®rst wave) and the relative velocity increases linearly with respect to time t according to :
: V u ˆ t V^ u^ =r. As a result, the value of ^ 0f at a given position after the transit of the ®rst wave can be obtained by integrating Eq. (18) with respect to time  : …r ^ f V^ u^  : ^ 0f ˆ f;i …V u†dt  ˆ 0:00011  10 3 K 2v1 0 which is consistent with the computed value (^ 0f =0.0001010

3

K from Fig. 14).

Fig. 12. History of ¯uid temperature at a depth of 1.00 cm, due to wave propagation: comparison between the results obtained with the linearised theory (i.e. without energy convection) and with the nonlinear analysis (i.e. with energy convection) for an applied ¯uid pressure pulse of p=10 MPa.

A. Gajo / International Journal of Plasticity 18 (2002) 313±344

337

Fig. 14 shows that for the intermediate permeability value (kD=2.010 12 m2), the ®nal value of 0f, and thus also the contribution of heat convection, is smaller than the one obtained in the case of high permeability (kD=2.010 9 m2). This is due to the lower relative velocity between the two phases (Fig. 15): in particular, before the arrival of the second wave, the relative velocity is completely nulli®ed by the high viscous coupling. Finally, for a low permeability (kD=2.010 18 m2), heat convection gives no contribution in energy ¯ux because the high viscous coupling completely nulli®es the relative velocity between the two phases. It is worth comparing the energy ¯uxes due to heat convection and conduction. The energy ¯ux due to heat conduction can be computed from the solution based on the transfer function method and can be compared with the energy ¯ux due to heat convection: this comparison was drawn in terms of 0f (which is directly connected to heat convection) and of the fraction of ¯uid temperature variation fx , which is due to heat conduction. The result is that, for an applied ¯uid pressure of 100 kPa and for a frequency content of about 25 kHz, f is about four orders of magnitude smaller than 0f. Said di€erence decreases as the frequency content increases. In other words, during wave propagation in saturated pervious porous media, the energy ¯ux due to heat convection is not negligible by comparison with the ¯ux due to heat conduction (as implicitly assumed in the linearised theory), it is much greater. However, since both heat convection and conduction are negligible, wave

Fig. 13. History of ¯uid reference temperature 0f at a depth of 1.00 cm in the framework of non-linear analysis: comparison between the results obtained with a small (p=100 kPa) and large (p=10 MPa) applied ¯uid pressure input pulse.

338

A. Gajo / International Journal of Plasticity 18 (2002) 313±344

propagation in a saturated pervious porous medium can be considered completely adiabatic, even when heat convection is considered and a large ¯uid pressure pulse is applied. Thus the reliability of the linearised analyses is demonstrated. Let us ®nally draw some practical conclusions which could be useful for the interpretation of dynamic tests. From the previous theoretical analysis, dynamic measurements can be reliably interpreted with the linearised theory, which implies an adiabatic wave propagation however. As a consequence, Biot's equations with adiabatic elastic parameters can be employed. Such adiabatic values can be obtained from a generalisation of the relationships for a single-phase medium. From the constitutive equations (Eqs. (9) and (10)] and from the energy balance equations [Eqs. (16)±(18)] with a null heat ¯ux, the following adiabatic constitutive relationships, which are relevant for a longitudinal wave (i.e. one-dimensional compression with a null lateral strain), can be obtained  :  (     1  ) :  d11 d12 11 21 c11 c12 11 12 u;x  ‡ T0 …41† : ˆ d21 d22 12 22 c21 c22 21 22 V;x  If water is the pore ¯uid and the solid skeleton consists of quartz, Eq. (41) can be approximated to satisfaction with the coecients dij computed from Eq. (25) using the adiabatic parameters of the solid skeleton and pore ¯uid …Ks †adiabatic  Ks ;

…Kf †adiabatic ˆ Kf ‡

T0 2  Kf f Cpf ff

…42†

Fig. 14. History of ¯uid reference temperature 0f at a depth of 1.00 cm in the framework of non-linear analysis, obtained with a small (p=100 kPa) applied ¯uid pressure input pulse.

A. Gajo / International Journal of Plasticity 18 (2002) 313±344

339

Fig. 15. Relative velocities between the ¯uid and solid phases at a depth of 1.00 cm in the framework of non-linear analysis, obtained with a small (p=100 kPa) applied ¯uid pressure input pulse.

as done by Gajo et al. (1997) [the adiabatic and isothermal sti€ness of water is (Kf)adiabatic=2.213 GPa and Kf=2.183 GPa, respectively]. Although the di€erence between adiabatic and isothermal properties is fairly small, it may nonetheless be important in the evaluation of dynamic tests (Gajo and MongiovõÁ, 1994). 7. Conclusions The momentum and energy balance equations for the non-isothermal dynamic behaviour of saturated porous media are reiterated, emphasising the terms related to energy convection due to ¯uid ¯ow. In particular, due to heat convection, the need for a di€erent pore ¯uid reference temperature with respect to the one for the solid phase is discussed and taken into account: this stems from the nature of the thermal coupling between the two phases and is usually disregarded in the literature, leading, for instance, to an erroneous evaluation of the velocity of propagation of the thermal wave in the problem of hot (or cold) ¯uid injection in a steady seepage ¯ow. A novel solution for linearised non-isothermal wave propagation is presented. This solution is based on the transfer functions and Fourier transforms method and the results are used to validate the numerical procedures developed to solve the nonlinear problem, for a wide range of viscous couplings. These procedures are based on the ®nite element method and spatial discretisation is achieved with the usual

340

A. Gajo / International Journal of Plasticity 18 (2002) 313±344

interpolation schemes and weighting functions of the discontinuous Petrov±Galerkin type. Six primary variables are considered. Finally the non-linear analysis of non-isothermal wave propagation is performed under a wide range of viscous couplings using the ®nite element method. The outcome is a great relative velocity between the ¯uid and solid phase, so the energy ¯ux due to heat convection is much greater than the one due to heat conduction, contrary to the implicit assumption of the linearised theory. Even when the permeability is great and the applied input pulse is large, the wave propagation can still be considered completely adiabatic, however, despite the heat ¯ux due to energy convection. The validity of the linearised analyses of non-isothermal wave propagation can therefore be considered as demonstrated and the previously-mentioned potential inconsistencies between the theoretical and the experimental results are probably related to other factors. Acknowledgements The author wishes to thank Professor B. Loret for his precious comments and support. Appendix A The following are the detailed expressions of the sub-matrices of Section 4, for a one-dimensional problem. The extension of these sub-matrices to two- or threedimensional problems can easily be achieved for all the terms (see, for instance, Gajo et al., 1994), except for the convective term (sub-matrices Di) for which reference should be made, for example, to Brooks and Hughes (1982). The symbol t^ represents the total boundary traction, whereas q^ s and q^ f are the boundary heat ¯uxes applied to the solid skeleton and pore ¯uid. … M11 ˆ NTu 11 Nu d



… M12 ˆ



… M22 ˆ



… Bˆ



NTu 22 Nu d

NTu bNu d

… C11 ˆ

NTu 12 Nu d



NTu c11 NT d

A. Gajo / International Journal of Plasticity 18 (2002) 313±344

… C12 ˆ



… C21 ˆ



… C22 ˆ



…

C22 ˆ



…

A11 ˆ



…

A12 ˆ



…

A21 ˆ



…

A22 ˆ



… Kˆ





… G2 ˆ



… Pˆ



NTT c22 NT d


WT nf Cpf NT d

NTT …a11 T0 †Bu d

NTT …a21 T0 †Bu d

NTT …a12 T0 †Bu d

NTT …a22 T0 †Bu d

BTu …

NTT

A11 ˆ

n†NT d

BTu nNT d

…

… A12 ˆ

NTT c21 NT d

BTu …l ‡ 2†Bu d

… G1 ˆ

NTT c12 NT d



1 NT d

Q

NTu a011 BT d

NTu a012 BT d

341

342

A. Gajo / International Journal of Plasticity 18 (2002) 313±344

… A21 ˆ



… A22 ˆ



… Eˆ



…

…

… D1 ˆ



… D2 ˆ

 0  a21 NT d

Q

NTT

 0  a22 NT d

Q

NTT kNT d

L11 ˆ L22 ˆ

NTT



BTT …1

n†ks BT d

BTT …nkf †BT d

h
 WT nf Cpf Nu V h NTT n…1

: i u BT d

 ff T0 †Nu V

: i u BT d

The nodal forces are given by … fs ˆ NTu …1 n†t^d … ff ˆ … qs ˆ … qf ˆ

NTu …n†t^d NTT q^ s d NTT q^ f d

In the previous relationships Bu ˆ LT Nu B T ˆ r T NT with LT= rT ˆ ‰@=@xŠ for a one-dimensional problem.

A. Gajo / International Journal of Plasticity 18 (2002) 313±344

343

References Bathe, K.J., 1982. Finite Element Procedures in Engineering Analysis. Prentice Hall, Englewood Cli€s, NJ. Berryman, J.G., 1980. Con®rmation of Biot's theory. Appl. Phys. Lett. 37 (4), 382±384. Biot, M.A., 1941. General theory of three-dimensional consolidation. J. Appl. Phys. 12, 155±164. Biot, M.A., 1956. Theory of propagation of elastic waves in a ¯uid saturated porous solid. J. Acoust. Soc. of America 28 (2), 168±191. Biot, M.A., 1962a. Mechanics of deformation and acoustic propagation in porous media. J. Appl. Phys. 33, 1482±1498. Biot, M.A., 1962b. Generalized theory of acoustic propagation in porous dissipation media. J. Acoust. Soc. of America 34, 1254±1264. Bowen, R.M., Chen, P.J., 1975. Waves in a binary mixture of linear elastic material. Journal de Mecanique 14, 237±266. Bowen, R.M., Reinicke, P.J., 1978. Planar progressive waves in a binary mixture of linear elastic material. J. Appl. Mech. 45, 493±499. Brooks, A.N., Hughes, T.J., 1982. Streamline upwind/Petrov±Galerkin formulations for convection dominated ¯ows with particular emphasis on the incompressible Navier±Stokes equations. Comp. M. Appl. Mech. Eng. 32, 199±259. Gajo, A., 1995. The in¯uence of viscous coupling in the propagation of elastic waves in saturated soil. J. Geotech. Eng., ASCE 121 (9), 636±644. Gajo, A., 1996. The e€ects of inertial coupling in the interpretation of dynamic soil tests. Geotechnique 46 (2), 245±257. Gajo, A., MongiovõÁ, L., 1994. The e€ects of measure accuracy in the interpretation of dynamic tests on saturated soils. In: Proc. Int. Symp. on Pre-failure Deformation Characteristics of Geomaterials, Hokkaido. Gajo, A., MongiovõÁ, L., 1995. An analytical solution for the transient response of saturated porous elastic solids. Int. J. Numer. Anal. Methods Geomech. 19, 399±413. Gajo, A., Saetta, A., Vitaliani, R., 1994. Evaluation of three and two ®nite element methods for the dynamic response of saturated soils. Int. J. Num. Meth. Eng. 37, 1231±1247. Gajo, A., Fedel, A., MongiovõÁ, L., 1997. Experimental analysis of the e€ects of ¯uid-solid couplings on the velocity of elastic waves in saturated porous media. Geotechnique 47 (5), 993±1008. Garg, S.K., Nayfeh, H., Good, A.J., 1974. Compressional waves in ¯uid-saturated elastic porous media. J. Appl. Phys. 45, 1968±1974. Hamdi, F., Taylor Smith, D., 1982. The in¯uence of permeability on the compressional wave velocity in marine sediments. Geoph. Prosp. 30, 622±640. Hiremath, M.S., Sandhu, R.S., Morland, L.W., W. E. Wolfe, W.E., 1988. Analysis of one-dimensional wave propagation in a ¯uid-saturated ®nite soil column. Int. J. Numer. Anal. Methods Geomech. 12, 121±139. Hueckel, T., Pellegrini, R., Del Olmo, C., 1998. A constitutive study of thermo-elasto-plasticity of deep carbonatic clays. Int. J. Numer. Anal. Methods Geomech. 22, 549±574. Hughes, T.J.R., 1987. The Finite Element Method. Prentice Hall, Englewood Cli€s, NJ. Johnson, D.L., Plona, T.J., Scala, C., 1982. Tortuosity and acoustic slow waves. Phys. Rev. Lett. 49 (25), 1840±1844. McTigue, D.F., 1986. Thermoelastic response of ¯uid-saturated porous rock. J. Geophys. Res. 91 ( B9), 9533±9542. Newmark, N.M., 1959. A method for computation for structural dynamics. Proc. A.S.C.E. 85 (EM1), 67±94. Nur, A., Byerlee, J.D., 1971. An exact e€ective stress law for elastic deformation of rock with ¯uids. J. Geophys. Res. 76 (26), 6414±6419. Pecker, C., Deresiewicz, H., 1973. Thermal e€ects on wave propagation in liquid-®lled porous media. Acta Mech. 16, 45±64. Plona, T.J., 1980. Observation of a second bulk compressional wave in a porous medium at ultrasonic frequencies. Appl. Phys. Lett. 36 (4), 259±261.

344

A. Gajo / International Journal of Plasticity 18 (2002) 313±344

Saada, R.A., Bonnet, G., Bouvard, D., 1996. Thermomechanical behaviour of casting sands: experiments and elastoplastic modelling. Int. J. Plasticity 12 (3), 273±294. Shao, J.F., 1997. A numerical solution for a thermo-hydro-mechanical coupling problem with heat convection. Int. J. Rock Mech. Min. Sci. 34 (1), 163±166. Van Der Grinten, J.G.M., Van Dongen, M.E.H., Van der Kogel, H., 1985. A shock-tube technique for studying pore-pressure propagation in a dry and water-saturated porous medium. J. Appl. Phys. 58, 2937±2942. Van Der Grinten, J.G.M., Van Dongen, M.E.H., Van der Kogel, H., 1987. Strain and pore pressure propagation in a water-saturated porous medium. J. Appl. Phys. 62, 4682±4687. Zienkiewicz, O.C., Shiomi, T., 1984. Dynamic behaviour of saturated porous media: the generalised Biot formulation and its numerical solution. Int. J. Numer. Anal. Methods Geomech. 8, 71±96. Zienkiewicz, O.C., Taylor, R.L., 1989. The Finite Element Method, Vol I., 4th Edition. McGraw-Hill.