## Abstract

This paper examines the coupled hydro-thermo-mechanical behaviour of a fluid-saturated porous sphere with a skeletal fabric that can exhibit either elastic or elasto-plastic mechanical behaviour. Analytical results for the thermo-poro-elastic response of the sphere subjected to transient heat transfer are complemented by computational results for the analogous thermo-poro-elasto-plastic problem. The results presented in the paper examine the influence of the permeability, thermal expansion properties of the pore fluid and the skeleton, and the elasto-plasticity effects of the porous skeleton on the time-dependent pore fluid pressure, displacement and stress within the sphere.

## 1. Introduction

An understanding of the thermo-hydro-mechanical (THM) behaviour of geomaterials is necessary for many applications in the environmental geosciences and geoenvironmental engineering. In the area of deep geological disposal of heat-emitting radioactive wastes (Gnirk 1993; Selvadurai & Nguyen 1997; Rutqvist *et al.* 2001; Stephansson *et al.* 2004; Alonso *et al.* 2005), earthquake hazards along fault zones (Rice 2006), geothermal energy extraction (Dickson & Fanelli 1995; Duffield & Sass 2003) and ground source energy extraction (Knellwolf *et al.* 2011; Laloui & Di Donna 2011), the coupled processes associated with heat transfer, mechanical deformations of the porous skeleton and fluid movement under thermal, hydraulic and mechanical effects are recognized as important considerations. Knowledge of the coupled THM processes in the fluid-saturated geological medium is also essential to accurately evaluate the pore-fluid pressure responses, where adverse fluid pressures can induce damage in the form of micro-fracturing to the porous skeletal structure. Investigations dealing with the modelling of THM processes in poro-elastic fluid-saturated media are quite extensive. The reader is directed to recent literature in geomechanics, particularly those conferences devoted to poromechanics (Selvadurai 1996; Thimus *et al.* 1998; Auriault *et al.* 2002; Abousleiman *et al.* 2005), articles with extensive reviews (Selvadurai 2007; Schanz 2009; Selvadurai & Selvadurai 2010) and the recent articles on nuclear waste management cited previously, for further information.

Conventional poro-elasticity and thermo-poro-elasticity problems have been examined in the literature, starting with the classical studies by Mandel (1953) and Cryer (1963) (see also Selvadurai & Shirazi 2005), and the amplification of the pore-fluid pressure during radial loading of a saturated clay was confirmed in the experimental studies conducted by Gibson *et al.* (1963). The thermo-poro-elastic problem for a hollow sphere subjected to a sudden rise in temperature and pressure on its inner wall was studied by Kodashima & Kurashige (1996). In their work, the heat transport equation included a nonlinear convective term; the analytical solution of the problem was obtained for the steady-state case, and the transient solution was obtained numerically. Similar studies were conducted by Rehbinder (1995), who considered problems with cylindrical and spherical symmetries, and the stationary solutions for the nonlinear thermo-poro-elastic problem were obtained using a perturbation technique. Belotserkovets & Prevost (2011) studied the transient response of a thermo-poro-elastic sphere subjected to an applied radial stress. The goal of their studies was to characterize the influence of the applied stress on the change in the temperature within the sphere.

The present paper includes the study of an idealized problem where the skeleton of the fluid-saturated material exhibits elasto-plastic effects. The treatment of such effects is not generally considered applicable to brittle geomaterials that have a greater susceptibility for damage development in the form of micro-cracks and micro-fissures that can alter the stiffness and fluid transmissivity properties (Selvadurai & Mahyari 1997; Mahyari & Selvadurai 1998; Selvadurai 2004; Selvadurai & Shirazi 2004, 2005; Shirazi & Selvadurai 2004). The consideration of plasticity effects is largely of interest to soft rocks that could display elasto-plastic constitutive behaviour in their skeletal responses, particularly in the small strain range (Desai & Siriwardane 1984; Darve 1990; Davis & Selvadurai 2004; Pietruszczak 2010). Thermo-poro-elasto-plasticity of clays has been discussed in connection with engineered clay barriers that have been proposed for high-level nuclear waste disposal endeavours (Hueckel & Borsetto 1990; Hueckel & Peano 1996; Laloui *et al.* 2005; Sanchez *et al.* 2008). Analytical investigations are rare; an example of the consolidation of a poro-elasto-plastic column is given by Pariseau (1999) and the Cam Clay plasticity model applied to the problem of thermal failure in saturated clays was considered by Hueckel *et al.* (2009). The Cam Clay model was also used to obtain an analytical solution for the undrained expansion of a spherical cavity (Silvestri & Abou-Samra 2011).

This paper presents a comparison of analytical and computational results for the idealized problem of the external heating of a fluid-saturated sphere with a skeletal response that can be described by elastic or elasto-plastic phenomena. The objectives of the computational studies presented in the paper are to develop a set of benchmark computational results that can be used to validate computational codes (such as COMSOL and ABAQUS) in their ability to provide inter-code validations.

## 2. Governing equations

