## Abstract

In this paper, we examine the coupled thermo-poromechanical behaviour of a fluid-saturated porous medium of infinite extent bounded internally by a fluid-filled cavity. The mechanical behaviour of the porous skeleton can either be Hookean elastic or elasto-plastic, with a constitutive response corresponding to a modified Cam Clay plasticity model. The fluid within the cavity can be subjected simultaneously to a temperature rise and a pressure pulse. The paper presents analytical results for the spherically symmetric thermo-poroelasticity problem and these are used to validate the thermo-poroelasticity module of a computational code. We proceed to examine the thermo-poroelasto-plasticity problem. Results presented in the paper illustrate the interaction between thermal and mechanical phenomena and their influence on the cavity fluid pressure and the skeletal stresses at the cavity boundary. The paper presents solutions that will be of value in benchmarking exercises.

## 1. Introduction

The thermo-hydro-mechanical (THM) behaviour of geomaterials is important to many areas of environmental geomechanics, including deep geological disposal of heat-emitting radioactive wastes [1–7], frictional heating during earthquake slip [8–10], geothermal energy extraction [11,12], freezing action in soils [13–16] and geomechanics considerations related to geological sequestration of carbon dioxide [17–19]. In many instances, the scope of geoenvironmental engineering problems in these areas is such that recourse must be made to computational methods to solve any realistic environmental geomechanics problems. It is imperative that the computational methods used for these purposes are able to provide solutions that are accurate. One approach that is extensively advocated is the validation of computational approaches through intercode calibrations [6,20]. An alternative is to use certain canonical problems, which can generate analytical solutions albeit largely for the case of thermo-poroelasticity. An example of a recent exercise along these lines is provided by Selvadurai & Suvorov [21], who examined the problem of the boundary heating of a fluid-saturated sphere with skeletal responses that can be either elastic or elasto-plastic. Constitutive theories that describe the mechanical behaviour of the porous skeleton of geomaterials are well documented in the literature [22–24]. The ability of an analytical solution to account for infinite domains is considered a requirement that computational approaches should be able to model accurately, particularly with regard to coupled THM modelling.

In this paper, we examine the poromechanics problem for a medium of infinite extent bounded internally by a fluid-filled spherical inclusion. We show that when the skeletal behaviour of the porous medium is elastic the solution to the thermo-poroelastic problem can be obtained in an analytical form. Examples of the analytical treatment of the isothermal problem of a pressurized cylindrical cavity and a spherical cavity subjected to a constant stress are given by Rice & Cleary [25] and Rice *et al*. [26]. The response of elastic and viscoelastic porous media bounded internally by a spherical cavity and subjected to a constant radial stress was obtained by Spillers [27]. The case of a rigid spherical heat source, with either a pervious or an impervious boundary, contained within a poroelastic medium was examined by Booker & Savvidou [28]. Wu *et al*. [29] obtained an analytical solution to the problem of a cylindrical cavity (wellbore) subjected to a constant fluid pressure and temperature rise on the cavity wall, and non-hydrostatic stresses applied remotely. The thermo-poroelastic problem for a hollow sphere subjected to a sudden rise in temperature and pressure on its inner wall was examined by Kodashima & Kurasighe [30] and the thermo-poroelastic and thermo-poroelasto-plastic problems for a solid sphere are given by Selvadurai & Suvorov [21].

The consideration of plasticity effects is of particular interest to soft geomaterials that could display elasto-plastic constitutive behaviour in their skeletal responses, especially in the small strain range [23,24,31,32]. The irreversible behaviour of the porous skeleton in fluid-saturated porous media subjected to heating that induces elasto-plastic yielding with hardening/softening, including the use of the modified Cam-Clay plasticity model is well documented in geomechanics literature [33,34]. The development of analytical solutions for transient problems in thermo-poroelasto-plasticity is, however, less well developed. The incremental nature of the conventional plasticity theories makes the analytical solution of the governing coupled equations less straightforward. An example of the isothermal consolidation of a poro-elasto-plastic column is given by Pariseau [35]. Further references in this area are given by Selvadurai [36,37] and Selvadurai & Suvorov [21].

The problem of damage development in the porous skeleton of brittle geomaterials and its influence on hydro-thermo-mechanical couplings are more relevant to many geomaterials and examples are given by Selvadurai & Mahyari [38], Mahyari & Selvadurai [39], Selvadurai [40], Selvadurai & Shirazi [41,42] and Shirazi & Selvadurai [43]. Knowledge of the coupled THM processes in a fluid-saturated geologic medium is also essential to accurately evaluate the pore fluid pressure responses, where adverse fluid pressures can induce damage in the form of microfracturing in the porous skeletal structure. Such processes can also contribute to the development of double porosity effects, and the description of THM processes in such materials requires alternative approaches [44]. Investigations dealing with the modelling of THM processes in poroelastic fluid-saturated media are quite extensive. For further information, the reader is directed to recent literature in geomechanics, particularly those conferences devoted to poromechanics and fluid transport [36,45–48], articles with extensive reviews [37,49–51] and articles related to nuclear waste management cited previously.

The focus of this paper is to present analytical results for the development of (i) pressure within a fluid-filled cavity located in a fluid-saturated poroelastic medium, and (ii) the skeletal stresses and displacements at the boundary of the cavity, when the cavity is subjected separately to pressurization and temperature rise. Similar investigations are conducted computationally to examine the response of a fluid-filled cavity contained in a porous medium with a porous skeleton that exhibits elasto-plastic constitutive behaviour. The paper presents a comparison of analytical and computational results for the spherical cavity problem, which can serve to validate the accuracy of computational results that can be obtained from general-purpose computational codes such as COMSOL and ABAQUS.

## 2. Governing equations

Attention is restricted to porous geomaterials that exhibit isotropic mechanical, fluid flow and heat conduction properties. The pore space is assumed to be fully saturated and the skeleton can exhibit mechanical and thermal deformations, associated with stresses and temperatures applied to the skeletal phase; we also assume that the failure of the porous skeleton can be described by an appropriate elasto-plastic model, characterized by a yield condition, a flow rule and a hardening law. The thermo-physical properties for the solid phase of the geomaterial are assumed to be constant, irrespective of whether it exhibits elastic or elasto-plastic responses. In addition, we assume that the fluid transport properties in the porous medium are the same in both the intact and yielded regions.

As the spherical cavity problem examined here exhibits a state of spherical symmetry, we can develop the governing equations of thermo-poroelasticity by considering the spherical polar coordinate system (*R*, *Θ*, *Φ*), where *R* is the radial distance, *Θ* is the polar angle and *Φ* is the azimuthal coordinate. The displacement vector is given by **u**={*u*(*R*,*t*),0,0}^{T}, where *u*(*R*,*t*) is the time-dependent radial displacement and the corresponding infinitesimal strain tensor ** ε**(

*R*,

*t*) and the total Cauchy stress tensor

**(**

*σ**R*,

*t*) are given by 2.1The development of the poroelastic part of the constitutive modelling adopts the procedures described by Selvadurai & Nguyen [52]. The constitutive equation governing the poroelastic response of the fluid-saturated porous medium is given by 2.2where

*p*(

*R*,

*t*) is pore fluid pressure,

**I**is the unit tensor,

*T*is the temperature;

*K*

_{D}and

*G*

_{D}are the bulk and shear moduli of the porous skeleton,

*α*

_{S}is the linear thermal expansion coefficient of the solid phase and

*ε*

_{V}is the volumetric strain. In (2.2),

*α*is the Biot coefficient defined by

*α*=1−

*K*

_{D}/

*K*

_{S}, where

*K*

_{S}is the bulk modulus for the solid phase. If , the Biot coefficient

*α*=1. In the absence of body forces and inertia effects, the equation of equilibrium can be expressed in terms of the total stress tensor in the form [53] 2.3which can be combined with constitutive relationship (2.2) to give 2.4where ∇ is the gradient operator and ∇

^{2}is Laplace's operator. For spherical symmetry, Laplace's operator takes the form 2.5The flow of water through the porous medium is described by the conventional form of Darcy's law applicable to an isotropic porous medium [54,55], i.e. the fluid velocity in the pore space

**v**

^{f}relative to the porous skeleton velocity

**v**

^{s}is given by 2.6where

*K*is the permeability and

*η*is the dynamic viscosity. The mass conservation equation for the fluid-saturated poroelastic geomaterial can be written as 2.7where

*α*

_{s}and

*α*

_{f}are the linear thermal expansion coefficients of the solid phase and pore fluid, respectively,

*K*

_{f}is the bulk modulus of the fluid and

*φ*is the porosity. Equation (2.7) also holds for elasto-plastic porous geomaterials with an incompressible solid phase that in itself behaves elastically.