We consider the class of porous materials in which the pore space is fluid-saturated and the skeleton can exhibit mechanical and thermal deformations associated with stresses and temperatures applied to the skeletal phase. The porous material undergoes failure, which can be described by an appropriate elasto-plastic model, characterized by an yield condition, a flow rule and a hardening law. It is assumed that the thermo-physical properties for the solid phase of the geomaterial are the same, irrespective of whether it exhibits elastic or plastic responses. Similarly, we assume that the fluid transport properties in the porous medium are the same in both the intact and yielded regions. We denote the total stress tensor in the pore skeleton by *σ*_{ij}(*x*_{i},*t*) and the pore-fluid pressure by *p*(*x*_{i},*t*). The development of the poro-elastic part of the constitutive modelling adopts the procedures described by Selvadurai & Nguyen (1995). Briefly, the constitutive equation governing the poro-elastic response of the fluid-saturated porous medium is given by
2.1where *δ*_{ij} is Kronecker’s delta and
2.2is the infinitesimal strain tensor expressed in terms of the displacement components *u*_{i}; *T* is the temperature; *K*_{D} and *G*_{D} are the bulk and shear moduli of the porous skeleton, *α*_{s} is the coefficient of linear thermal expansion of the solid phase and the comma denotes the partial derivative with respect to the spatial variable *x*_{i}. In (2.1), *α* is Biot’s coefficient defined by
2.3where *K*_{s} is the bulk modulus for the solid phase. In most geomechanical applications, , with the result that *α*=1. In the absence of body forces and gravity effects, the equation of equilibrium expressed in terms of the total stress
2.4can be combined with the constitutive relationship (2.1) to give
2.5where ∇^{2} is Laplace’s operator. The flow of water through the porous medium is described by the conventional form of Darcy’s law applicable to an isotropic porous medium: i.e. the fluid velocity in the pore space ( relative to the porous skeleton velocity ( is given by
2.6where *K* is the permeability and *η* is the dynamic viscosity. The mass conservation equation for the fluid-saturated poro-elastic geomaterial can be written as
2.7where *α*_{s} is the coefficient of linear thermal expansion of the solid phase, *α*_{f} is the coefficient of linear thermal expansion of the pore fluid, *u*_{k,k}(=*ε*_{vol}) is the volumetric strain of the skeleton, which can include both elastic and plastic components, *K*_{f} is the bulk modulus of the fluid and *φ* is the porosity. Equation (2.7) also holds for elasto-plastic porous materials with incompressible and elastic solid phases.

When dealing with quasi-static thermo-poro-elasto-plastic materials, we can assume that neither the deformations of the medium nor the fluid flow process results in heat generation, and therefore does not contribute to any change in temperature within the medium. This considerably simplifies the equations governing the heat transfer process, which can be expressed as a classical heat conduction equation of the form
2.8where *k*_{c}, *ρ* and *c*_{p} are, respectively, the effective values of the thermal conductivity, mass density and heat capacity of the saturated porous medium. These can be expressed as the averages of the properties of the solid (subscript s) and fluid (subscript w) phases: *k*_{c}=*φk*_{cw}+(1−*φ*)*k*_{cs} and *ρc*_{p}=*φ*(*ρc*_{p})_{w}+(1−*φ*)(*ρc*_{p})_{s}.

In addition to the elastic response associated with the porous skeleton, we assume that the skeleton can exhibit failure upon attainment of a stress state prescribed by a yield condition. The choice of a yield condition is dependent on the type of geomaterial and, for the purposes of illustration, we consider the modified Cam-Clay model (Schofield & Wroth 1968; Desai & Siriwardane 1984; Davis & Selvadurai 2004; Pietruszczak 2010). The yield surface is given by
2.9where *q* is the von-Mises stress, *a* is the radius of the yield surface, *M* is the slope of the critical state line and is the isotropic stress applicable to the skeletal stress: i.e.
2.10The centre of the yield surface (*a*,0) in the plane can be expressed as , where is the initial yield stress for the isotropic compression stress state, and
2.11is the hardening law that prescribes the dependence of the isotropic stress on the volumetric plastic strain.

The yield condition (2.9) can be written as
2.12To determine the incremental plastic strains, we specify an associated flow rule such that
2.13The plastic multiplier d*λ* can be found from the consistency condition
2.14Substituting the flow rule (2.13) into the consistency condition (2.14), the plastic multiplier can be obtained as
2.15where *D*_{ijkl}=(*K*_{D}−2/3*G*_{D})*δ*_{ij}*δ*_{kl}+*G*_{D}(*δ*_{ik}*δ*_{jl}+*δ*_{il}*δ*_{jk}) is the elastic constitutive tensor.

The incremental plastic strains can now be obtained as
2.16The incremental form of the elasto-plastic constitutive tensor defined as can be expressed as
2.17The hardening rule adopted here assumes that
18where *H* is a positive constant. If the volumetric strain is negative (compaction), hardening takes place and the yield stress also increases. In this case, the loading path in the -plane intersects the yield surface when the effective pressure is larger than the radius of the yield surface. On the other hand, if the volumetric strain is positive (dilatation), the yield stress decreases. In this case, the loading path in the -plane intersects the yield surface when the effective pressure is smaller than the radius of the yield surface.

## 3. Analysis of the thermo-poro-elastic sphere

Ideally, the results derived from the computational modelling should be validated through comparisons with either exact or approximate mathematical solutions of the governing linear or nonlinear partial differential equations applicable for both the thermo-poro-elastic and thermo-poro-elasto-plastic material responses. Such solutions are either unavailable or available in forms that are not amenable for convenient calibration of the computational codes. Some progress can be made in the analysis of the thermo-poro-elasticity problem related to *boundary* heating of a fluid-saturated poro-elastic sphere and the problem of *uniform* heating of a poro-elasto-plastic sphere. The basic approach for the poro-elastic sphere will be described in this section.

We consider the problem of a fluid-saturated poro-elastic sphere of radius *R*_{0} that is subjected to boundary heating, i.e. a prescribed temperature change at the boundary. Because the problem related to the heating of the sphere exhibits spherical symmetry, the only non-zero primary-dependent variables are the radial displacement *u*(*R*,*t*), the pore-fluid pressure *p*(*R*,*t*) and the temperature *T*(*R*,*t*). The solution of the problem of the sphere subjected to the temperature change at the boundary will allow us to derive solutions for other (related) problems such as a sphere subjected to radial stress at the boundary, and also to the uniform temperature change applied to the entire sphere.

The equilibrium equation (2.5) in the spherical coordinate system takes the form
3.1The equilibrium equation (3.1) can also be written in terms of the volumetric strain *e*=(∂*u*/∂*R*)+(2*u*/*R*) as
3.2From (3.2), the pressure can be expressed as
3.3where *f*(*t*) is an arbitrary function of time.

The fluid flow equation (2.7) is given by
3.4where *M*_{B}=[(*φ*/*K*_{f}+(*α*−*φ*)/*K*_{s})]^{−1} is the Biot modulus.

Now assume that the fluid and solid phases are incompressible. By taking the divergence of the equilibrium equation (3.2) and substituting ∇^{2}*p* into the fluid flow equation (3.4), we can obtain the governing equation for the volumetric strain,
3.5In the case of compressible phases, the equation (3.5) will be coupled with the fluid pressure (Belotserkovets & Prevost 2011).

The temperature satisfies the heat conduction equation 3.6A solution to the governing equations (3.5) and (3.6) for the volumetric strain and temperature should satisfy the boundary and initial conditions applicable to a problem.

The boundary of the sphere is subjected to the temperature change *T*_{0}, the fluid pressure and the radial stress at the boundary are maintained at zero values: i.e.
3.7It was shown by Mason *et al.* (1991) that the stress boundary condition can be written in terms of the volumetric strain as
3.8At time *t*=0, all the dependent variables within the sphere are assumed to be zero.

The solution of the heat conduction equation (3.6) can be obtained by using a Laplace transform technique; the Laplace transform of the temperature field can be found as
3.9By using the Cauchy residue theorem, the expression (3.9) can be inverted and the temperature field can be obtained in the form
3.10In the limit as *R*→0, .

By taking the Laplace transform of the governing equation for the volumetric strain (3.5), we obtain
3.11where the transform of the temperature field (3.9) has been used. The unknown function *A*(*s*) can be found by enforcing the zero-traction boundary condition (3.7), and the result can be written as
3.12The Laplace transform of the pressure can be found by using (3.3) and a null boundary condition for the pressure (3.7)
3.13The radial displacement can be found using the formula (Mason *et al.* 1991)
3.14which gives
3.15Using the final value theorem, it can be shown that the displacement distribution within the sphere in the long term is given by
3.16which is independent of the thermal expansion of the fluid phase. When applying the final value theorem, the useful relations are and .

In order to find the inverse Laplace transform of , it is important to note that the zeros of *F*(*sω*), i.e. *F*(*s*_{n}*ω*)=0, are among the poles of the first term of the expression for ; these coincide with the poles found by Mason *et al.* (1991) for the problem of a sphere subjected to a *mechanical load* and can be found by solving the following equation
where *ν*_{D} is the skeletal Poisson’s ratio. When , *x*_{n}→(*nπ*)^{2}. The other poles are *s*_{n}=−(*n*^{2}*π*^{2}/*κ*) and *s*=0. By applying the Cauchy residue theorem, the expression for the volumetric strain can be obtained as
3.17where
3.18The fluid pressure can be obtained from (3.13) using the volumetric strain given in (3.17) and the temperature field (3.10),
3.19

### Remark 3.1

From the stress–strain constitutive equations, the mean effective stress is equal to
3.20On the external surface, *R*=*R*_{0}, the radial stress and the fluid pressure are equal to zero, and therefore the effective hoop stress can be obtained from (3.20) as
3.21Using the stress–strain relationships, we obtain
3.22and the strains on the external surface can be expressed in terms of the volumetric strain *e* as
3.23Because *u*(*R*,*t*)=*ε*_{φφ}(*R*,*t*)*R*, the expression for the hoop strain *ε*_{φφ}(*R*_{0},*t*) enables us to find the surface displacement in terms of the volumetric strain,
3.24The *initial* volumetric strain on the external surface *R*=*R*_{0} can be derived from (3.11), using the initial value theorem,
3.25Substitution of this value into (3.23) allows us to find the initial strains on the external surface,
3.26This shows that the initial value of the surface displacement *u*(*R*_{0},0) is zero. In the long term , the sphere is stress-free, and the displacement of the poro-elastic sphere is equal to that of the elastic sphere, i.e. *u*(*R*)=*α*_{s}*T*_{0}*R*.

Again, using the expression for the initial volumetric strain (3.25) in (3.20), the initial value of the hoop stress on the external surface can be obtained as 3.27The initial skeletal hoop stress on the external surface is compressive if the temperature change is positive. In the long term, the sphere is stress-free.

### Remark 3.2

We now consider a problem when the initial temperature is non-zero. At first, it is sufficient to assume that *T*(*R*,*t*)≡*T*_{C}, constant, for all *R* and *t*, and then use the superposition principle, if necessary. For this particular case, the volumetric strain can be obtained from (3.17) by taking the limit *κ*→0 and setting *T*_{0}=*T*_{C},
3.28The fluid pressure can still be found from (3.19) but with *T*(*R*,*t*)=*T*_{C}. It can be shown that the initial value of the fluid pressure *p*(*R*,0)=*φK*_{D}[3*α*_{f}−3*α*_{s}]*T*_{C}, 0≤*R*<*R*_{0}, derived from (3.19), is in fact equal to the pressure in the undrained sphere, when the change in the mass of the fluid is zero. However, for a sphere with a drained boundary condition, the same response can be observed initially at *t*=0 for 0≤*R*<*R*_{0}, if the temperature inside the sphere *instantaneously* becomes equal to *T*_{C}. For this case, the pressure distribution inside the sphere exhibits a Mandel–Cryer-type effect (Mason *et al.* 1991). (In general, the Mandel–Cryer-type effect can take place in a poro-elastic sphere when there is a sudden change in the fluid-flow boundary condition on the external surface from an initially undrained state to a drained state.)

Let us find initial values of strains and stresses inside a sphere that is heated uniformly and instantaneously. From (3.5), the initial value of the volumetric strain is given by
3.29Because the heating is uniform, we can assume that the displacement depends linearly on *R*, and therefore the components of the strain tensor are equal, i.e.
3.30From the stress–strain constitutive relations, analogous to (3.22), the components of the stress tensor can be obtained as
3.31At the boundary *R*=*R*_{0}, the initial hoop strain is given by the same expression owing to continuity of the displacement, i.e.
3.32The radial component of the strain tensor can be found from the constitutive relations (3.22) by setting the radial stress to zero, . This gives
3.33We can then compute the initial hoop stress from (3.22)
3.34The initially drained condition is another admissible set of *initial* conditions. For this case,
3.35

### Remark 3.3

Suppose that the thermal expansion coefficients and temperature change in (3.28) are *K*_{D}3*α*_{s}*T*_{C}=−*P*_{0}; *φ*/(1−*φ*)*K*_{D}3*α*_{f}*T*_{C}=*P*_{0}, where *P*_{0} is a positive constant. Then the expression (3.28) can be simplified to
3.36where *ν*_{D} is the Poisson’s ratio for the undrained material. In deriving (3.36), we have used the equation for the roots *x*_{n}. Equation (3.36) corresponds to the equation derived by Mason *et al.* (1991) for a poro-elastic sphere subjected to a compressive radial stress −*P*_{0} at the boundary *R*=*R*_{0}.

### Remark 3.4

It is instructive to solve the problem for the undrained poro-elastic sphere subjected to the temperature change *T*(*R*_{0},*t*)=*T*_{0} on the external surface. The volumetric strain for the undrained sphere can be obtained from (3.5) by setting the permeability equal to zero,
3.37The fluid pressure on the external surface can be found from (3.8)
3.38Using the expression (3.37) for the volumetric strain, the fluid pressure on the external surface *R*=*R*_{0} can be evaluated from (3.38) as
3.39This expression can be substituted into (3.19) to obtain the fluid pressure for all points within the sphere.

The fluid pressure in the undrained sphere subjected to a uniform and constant temperature change, i.e. *T*(*R*,*t*)≡*T*_{C}, can be found from (3.39) by setting *T*_{0}=*T*_{C} and ,
3.40It is clear that if the temperature *T*_{C}>0, the thermal expansion of the solid phase induces negative fluid pressure, and the thermal expansion of the fluid phase results in positive fluid pressure.

### Remark 3.5

Consider the case where the poro-elastic sphere with incompressible constituents is subjected to a non-zero fluid pressure and zero temperature change on the external surface. For this type of loading, the resulting strains will be zero and no change in the mass of the fluid will occur. The fluid pressure within the sphere will be equal to the applied pressure at the boundary *at all times*.

### Remark 3.6

The inverse Laplace transform can also be found numerically by using an expansion involving Legendre polynomials. Suppose that *F* is the Laplace transform of a function *f*, i.e.
3.41Then the original function *f*(*t*) can be found using the Legendre polynomial expansion method,
3.42where *P*_{n} is the Legendre polynomial of order *n*. The first five Legendre polynomials are: *P*_{0}=1, *P*_{1}(*x*)=*x*, *P*_{2}(*x*)=1/2(3*x*^{2}−1), *P*_{3}(*x*)=1/2(5*x*^{3}−3*x*), *P*_{4}(*x*)=1/8(35*x*^{4}−30*x*^{2}+3) and *P*_{5}(*x*)=1/8(63*x*^{5}−70*x*^{3}+15*x*).

Special attention is required when two time constants *ω* and *κ* are equal, i.e. *ω*=*κ*, because the particular solution of the inhomogeneous equation (3.5) needs to be modified. However, a good approximation to the actual solution in this case can be obtained by setting *ω*≈*κ*.

## 4. Results for the thermo-poro-elastic sphere

To validate the analytical solution for the thermo-poro-elastic sphere presented in §3, we consider a thermo-poro-elastic sphere of radius *R*_{0}=10 m subjected to a temperature rise of 100^{°}C on the external surface. We implicitly assume that although the temperature of the fluid is set to 100^{°}C, we do not account for any phase transformations during heating. The initial temperature is assumed equal to zero. The fluid pressure and the radial stress at the boundary are set to zero. The results were evaluated analytically (using the closed-form analytical expressions (3.17)–(3.18) for the volumetric strain, (3.19) for the fluid pressure, (3.21) for the hoop stress and (3.24) for the surface displacement) and numerically (using the finite-element program COMSOL). The problem can be solved using a two-dimensional axisymmetric formulation. Figure 1 shows the projection of the sphere onto the *x*_{1}−*x*_{2} plane of the Cartesian coordinate system and the finite-element mesh generated by COMSOL. Appropriate zero-boundary conditions were imposed on the displacements, fluid velocity and tractions at the relevant boundaries of the domain to ensure symmetry.

The mechanical and physical parameters of the thermo-poro-elastic problem are specified as follows: Young’s modulus of the drained skeleton *E*=60×10^{9} Pa; Poisson’s ratio *ν*=0.30; viscosity of the pore fluid *η*=0.001 Pa s; permeability *K*=3×10^{−19} or 3×10^{−20} m^{2}; thermal conductivity *k*_{c}=3.15 W m^{−1} ^{°}C^{−1}; specific heat of the solid material *c*_{s}=700 J kg^{−1} ^{°}C^{−1}; density of the solid material *ρ*_{s}=2700*kg* m^{−3}; specific heat of the fluid *c*_{w}=4190 J kg^{−1} ^{°}C^{−1}; density of the fluid *ρ*_{w}=1000 kg m^{−3}; coefficient of linear thermal expansion of dry skeleton *α*_{s}=8.3×10^{−6} ^{°}C^{−1}; coefficient of linear thermal expansion of fluid *α*_{f}=69×10^{−6} ^{°}C^{−1}; porosity *φ*=0.25. The solid material composing the porous skeleton and the fluid are assumed to be incompressible.

The time-dependent evolution of the temperature in the sphere can be obtained using the analytical expression (3.10). The temperature is initially equal to zero at all points, except at the external surface, where it is prescribed equal to *T*_{0}=100^{°}C, and it then gradually increases reaching the maximum value of 100^{°}C.

The time-dependent evolution of the volumetric strain inside the poro-elastic sphere is shown in figure 2. The volumetric strain is plotted for the centre, the mid-radius and the external surface. The permeability of the sphere is set to *K*=3×10^{−19} m^{2}. The initial value of the volumetric strain is zero at all points inside the sphere except at the boundary, where it is equal to 3*α*_{s}*T*_{0}*K*_{D}/(*K*_{D}+(4/3)*G*_{D})=0.0015. The long-term value of the volumetric strain is equal to that of the elastic sphere, i.e. 3*α*_{s}*T*_{0}=0.0025, for all points.

Figure 3 shows the volumetric strain at the boundary of thermo-poro-elastic sphere for two values of permeability, *K*=3×10^{−19} m^{2} and *K*=3×10^{−20} m^{2}. If the permeability is *K*=3×10^{−19} m^{2}, the volumetric strain monotonically increases as a function of time. If the permeability is an order of magnitude lower (i.e. *K*=3×10^{−20} m^{2}), the evolution of the volumetric strain is not monotonic and the transient value of the strain exceeds its steady-state value, i.e. 3*α*_{s}*T*_{0}=0.0025. This behaviour does not occur in the thermo-elastic sphere and is a consequence of the thermal expansion of the fluid that flows out of the sphere very slowly because of the low permeability.

Figure 4 displays the evolution of the hoop stress at the boundary of the thermo-poro-elastic sphere. As before, two values of permeability are used, *K*=3×10^{−19} m^{2} and *K*=3×10^{−20} m^{2}. From (3.27), the initial value of the hoop stress is given by −6*α*_{s}*T*_{0}*K*_{D}*G*_{D}/(*K*_{D}+(4/3)*G*_{D})=−71.1 MPa. The hoop stress is a monotonic function of time and entirely compressive if the permeability is *K*=3×10^{−19} m^{2}. However, if the permeability is *K*=3×10^{−20} m^{2}, the hoop stress is not a monotonic function and can become tensile for some period of time during the deformation. This phenomenon can be explained by the thermal expansion of the fluid and cannot be captured by the thermo-elastic model. The steady-state value of the stress is zero.

The surface displacement of the poro-elastic sphere is depicted in figure 5 for two permeability values, *K*=3×10^{−19} m^{2} and *K*=3×10^{−20} m^{2}. The initial displacement is zero and the final displacement is equal to that of the elastic sphere, *α*_{s}*T*_{0}*R*_{0}=0.0083 m. If the permeability of the sphere is small, the surface displacement is not a monotonic function of time and exceeds its steady-state value for some period of time.

A striking feature of the poro-elastic solution is the presence of the fluid pressure as a dependent variable. Figure 6 shows the fluid pressure at the centre of the poro-elastic sphere for two permeabilities *K*=3×10^{−19} m^{2} and *K*=3×10^{−20} m^{2}. The fluid pressure is shown to be smaller in the sphere with the greater permeability. The fluid pressure is a non-monotonic function of time for both values of the permeability. The initial and the final values of the fluid pressure are zero, consistent with the diffusive character of the pore pressure variation.

Figure 7 shows the influence of the thermal expansion coefficient of the solid and fluid phases on the fluid pressure at the centre of the poro-elastic sphere. The thermal expansion of the solid phase is *α*_{s}=0 or *α*_{s}=8.3×10^{−6} ^{°}*C*^{−1}, and the thermal expansion of the fluid phase is . The permeability is set to *K*=3×10^{−20} m^{2}. By comparing the two graphs, it is apparent that the thermal expansion of the solid phase gives rise to negative fluid pressure, while the thermal expansion of the fluid phase leads to positive pressure.

Figure 8 shows the fluid pressure evolution at the centre of two poro-elastic spheres subjected to a *uniform* temperature change *T*(*R*,*t*)=100^{°}C for all *R* and *t*. It is assumed that this temperature change is applied instantaneously as a Heaviside step function of time. For the first sphere, the fluid pressure at the boundary is zero (drained sphere), and for the second sphere, the velocity at the boundary is zero (undrained sphere). For the latter case, from (3.40), the pressure is constant and equal to 3*φK*_{D}(*α*_{f}−*α*_{s})*T*=227 MPa, if the porosity is *φ*=0.25. This expression shows that the thermal expansions of the phases result in fluid pressures of opposite signs. In the short term, the drained sphere experiences the Mandel–Cryer effect (Mason *et al.* 1991) because the fluid pressure inside the drained sphere in fact turns out to be larger, in absolute value, than the fluid pressure inside the undrained sphere.

## 5. Results for the thermo-poro-elasto-plastic sphere

In this section, we examine the time-dependent response of a poro-elasto-plastic sphere of radius *R*_{0} subjected to a sudden positive temperature change *T*_{0} on the external surface. The fluid pressure and applied radial stress at the boundary are set to zero.

The elasto-plastic behaviour of the skeleton of the geomaterial is described by the isothermal-modified Cam-Clay model. Since this model is not implemented in COMSOL, the finite-element code ABAQUS was used instead, since the modified Cam-Clay plasticity model is in the library of plasticity models implemented in ABAQUS.

The initial yield stress or initial size of the yield surface is assumed to be or 50 MPa; further evolution of the yield surface is controlled by the specified hardening rule GPa; the parameter *M* in (2.9) is set to 1.

Figure 9 displays two yield surfaces of the modified Cam-Clay model for the values of the initial yield stress equal to 40 and 50 MPa. The yield surfaces are plotted in the effective stress space and in fact correspond to any axisymmetric stress state irrespective of the values of the temperature or fluid pressure. The equation for the yield surface illustrated in figure 9 can be obtained from (2.12) as 5.1We immediately observe from figure 9 that the material cannot withstand large tensile stresses but the range of admissible compressive stresses is much broader. In fact, it can be shown that the tensile radial stresses lie in the range [0, ], but the compressive radial stresses can change from to 0. These bounds can of course be changed owing to material hardening or softening but we do not expect these changes to be significant. Therefore, we can place a bound on the admissible compressive and tensile radial stresses that can be applied to the surface of the elasto-plastic sphere.

The yield surface always passes through the origin (0,0) regardless of the changes in the yield stress; thus we can consider this point to be stationary. The yield surface intersects the hoop stress axis at two points and .

For calibrating a numerical solution of the problem of *transient* heat transfer, it is useful to consider first the porous elasto-plastic sphere subjected to the *uniform* temperature change [0,*T*_{C}] applied instantaneously to the entire sphere because, for this case, the *analytical solution* can be obtained. For this case, according to (3.31) and (3.34), the initial values of stresses for the poro-elastic sphere are
5.2It is clear that the loading path for the elastic sphere makes an acute angle with the outward normal to the yield surface at the zero stress point if
5.3From (5.2), this inequality holds if *α*_{f}>*α*_{s}. Therefore, given that *α*_{f}>*α*_{s}, the loading path for the elasto-plastic sphere cannot be within the yield surface and plastic deformations must be involved; on the other hand, the loading path cannot go outside the yield surface because the zero stress point is stationary and there is no expansion or contraction of the yield surface at this point. Hence, the loading path for the elasto-plastic sphere must either remain at the zero stress point or continue to move along the yield surface. Let us assume that and show that this is in fact the correct solution if *α*_{f}>*α*_{s}.

Consider a porous elasto-plastic sphere with *α*_{f}>*α*_{s} subjected to the uniform temperature change [0,*T*_{C}] applied instantaneously. Suppose that the effective stresses are zero, i.e. . From the equilibrium equations (2.4), it follows that the fluid pressure must be uniform; however, since *p*=0 at the boundary, the fluid pressure must be zero throughout. The traction boundary condition *σ*_{RR}(*R*_{0},*t*)=0 is trivially satisfied. The fluid flow equation (2.7) can be satisfied if we set the volumetric strain equal to
5.4Owing to the uniform character of the stress state, the strain components are equal and can thus be obtained from (5.4) as
5.5The elastic components of the strain tensor can be derived from the stress–strain relationships (3.22), using the fact that the effective stresses are zero. This gives
5.6Therefore, the plastic strains are given by
5.7If, however, *α*_{f}≤*α*_{s}, the plastic strains are initially equal to zero, and the loading path must lie within the yield surface; the solution sought has thus been obtained. Clearly, this solution is not dependent on the value of the yield stress and hardening.

Consider now the original problem of the porous elasto-plastic sphere subjected to the temperature change *T*_{0} prescribed on the external surface. The radial stress and fluid pressure *p* are zero at the boundary. It is clear that the loading path for the external surface will always lie on the axis. The initial hoop stress for the elastic range is given by (3.27), and the temperature that corresponds to the onset of yielding on the external surface can be derived from
5.8which for an initial yield stress of 40 MPa, gives
5.9Figure 9 shows the loading path for the external surface when the applied temperature *T*_{0}=20 ^{°}C. Because this temperature is smaller than the critical temperature (5.9), the loading path starts at a point within the yield surface and no yielding will occur in the short term. For *t*>0, the hoop stress starts to decrease until the loading path hits the zero stress point. From this time on, deformation on the external surface takes place at the constant zero stress state, i.e. , and accumulation of the plastic strain takes place at this point. Note that, at this point, the normal to the yield surface always measures 45^{°} to the horizontal axis and therefore, the radial and circumferential components of the plastic strain must be equal, i.e. .

Figure 10 shows the hoop stress–hoop strain curve for the external surface of the porous sphere subjected to the temperature change *T*_{0}=20 ^{°}C at the boundary *R*_{0}=10 m. At time *t*=0, the hoop stress is −14.22 MPa as follows from (3.27). The permeability of the sphere is *K*=3×10^{−19} m^{2}, and therefore, for the analagous poro-elastic sphere, as was shown in figure 4, the hoop stress is a monotonic function. In the long term, this stress reduces to zero and the value of the hoop strain is *α*_{s}*T*_{0}=1.66×10^{−4}. For the poro-elasto-plastic sphere, however, the accumulation of the plastic strain takes place at the point of zero stress, which results in the development of the horizontal segment on the graph. Because the stress is zero on this part of the curve, the increment of the elastic strain is also zero. The final value of the strain for the poro-elasto-plastic sphere is approximately equal to 4.9*E*−4. This result was obtained using the ABAQUS finite-element program. Dashed lines indicate the position of the initial yield surface for .

Figure 11 shows the evolution of the surface displacement of the poro-elasto-plastic sphere subjected to the temperature change *T*_{0}=20^{°}C at the boundary. The permeability of the sphere is *K*=3×10^{−19} m^{2}. Results are shown for two values of the thermal expansion coefficient of the solid phase, *α*_{s}=8.3×10^{−6}^{°}C^{−1} and *α*_{s}=0. The thermal expansion of the fluid phase is *α*_{f}=6.9×10^{−5}^{°}C^{−1}. As was shown in figure 10, deformation on the external surface of the sphere with *α*_{s}=8.3×10^{−6}^{°}C^{−1} is initially elastic until the hoop strain reaches *ε*_{φφ}=*α*_{s}*T*_{0}, and the loading path reaches the zero stress point. After the loading path has reached the zero stress point, the plastic strain begins to accumulate with no further increment in the elastic strain. Therefore, in the long term, the total surface displacement can be written as
5.10where is the plastic hoop strain.

For the elasto-plastic sphere with *α*_{s}=8.3×10^{−6}^{°}C^{−1}, the surface displacement in the long term is and for the elasto-plastic sphere with *α*_{s}=0, . These values can be compared with the displacement of the elasto-plastic sphere heated uniformly and instantaneously to the temperature *T*_{0}=20^{°}C. From (3.32), the surface displacement of such a sphere is equal to *u*(*R*_{0},*t*)=*R*_{0}[*φα*_{f}+(1−*φ*)*α*_{s}]*T*_{0} which gives the value of 0.0047 m if *α*_{s}=8.3×10^{−6}^{°}C^{−1} and 0.00345 m if *α*_{s}=0. Although these displacements are very close to those of the elasto-plastic sphere heated at the boundary, it is worth noting that in the case of the elasto-plastic sphere subjected to a uniform temperature change, the strain is entirely plastic, whereas in the case of the sphere heated at the boundary, part of the total strain, namely *ε*_{φφ}(*R*_{0},*t*)=*α*_{s}*T*_{0}, is elastic. The present results do not depend on the value of the hardening modulus because all plastic deformation is accumulated at the state of the constant zero stress.

Figure 12 shows the fluid pressure at the centre of the elasto-plastic sphere as a function of time. As before, the sphere is subjected to the heating *T*_{0}=20^{°}C at the boundary. The permeability of the sphere is *K*=3×10^{−19} m^{2}. The results are depicted for two values of the thermal expansion of the solid phase, *α*_{s}=8.3×10^{−6}^{°}C^{−1} and *α*_{s}=0. The initial negative values of the fluid pressure are larger for the sphere with *α*_{s}=8.3×10^{−6}^{°}C^{−1}, which supports our assumption that the thermal expansion of the solid phase causes negative fluid pressure (compare this result with the elastic sphere in figure 6). After some time, owing to the thermal expansion of the fluid phase, the fluid pressure becomes positive. The value of the positive fluid pressure for the elasto-plastic sphere is expected to be lower than that for the elastic counterpart owing to the loss of stiffness of the skeletal structure of the material during plastic deformation. In fact, if *α*_{s}=8.3×10^{−6}^{°}C^{−1}, the maximum value of the positive fluid pressure is 4.17 MPa for the *elastic* sphere (this can be deduced from figure 6), and only 0.692 MPa for the *elasto-plastic* sphere (figure 12). The long-term value of the fluid pressure is zero. As before, these results do not depend on hardening but depend on the value of the initial yield stress.

It is important to recognize that the plastic deformations always take place in the given elasto-plastic sphere subjected to heating, no matter how small the applied temperature is. If the temperature is smaller than the critical temperature (5.9), the plastic strain is confined in the short term to a region that does not include the external surface and the centre of the sphere. In the long term, however, the plastic strain becomes non-zero in the entire sphere. The maximum of the plastic strain occurs some distance away from the external surface, approximately equal to 1/4 of the radius. The plastic strain is minimal on the external surface. The plastic deformations can be avoided if the compressive stress is applied on the external surface along with the temperature change.

## 6. Concluding remarks

The coupling of thermo-poro-elasto-plastic effects in fluid-saturated porous media continues to be of interest to many areas of geosciences and geomechanics. This paper presents a benchmark problem that gives a comparison of results for the *transient* heating of poro-elastic and poro-elasto-plastic spheres derived using analytical techniques and two computational codes COMSOL and ABAQUS. The results compare favourably.

The transient solution obtained for the fluid-saturated poro-elastic sphere heated at the boundary shows that the consideration of thermo-poro-elastic effects is important. First of all, poro-elasticity of the fluid-saturated material gives rise to the pore fluid pressure as a dependent variable, which is absent in the analogous thermo-elastic problem. Also, certain parameters, such as the permeability and thermal expansion of the fluid, excluded from the standard thermo-elastic analysis, have a considerable influence on the transient values of the fluid pressure, hoop stress and surface displacement during the heating process.

For the present problem, we showed that the fluid pressure is dependent on the permeability and thermal expansion coefficients of the phases. As expected, for low permeabilities, the fluid pressure is significantly larger. Furthermore, the thermal expansion of the solid phase in general leads to a negative fluid pressure in the short term but the thermal expansion of the fluid results in positive fluid pressure over the longer term. We also observed that for sufficiently small values of permeability, the surface displacement can exceed the steady-state value that is obtained by performing standard thermo-elastic analysis. Tensile hoop stress can also be created on the external surface of the sphere if the permeability is sufficiently small.

The importance of the elasto-plastic effects is emphasized by considering the problem of a porous elasto-plastic sphere heated at the boundary and also the problem of instantaneous uniform heating over the entire domain. The modified Cam-Clay plasticity model is used for the elasto-plastic analysis. For the problem of the uniform heating of a sphere, a simple analytical solution can be constructed, which by itself can serve as a benchmark for the computational codes that include plasticity models. The yield surface plotted in the radial–hoop stress space shows that the range of admissible tensile radial stresses that can be applied safely at the boundary of the elasto-plastic sphere is much smaller than the range of admissible compressive stresses. For the problem of the sphere heated at the boundary, we showed that the loading path for the external surface can initially lie inside the yield surface but then it intersects the yield surface at the stationary zero stress point in which accumulation of the plastic strain takes place. For the elasto-plastic sphere, owing to the thermal expansion of the skeleton and positive volumetric plastic strain, a negative fluid pressure is created in the short term, and, because of the thermal expansion of the fluid, a positive fluid pressure is induced subsequently. The magnitude of the positive fluid pressure in the poro-elasto-plastic sphere is, however, smaller than that in the poro-elastic sphere because the solid phase becomes less stiff in the elasto-plastic sphere.

## Acknowledgements

The work described in this paper was supported by a Natural Sciences and Engineering Research Council of Canada *Discovery Grant* awarded to the first author. The authors are grateful to the comments of the reviewers that resulted in the new developments presented in the paper. The use of the computational codes COMSOL and ABAQUS is purely for demonstration purposes only. Neither the authors nor the research sponsors advocate or recommend the use of these codes without conducting suitable validation procedures that test their accuracy in a rigorous fashion.

- Received January 23, 2012.
- Accepted April 12, 2012.

- This journal is © 2012 The Royal Society