When dealing with quasi-static deformations in thermo-poroelasto-plastic materials with relatively low permeability (e.g. dense sandstones, limestones and granite; [56–59]), we can assume that neither the deformation rates of the medium nor the fluid flow rates result in heat generation, and therefore do 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 [60]
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 weighted averages of the properties of the solid (subscript ‘s’) and fluid (subscript *w*) phases: i.e. *k*_{c}=*φk*_{cw}+(1−*φ*)*k*_{cs}, *ρ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 yield condition depends on the type of geomaterial and, for the purposes of illustration, we consider the modified Cam-Clay model [22–24]. The yield surface is given by
2.9where *a*_{y} is the radius of the yield surface, *M*_{y} is the slope of the critical state line, is the *effective pressure* and *q* is the *von Mises stress* defined as
2.10Here, is the effective stress and tr denotes the trace of the tensor. The centre of the yield surface (*a*_{y},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.

Using (2.11), yield condition (2.9) can be rewritten 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 flow rule (2.13) into consistency condition (2.14), the plastic multiplier can be obtained as
2.15where **D**^{E} is the elasticity tensor of isotropic materials [53].

The incremental total 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
2.18where *H* is a positive constant. When the volumetric strain is negative (compaction), hardening will take place and the yield stress will increase. In such situations, the load path in the -plane will intersect the yield surface when the effective pressure is larger than the radius of the yield surface. Alternatively, if the volumetric strain is positive (dilatation), then the yield stress will decrease, and the load path in the -plane will intersect the yield surface when the effective pressure is smaller than the radius of the yield surface.

Hardening law (2.18) can be obtained as an approximation to the exponential hardening law
2.19where *e*_{0}=*φ*/(1−*φ*) is the initial void ratio, *λ* and *κ* denote, respectively, the slope of the normal consolidation line and unload–reload compression lines in the isotropic compression test. Typical values of these parameters are *λ*=0.2 and *κ*=0.05.

## 3. Thermo-poroelasticity of the fluid-filled cavity

We examine the problem of a fluid-saturated poroelastic medium of infinite extent, which is bounded internally by a spherical cavity of radius *a*. The cavity is filled with an inviscid fluid identical to that saturating the pore space of the poroelastic medium. Poroelastic effects in the system can be induced by pressurization of the fluid cavity; thermo-poroelasticity effects can be induced by heating the fluid cavity. In each of these loading scenarios, the information of interest relates to the development and decay of fluid pressure *p* within the cavity, the circumferential stress *σ*_{ΦΦ}(*a*,*t*) and the radial displacement *u*(*a*,*t*) at the boundary of the cavity.

For spherically symmetric problems, the governing thermo-poroelasticity equations are
3.1where *M*=(*φ*/*K*_{f}+(*α*−*φ*)/*K*_{s})^{−1} is the Biot modulus and *φ* is the porosity.

Boundary and initial conditions are given as follows:
3.2The boundary conditions at infinity always correspond to the undrained state, which means that the applied stress , the fluid pressure and the volumetric strain at infinity must be related by
3.3where *K*_{u}=*K*_{D}+*α*^{2}*M* is the undrained bulk modulus.

The initial condition at *t*=0 corresponds to the *internal* fields, generated by the compressive stress suddenly applied at infinity at *t*=0. ‘Initial’ conditions at *t*=*t*_{0} for 0<*R*<*a* imply that sudden pressure and temperature increases are applied at time *t*=*t*_{0} to the fluid volume of the spherical cavity.

Let the initial fluid pressure in the cavity, generated by the applied stress , be denoted by as in (3.2). This pressure can be related to the applied stress if it is assumed that the initial state in the cavity is undrained. This condition gives , where *u*_{a} is the displacement of the cavity wall. As shown below (see (3.11)), the displacement of the cavity wall is related to the applied stress at infinity and the fluid pressure in the cavity by . From these relationships, we can obtain the initial fluid pressure in the cavity as
3.4Note that for very stiff geomaterials, the elastic moduli and , and thus the initial fluid pressure in the cavity reduces to zero.

The fluid flow equation in (3.1) can be simplified if we rewrite (3.1_{1}) as,
3.5Integrating this equation with respect to *R* and using the fact that at infinity the dependent variables are constant, we obtain
3.6Using this relationship, the fluid flow equation in (3.1) can be rewritten as
3.7The fluid flow equation now contains only two unknowns, the fluid pressure and the temperature, and thus becomes *decoupled* from the elasticity equation. It is convenient to introduce a new storage coefficient and a new thermal expansion of the solid phase, such that this equation becomes formally equivalent to the diffusion equation for the fluid pressure, used in ground water flow calculations [55,61]
3.8Equation (3.7), valid for the region of the geomaterial *R*≥*a*, must now be solved together with the heat conduction equation for the temperature field, i.e.
3.9To obtain an analytical solution of the problem in a closed form, it is also advantageous to replace conditions prescribed inside the spherical cavity with boundary conditions at the cavity wall and solve the problem for the region *R*≥*a* only. The boundary conditions at the cavity wall are
3.10where *T*_{A} is the volume average of the temperature in the cavity, *u* is the radial displacement. The first boundary condition is derived in appendix A. The displacement *u* can be eliminated from boundary condition (3.10_{1}) if we use the following relationships between the displacement at the cavity wall and the fluid pressure in the cavity
3.11The last relationship in (3.11) was already used for deriving (3.4). In turn, equation (3.11) can be derived using the traction boundary condition *σ*_{RR}=−*p* at the cavity wall and the relationship (3.6) between the volumetric strain and the fluid pressure. Namely, from traction boundary condition (3.10_{2}) and equation (3.6), we have a system of two linear equations for *u* and ∂*u*/∂*R*
3.12The first equation of (3.12) takes into account the fact that the pressure may be discontinuous at the cavity wall at time *t*=*t*_{0} when the sudden pressure increase (pulse) is applied, i.e. *p*(*R*=*a*+0)≠*p*(*R*=*a*−0). By solving (3.12) for *u* and ∂*u*/∂*R* we can derive (3.11).

Next, using the connection between the displacement of the cavity wall and fluid pressure (3.11), the boundary condition at the cavity wall (3.10_{1}) can be rewritten as
3.13Let us introduce a *modified* bulk modulus of the fluid in the cavity
3.14With this notation, the boundary condition at the cavity wall finally becomes
3.15In summary, we derived the equation for fluid pressure (3.9_{1}) decoupled from the elasticity equation and the corresponding boundary condition at cavity wall (3.15), also decoupled from the mechanical variables.

The solution of the differential equation for the fluid pressure and temperature can be obtained by using a Laplace transform technique. First, we derive a solution for *t*<*t*_{0} when the stress is applied at infinity. For this case, *T*=0. Applying a Laplace transform to differential equation (3.9_{1}) gives
3.16where *s* is the Laplace transform variable. Applying a Laplace transform to boundary condition (3.15) at the cavity wall gives
3.17The solution of (3.16) for the pressure field can be represented by
3.18where *c*^{2}=(*η*/*K*)*sS*. The constant *A* can be found from boundary condition (3.17). It gives
3.19The resulting fluid pressure at the cavity wall becomes
3.20The initial fluid pressure in the cavity is given by (3.4).

Now consider the solution at *t*>*t*_{0} after a sudden increase in the fluid pressure (pulse) and/or the temperature. Owing to the linearity of the problem, we can assume that the applied compressive stress at infinity is zero, and then use the principle of superposition to take into account the solution at *t*<*t*_{0}, which is caused by the compressive stress. In addition, to simplify the derivation, we can set *t*_{0}=0.

The sudden change in the fluid pressure *p*_{0} can be caused by rapid fluid injection, leakage and/or a rapid increase in temperature *T*_{0}. Consider first the case when *T*=0; here, the Laplace transform of the fluid pressure can be obtained in a similar manner to the case of the applied compressive load (see (3.18) and (3.19)). The Laplace transform of the pressure is given by
3.21Consider now the case where a rapid increase in temperature *T*_{0}≠0 is applied to the volume of the cavity. Let us first obtain the solution to the heat conduction equation (3.9_{2}). In order to simplify the analytic derivation, we assume that the thermal properties of the inclusion and the surrounding geomaterial are identical. The Laplace transform version of heat conduction equation (3.9_{2}) is given by
3.22The solution can be represented by
3.23where *γ*^{2}=*sρc*_{p}/*k*_{c}. Using the standard temperature and heat flux continuity conditions at *R*=*a*, we find that
3.24The volume average of the temperature field in the cavity can be found by integration of (3.24). This gives
3.25The moment the temperature increase is applied, there is no change in the fluid content inside the cavity, i.e. the cavity is ‘initially’ undrained. The initial pressure *p*_{0} in the cavity is
3.26Again, using the connection between the displacement of the cavity wall *u*_{a} and the fluid pressure given by (3.11), we can show that
3.27The pressure field in the geomaterial, as a function of time, can be found by solving differential equation (3.9_{1}). Recall also the boundary condition at cavity wall (3.15); applying the Laplace transform to differential equation (3.9_{1}) gives
3.28and applying the Laplace transform to the boundary condition at the cavity wall (3.15) yields
3.29owing to the non-zero initial values of the pressure and temperature.

The solution of (3.28) can be represented by
3.30where, as before, *c*^{2}=(*η*/*K*)*sS*. The constant can be found by substituting this solution into differential equation (3.28) and using (3.24), resulting in
3.31The remaining constant *A* is found by enforcing boundary condition (3.29) giving
3.32With the use of (3.31) and (3.32), the result for the transformed expression for the pressure (3.30) at the cavity boundary, *R*=*a*, is simplified to
3.33Note that this expression is valid for all points within the cavity, *R*≤*a*, because the pressure within the cavity is *uniform*.

The inverse Laplace transform for the pressure (3.33) and temperature (3.24) can be found numerically by using an expansion involving Legendre polynomials [21]. Alternatively, the pressure and temperature fields can be found by evaluation of the Bromwich integral, the details of which are presented in Rice *et al*. [26]. Using the latter approach, the temperature at the cavity wall and at the centre of the cavity can be expressed as
3.34Note that the temperature field is *not uniform* within the cavity.

The resulting fluid pressure in the cavity can be expressed as
3.35where
3.36The first contribution to the pressure field is made by the initial pressure increase *p*_{0}. The second part is due to the thermal expansion of the fluid *inside* the cavity, whereas the third part is due to the thermal expansion of the geomaterial constituents. For convenience, a further change of variables *y*=*s*^{2} can be made in (3.36).

It is worth noting that the denominator [(*s*−*D*)^{2}+*D*^{2}*c*^{2}*a*^{2}] has complex roots if 3K_{f}*S*<4. This constraint is satisfied by most geomaterials, particularly commonly occurring rocks. Therefore, the integrands have no singularities on the interval of real numbers and the integrals can be evaluated numerically. In the case of *real* roots of the denominator, the integral in *p*_{M}(*R*,*t*) can also be evaluated analytically [62].

Now, we can evaluate the stresses at the cavity wall. The effective radial stress can be directly evaluated using the traction boundary condition at the cavity wall (3.10_{2}) and the definition of the effective stress, i.e.
3.37where *p*_{−}=*p*(*a*−0) and *p*_{+}=*p*(*a*+0), and is the effective stress. Here, we take into account the fact that, at the instant when the pressure pulse or sudden temperature change is applied, there is a pressure jump at the interface. For all other times, the pressure is continuous, i.e. *p*(*a*)=*p*(*a*−0)=*p*(*a*+0).

Using the expression for the radial stress at the cavity wall (3.12_{1}), the expression for the volumetric strain (3.6) and the stress–strain relationship for the effective hoop stress , the hoop stress at the cavity wall can be expressed as,
3.38or, alternatively,
3.39where *T*_{+}=*T*(*a*+0). The effective pressure at the cavity wall can now be evaluated from (3.37) and (3.38) as
3.40

## 4. Results for the thermo-poroelasticity problem

In this section, we present results pertaining to the thermo-poroelastic problem. Consider a fluid-filled cavity of radius *a*=0.02 m embedded in a geomaterial of infinite size. The porosity of the geomaterial is *φ*=0.01, and the Biot coefficient is *α*=0.7. The elastic properties of the geomaterial are bulk modulus *K*_{D}=5 GPa, and Poisson's ratio *ν*=0.3. The permeability of the geomaterial is *K*=3×10^{−19} m^{2}, and the dynamic viscosity of the pore fluid is *η*=0.001 Pa s. The bulk modulus of the fluid is *K*_{f}=2.2 GPa. The thermal properties are thermal expansion of the fluid *α*_{f}=69×10^{−6} 1/C and the thermal expansion of the solid phase *α*_{s}=8.3×10^{−6} 1/C. The thermal conductivity of both the fluid cavity and the geomaterial is *k*_{c}=3.15 Wm^{−1} C^{−1}, whereas the specific heat multiplied by the density is *ρc*_{p}=1912913 J m^{−3} C^{−1}.

As the elastic problem is linear, it makes sense to separate the effects caused by the compressive stress applied at infinity, the sudden change in the fluid pressure and/or the temperature applied to the volume of the cavity. Each of these loadings will thus be considered separately. The computational modelling was performed using the COMSOL Multiphysics code. The domain used in the finite-element modelling of the infinite region was a two-dimensional sector bounded externally at *R*=20*a*, where *a* is the radius of the cavity, and contained by planes *Θ*=0 and *Θ*=*π*/2. The element used in the finite-element modelling was of the six-node Lagrangian type; the displacement, the fluid velocity and the heat flux were all prescribed to be zero in the normal direction, on all nodal points of the bounding planes *Θ*=0 and *Θ*=*π*/2. The boundary conditions on the outer region of the finite-element domain corresponded to the applied stress at infinity and zero fluid velocity in the radial direction to represent undrained conditions. The computational results obtained for the thermo-poroelasticity problems, from the computational code COMSOL, are indicated in figures 1–4 and they confirm the validity of the modelling of infinite domains with appropriate choices for the positioning of the far-field boundaries and boundary conditions. The results presented in this paper generally support the conclusions arrived at by the authors in their previous investigations [21].

Figure 1 shows the evolution of the fluid pressure in the cavity when a compressive radial stress of magnitude 1 MPa is applied to the geomaterial as a step function at infinity. The corresponding fluid pressure in the geomaterial at infinity is evaluated from (3.3) as
4.1assuming that the geomaterial is undrained at infinity at all times. This gives us MPa. The initial fluid pressure in the cavity is evaluated from (3.4) as
4.2assuming that the fluid cavity is initially undrained. The result is . Figure 1 shows the fluid pressure obtained using both the coupled theory of poroelasticity, in which the hydro-mechanical (HM) interaction between the unknown variables is taken into account, and the uncoupled theory. Using the term ‘uncoupled solution’, we imply that the solution is obtained by neglecting both the volumetric strain *ε*_{V} in the fluid pressure equation and the displacement of the cavity wall *u*_{a} in boundary condition (3.10) at the cavity wall. The initial conditions for the coupled and uncoupled problems are assumed to be the same. We can see that the fluid pressure obtained by solving the coupled system of equations of poroelasticity shows a slower change. In general, the rate of change in the fluid pressure depends not only on the permeability but also on the ratio *K*_{f}/*K*_{D}. When this ratio decreases (*K*_{D} increases), the rate of change in the fluid pressure for the coupled solution will be more rapid, and the difference between the two solutions, coupled and uncoupled, will become less notable.

Consider now the effect of a sudden change in the fluid pressure inside the cavity. Figure 2 shows the evolution of the fluid pressure in the cavity when there is a sudden increase in the fluid pressure (pressure pulse) of magnitude *p*_{0}=1 MPa applied to the volume of the cavity at *t*=0. This increase in the fluid pressure can be caused, for example, by fluid injection or a sudden temperature change if there is a large thermal resistance that does not permit heat flow across the cavity wall. The applied compressive stress at infinity is set to zero. Again, we present two solutions corresponding to the coupled and uncoupled system of equations. We observe that the rate of change of the fluid pressure for the coupled solution is slower than that for the uncoupled solution.

Figure 3 shows the effect of the permeability of the geomaterial on the evolution of the fluid pressure in the cavity. Three values for the permeability were used: the reference value *K*_{0}=3×10^{−19} m^{2}, and values three times smaller and larger than the reference value. Again, the sudden increase in the fluid pressure (pressure pulse) of magnitude *p*_{0}=1 MPa is applied to the volume of the cavity. Figure 3 presents the solution obtained by solving the coupled system of equations of poroelasticity. We can see that the rate of change in the fluid pressure in the cavity depends on the permeability of the geomaterial, and thus the application of the pressure pulse can potentially be used in experiments to measure the permeability of geomaterials.

Figure 4 illustrates the effect of a sudden (or initial) temperature change *T*_{0} in the volume of the fluid cavity. The applied compressive stress at infinity is assumed to be zero, and the thermal resistance at the cavity wall is absent. The initial temperature increase within the cavity is taken as *T*_{0}=11^{°}C. It is assumed that the cavity is initially undrained, i.e. no injection or leakage of the fluid occurs initially in the cavity, unlike the previous case. This condition allows us to obtain the initial pressure in the cavity as (see (3.27))
4.3This result gives *p*_{0}=2.92 MPa. From figure 4, we can see that the dissipation of the fluid pressure is faster compared with the case of the ‘pure’ pressure pulse, i.e. without any temperature change. This can be explained by the fact that the dissipation of the initial fluid pressure (4.3) is caused not only by the fluid flow across the cavity wall, but also by the heat dissipation. We can see that the fluid pressure can become negative. In figure 4, we present the solutions to both the thermo-hydro-mechanically coupled and uncoupled systems of equations. In the latter case, the effect of the volumetric strain of the porous skeleton *ε*_{V} in fluid flow equation (2.7) and the skeletal displacement of the cavity wall *u*_{a} in boundary condition (3.10) are neglected. The initial pressure in the cavity for both solutions is assumed to be the same, i.e. *p*_{0}=2.92 MPa. We can see that in the case of the coupled solution, the fluid pressure dissipation is much slower, and the magnitude of the negative pressure is small. It should be noted that the fluid pressure in the cavity could be generated not only by heating the fluid within the cavity but also by *injecting* a preheated fluid into an empty cavity. In this case, the volume of fluid that must be injected to create the same pressure as (4.3) is derived in appendix B.

Figure 5 shows the influence of the permeability of the geomaterial on the evolution of the fluid pressure in the cavity subjected to a sudden increase in the temperature *T*_{0}=11^{°}C. As before, three values of the permeability of the geomaterial are used: *K*=3×10^{−19}, *K*=1×10^{−19} and *K*=9×10^{−19} m^{2}. Owing to the relatively rapid dissipation of the fluid pressure in the cavity, it is more difficult to distinguish among the three solutions and to use the sudden temperature increase in an experiment for measuring the permeability of a geomaterial. However, we can observe that the time required for the pressure to intersect the horizontal axis (zero pressure level) clearly depends on the permeability of the geomaterial and this fact can potentially be useful in experiments.

The fluid pressures shown in figures 1–5 have been obtained by inversion of transformed solution (3.33) using the numerical inversion algorithm [21]. Alternatively, the integrals in (3.35), (3.36) can be evaluated numerically. It was confirmed that these solutions compare very accurately with results derived using the finite-element method. In the finite-element method, boundary condition (3.10) at the cavity wall can be difficult to implement; instead, the fluid cavity is included in the model and represented as a porous material for which the permeability *K* is very large, whereas the elastic stiffness *K*_{D} is very small. Both the porosity and the Biot coefficient for such a porous material must be set to unity. The results of finite-element modelling obtained using the COMSOL software are indicated in dotted lines in figures 1–4.

## 5. Results for the thermo-poroelasto-plasticity problem

In this section, we present the computational results pertaining to the thermo-poroelasto-plastic problem. The elastic and thermal parameters for the geomaterial and fluid cavity are the same as given in §4. The initial yield stress or initial size of the yield surface is . Further evolution of the yield surface is controlled by the hardening rule , and the parameter *M*_{y} in (2.9) is set to unity.

The geomaterial region containing the fluid-filled cavity is first subjected to a far-field compressive radial stress applied at infinity. Once a steady state is reached (0.9726 MPa), after approximately 2000 s, a sudden increase in the pressure or temperature is applied to the fluid cavity. This increase in the fluid pressure can be caused by fluid injection and/or heating.

Figure 6 shows the fluid pressure evolution in the cavity when, at time *t*=0, the geomaterial is subjected to the compressive radial stress and, at time *t*_{0}=2000 s, a sudden increase in the fluid pressure is applied to the volume of the cavity. Solutions are given for both poroelastic and poroelasto-plastic geomaterials. For *t*<2000 s, these solutions coincide, i.e. the applied compressive radial stress does not cause any plastic strains. When the pressure increase (pulse) is applied to the volume of the cavity, plastic strains occur. The maximum plastic strain is reached at the time of application of the pressure pulse, i.e. at *t*=2000 s; the plastic strains remain constant after *t*=2000 s. If the same amount of fluid is injected to the volume of the cavity, the instantaneous increase in the fluid pressure at *t*=2000 s will be smaller for the poroelasto-plastic geomaterial. Therefore, to achieve the same value of the instantaneous pressure increase in the cavity, more fluid must be injected into the spherical cavity contained within a poroelasto-plastic geomaterial. The maximum value of the fluid pressure is chosen to be 21.9 MPa, both for the poroelastic and poroelasto-plastic geomaterials. It is evident from figure 6 that, although the maximum fluid pressure values are equal for the poroelastic and poroelasto-plastic geomaterials, the rate of fluid pressure dissipation is greater when the geomaterial is poroelasto-plastic.

Figure 7 illustrates the fluid pressure evolution in a cavity subjected to a temperature increase of 42^{°}C at time *t*_{0}=2000 s. As before, the geomaterial is subjected to a far-field compressive radial stress applied at infinity at *t*=0 s. Solutions are given for both poroelastic and poroelasto-plastic geomaterials. A corresponding increase in the fluid pressure at time *t*_{0}=2000 s is observed assuming that the cavity is undrained. For *T*=42^{°}C, the instantaneous increase in the pressure is equal to 11.1527 MPa if the geomaterial is elastic. We note that this instantaneous increase in the fluid pressure must be smaller for the poroelasto-plastic geomaterial, but this is not apparent in figure 7 due to the fact that the plastic strains are small for the selected temperature increase. If the magnitude of the applied temperature is chosen to be larger, then the difference between the instantaneous increases in the fluid pressure for the poroelastic and poroelasto-plastic geomaterials will become more evident. As can be seen in figure 7, the decrease in the fluid pressure is faster for the poroelasto-plastic geomaterial, and the extremum point in the fluid pressure evolution curve is more apparent.

Figure 8 shows the radial displacement of the surface of a cavity located either in a poroelastic or poroelasto-plastic geomaterial region. The geomaterial is first subjected to a compressive radial stress applied at infinity and, after reaching a steady state, a pressure pulse is applied within the cavity. The maximum resulting fluid pressure in the cavity immediately after application of the pressure pulse is 21.9 MPa (the same conditions that were used in figure 6). Figure 8 demonstrates that when the geomaterial is subjected to the compressive stress applied remotely, the radial displacement of the cavity is small but negative (contraction). Fluid flows towards the cavity centre as the fluid pressure in the cavity is initially smaller than that in the geomaterial (figure 1). After application of the pressure pulse, the cavity expands to accommodate the fluid that was injected. As the injected fluid flows from the cavity into the geomaterial, the radial displacement of the cavity gradually decreases. If the geomaterial is poroelastic, all the fluid that was injected is transferred to the geomaterial in the long term. In the case of the poroelasto-plastic geomaterial, however, there is a residual cavity displacement, which suggests that not all the fluid that was injected initially is able to flow into the geomaterial.

Finally, it should be mentioned that the poroelasto-plastic solutions were obtained using the finite-element code ABAQUS which is well suited for solving the poroelasto-plasticity problems. As shown before, the fluid cavity was modelled as a poroelastic medium with porosity equal to unity, with a comparatively large permeability and a relatively small value for the elastic stiffness. The heat conduction problem was solved first; the resulting temperature field was used as an input to the second problem, for which the unknowns were the displacements and the fluid pressure. Similar to the computational modelling of the poroelastic THM problem, the domain used in the poroelasto-plastic finite-element modelling of the infinite region was a sector bounded externally at *R*=20*a*, where *a* is the radius of the cavity, and contained by the planes *Θ*=0 and *Θ*=*π*/2. The element used in the finite-element modelling was of the six-node Lagrangian type; zero displacement, zero fluid velocity and zero heat flux were prescribed in the normal direction on all nodal points of the planes *Θ*=0 and *Θ*=*π*/2. The boundary conditions on the outer boundary of the finite-element domain corresponded to the applied stress at infinity and zero fluid velocity in the radial direction.

## 6. Concluding remarks

The coupling of thermo-poroelasto-plastic effects in fluid-saturated porous media continues to be of interest to many areas in the geosciences and geomechanics. This paper presented a benchmark problem that gives a comparison of results for the transient response of poroelastic and poroelasto-plastic geomaterial domains of infinite extent that are bounded internally by a fluid-filled spherical cavity. The analytical results were derived using a Laplace transform technique and analytical results are used to calibrate the results obtained from the two finite-element codes COMSOL and ABAQUS. The computational results compare favourably with the analytical solutions.

Analytical solutions are obtained for several loading cases, namely a far-field compressive radial stress, a sudden pressure change (pulse) applied to the volume of cavity and a sudden temperature change applied to the cavity. The evolution of the fluid pressure generated in the cavity, or equivalently at the cavity/geomaterial interface, was studied for each of these loading cases.

Application of a far-field compressive stress generates a positive pressure inside the cavity of a magnitude that is initially smaller than the far-field value of the fluid pressure. The fluid pressure in the cavity then increases with time until it reaches a steady-state value equal to the far-field fluid pressure. In addition, we note that the circumferential stress created at the cavity wall for this type of loading is compressive and it reduces or eliminates the plastic strains in the geomaterial if the pressure pulse and/or temperature increase are applied afterwards.

The pressure pulse in the cavity can be generated by a sudden fluid injection or leakage, and also by a sudden temperature change if there is a high thermal resistance at the cavity/ geomaterial interface. As the results indicate, the fluid pressure in the cavity when subjected to the pressure pulse decreases with time. The rate of the pressure dissipation clearly depends on the permeability of the geomaterial, which allows us to use these results in experiments for measuring permeability. The displacement of the cavity wall, given in terms of the fluid pressure, *u*_{a}=*pa*/(4*G*_{D}), enables us to find the volume of the fluid that is currently residing in the cavity. The results presented in the paper also point out the distinctions that arise in the responses when full THM coupling is incorporated.

The fluid pressure caused by a sudden temperature increase in the cavity is initially positive, but as time progresses the fluid pressure reaches zero and becomes negative until it returns to a steady-state zero value. The time required for the pressure to first reach zero depends on the permeability of the geomaterial and does not depend on the value of the applied temperature, if the problem is linear. This fact can also be used in experiments for measuring permeability of geomaterials. In addition, the expression for the value of the initial pressure in the heated cavity, derived under the assumption that the cavity is initially undrained, can be used to infer the value of the shear modulus of the geomaterial (see (3.27)). It is also important to note that, in the case of temperature increase, the rate of change in the pressure is faster compared with the case of a ‘pure’ pressure pulse, owing to heat dissipation into the geomaterial. The circumferential stress at the cavity wall generated by the ‘pure’ pressure pulse and/or by a temperature increase is tensile, which induces plastic strains in the poroelasto-plastic geomaterial.

For poroelasto-plastic geomaterials, the loading history consisted of the applied compressive radial stress at infinity and, after 2000 s when a steady state was reached, the application of a pressure pulse or sudden temperature increase in the cavity. The applied stress at infinity was chosen to be sufficiently small that it did not induce any plastic strains, but the subsequent application of a pressure pulse or temperature increase generated plastic strains. The magnitude of these plastic strains was, however, reduced because of the compressive radial stress applied at infinity. This can be explained by the fact that the hoop stresses at the cavity wall generated by the far-field compressive radial stress and by the pressure pulse or temperature increase have different signs. A comparison of the fluid pressure evolution in the cavity for the poroelastic and poroelasto-plastic geomaterials suggests that, in general, the fluid pressure in the elasto-plastic geomaterial is smaller and the rate of change in the fluid pressure is faster. The results presented are relevant to geomaterials where the skeletal elasto-plastic behaviour can be modelled by the modified Cam-Clay plasticity model. The procedures can be extended to include other plasticity models appropriate for soft geomaterials and heavily overconsolidated argillaceous geomaterials.

## Disclaimer

The use of the computational codes COMSOL and ABAQUS is purely for demonstration purposes only. The authors do not advocate or recommend the use of these codes without conducting suitable validation procedures that test their accuracy in a rigorous fashion.

## Funding statement

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 kind remarks and constructive suggestions of the reviewers are gratefully acknowledged.

## Appendix A

Here, we derive the boundary condition (3.10_{1}) at the cavity wall *R*=*a*. It can be derived using the divergence theorem; however, it is more instructive to derive it directly from the system of differential equations for a porous material applied to the fluid cavity. In fact, the fluid cavity can be treated as a porous material with porosity *φ*=1, Biot coefficient *α*=1, the small elastic moduli of the skeleton *K*_{D}, *G*_{D}→0, large permeability and zero thermal expansion of the solid phase *α*_{s}. For this special case of a porous material, *M*^{−1}=1/*K*_{f}.

Differential equations (3.1) applied to the case of such a cavity become
A 1where *v*_{f} is the radial fluid velocity in the cavity. As the elastic moduli are small, the equilibrium equation implies that *p*_{,R} is equal to zero, i.e. the fluid pressure *p* is uniform, and thus only a function of time. Assume that the fluid velocity and the displacement in the cavity are linear functions of the radius *R*, i.e.
A 2where *v*_{fa} is the fluid velocity at *R*=*a* and *u*_{a} is the displacement of the cavity wall, *R*=*a*.

Then, the volumetric strain in the cavity becomes
A 3and is, therefore, uniform. Substituting expressions for the fluid velocity (A 2) and the volumetric strain (A 3) into (A 1) gives
A 4The fluid velocity at the cavity wall is given by
A 5Thus,
A 6We note that this equation cannot be satisfied exactly because, on the left-hand side, we only have a function of time, but on the right-hand side of the equation we have a function of time and radial position, i.e. *T*(*R*,*t*). However, we can satisfy this equation in a weak sense, i.e. considering the average of the temperatures in the cavity, i.e.
A 7where
A 8This completes the derivation of boundary condition (3.10_{1}). Boundary condition (A 6), but without a temperature contribution, also appears in De Josselin de Jong [62].

## Appendix B

Instead of heating the fluid already resident in the cavity, we can initially pressurize the cavity by injecting preheated fluid into the empty cavity. If the temperature rise is equal to *T*_{0}, the volume of the heated fluid that must be injected can be evaluated as *V*^{inject}_{f}=*V*_{0}(1+3*α*_{f}*T*_{0}), where *V*_{0}=4/3*πa*^{3} is the volume of the fluid at the reference temperature state, equal to the cavity volume. From (4.3), it follows that the resulting pressure in the cavity will be equal to
B 1Eliminating the thermal expansion coefficient from these two expressions, we obtain
B 2Therefore, the volume of the fluid that must be injected into the *empty* cavity is
B 3

- Received September 24, 2013.
- Accepted November 26, 2013.

- © 2014 The Author(s) Published by the Royal Society. All rights reserved.