Consolidation around a point heat source (correction and verification)

We provide a brief review of the results obtained by Booker and Savvidou in their original 1985 paper describing the consolidation around a point heat source embedded in a fluid‐saturated porous medium. Although it gained a lot of attention in recent years, mainly as a benchmark problem for numerical codes, the analytical expressions for normal and shear stresses could not be confirmed to be correct. In this manuscript, we verify the expressions for temperature, pressure, and displacement and derive new expressions for the stress tensor withstanding numerical verification using the scientific open‐source thermo‐hydro‐mechanical‐chemical (THMC) simulator OpenGeoSys.

solid skeleton u i . A rising temperature causes a distinct response of the fluid and solid skeleton. The resulting pressure gradient causes the fluid to flow, reaching equilibrium as the time passes and pore pressure dissipates. The consolidation depends mainly on the interplay between hydraulic and mechanical properties of the porous medium. One example for the use of this solution in code verification of different modeling teams is Task E in the Decovalex 2019 phase, * which is aimed at simulating THM processes in Callovo-Oxfordian clay rock. 10 Although the mentioned analytical solution has been extensively used in the recent past (193 citations † at the time of writing this article) and serves as a benchmark to validate numerical solvers, there are a number of inaccuracies, notational inconsistencies, and, in case of the stress solution, incorrect results. Given the usefulness of this paper in recent research, this article is an attempt to identify and correct these inaccuracies and provide a clean solution that is also verified numerically using the open-source THMC simulation platform OpenGeoSys. 8,11 In Section 2, we give the general problem statement and list the system of equations governing a typical THM problem. In Section 3, the complete analytical solution is listed. We divide the solution into two parts. The first part, which is similar to the original work, is referred to as the verified part, and the second part contains the corrected solution for stresses and, thus, is referred to as the corrected part. All results are presented in a Cartesian coordinate system. Furthermore, in Section 3.3, we also list the set of equations in spherical coordinates. Section 4 contains the numerical verification of the solution using OpenGeoSys. Throughout the paper, we stick to the nomenclature from the original publication unless consistency demands otherwise.

GENERAL PROBLEM STATEMENT WITH GOVERNING EQUATIONS
It appears that the sign convention used in the original work is not consistent. In this work, we use the solid mechanics sign convention (tensile stress being positive). We use comma convention to indicate the partial derivative, i.e, (·) , i = (·)∕ x i . A repeated index in case of subscripts i, j, k, l will always imply summation (1 to 3) unless explicitly stated otherwise.
As the name indicates, a coupled THM model is based on a thermal part, a hydraulic part, and a mechanical part, which are not independent of each other. The thermal part in this case is based on the energy balance equation, which is given as 12 m where q T is a heat source per unit volume, and m, K, and v i (Darcy velocity) are given as The hydraulic part is based on the mass balance equation, which is given as where q H is the source term for the fluid, while a u is given as The mechanical part is based on the balance of momentum of the mixture and is given as where = w + (1 − ) s and ij is the total stress, which is given as where ij refers to Kronecker delta and ′ i is the effective stress tensor, which is given as where C ijkl and kl are the stiffness and strain tensors, respectively. For isotropic case, the above equation can be rewritten as

THM CONSOLIDATION AROUND A POINT HEAT SOURCE IN AN INFINITE MEDIUM-ANALYTICAL SOLUTION
Before proceeding further, it is appropriate to mention the assumptions made while deriving the analytical solution in the original work. No advection is considered in the energy balance, so the second term on the left-hand side of Equation (1) vanishes. The energy balance thus becomes uni-laterally coupled to the HM part of the model. Also, the gravitational force is neglected, so the last term on the right-hand side of Equation (4) and the left-hand side of Equation (7) vanishes. The solid and fluid phases are both assumed to be intrinsically incompressible, which results in B = 1 and = 0. As a consequence, the volume change of the phases themselves is caused only by the change in temperature. No external fluid source or sink is present; ie, q H = 0. The following is the list of derived parameters used in the original work: It should be noted here that the symbol was wrongly used in the expression for c (coefficient of consolidation) in the original work. Furthermore, we avoid the use of expression f * and derive the relations for stress from the expression for displacement, which was found to be correct in the original work in combination with the constitutive law for thermo-elasticity. For this purpose, we need the following expressions:

Verified solution
The solutions for T, p, and u i are found to be correct in the original work and therefore listed here without any change.

Corrected solution
The corrected solution for stresses is derived by hand using Equation (10) and verified numerically in Section 4. The following are the expressions for the normal effective stresses and the shear effective stresses.
The total stresses follow from Equation (8).

Spherical coordinates
Since the THM consolidation around a point heat source represents a spherically symmetric problem, it is convenient to express the corresponding analytical solution in spherical coordinates. In this work, we used the displacement solution in Cartesian coordinates to obtain the stress solution as given below. It should be noted here that the same results should be obtained by using x 1 = r, x 2 = x 3 = 0 or any other equivalent expression in the displacement and stress solutions in  Cartesian coordinates given in the previous section.

Computational details
The verification of the analytical model was performed using version 6.2.0-44 of the scientific open-source THMC simulator OpenGeoSys ‡ on a two-dimensional 20 m × 20 m model region with rotational axial symmetry representing a half-space configuration of the spherically symmetrical problem. The modeling domain was subdivided into 3830 quadrilateral mesh elements with an element size growing from 10 −3 to 10 0 m. Along the inner (symmetry) boundaries, we constrain the axis-normal displacements to zero, while we require to have no effects on pore pressure (p = 0) and temperature (T = T 0 = 273.15 K) along the outer boundaries. Starting from an initial pressure of p = 0 and T 0 = 273.15 K, we propagated the heat coming from a point heat source of Q = 300 W located at the origin through the medium with THM coupling as implemented in the THERMO_HYDRO_MECHANICS process in OpenGeoSys using a fixed time-stepping scheme with Δt = 5000 seconds. This process solves the equations listed in Section 2 by choosing pressure, temperature, and displacements as primary variables. Isoparametric quadratic elements were used for the displacement field, while linear elements were used for temperature and pressure. Nonlinearities are resolved with an iterative-incremental Newton-Raphson solution strategy, while time integration is fully implicit based on a backward Euler scheme. The tolerances for the nonlinear solver were set to 10 −10 for temperature and displacement components and 10 −5 for the pressure.
The linear system was solved using the direct SparseLU solver from the Eigen library. §

Numerical verification of the analytical solution
To provide further support for the verified and corrected analytical solutions (Section 3), we evaluated the expressions together with the numerical model as function of time at an arbitrarily chosen point P = (1.30m, 0.682m, 0.0) as well as at t = 5 × 10 5 seconds as a function of r for a given parameter set ( Table 1). The results are presented in Figure 1 for the temporal and in Figure 2 for the spatial dependencies of (a) temperatures, (b) pressures, (c) displacements, and (d) stresses, respectively. We find for the verified expressions of Booker and Savvidou 9 as well as for the stress expressions derived in Section 3.2 of this paper an excellent agreement between the analytical and numerical models. Greater deviations from the analytical model can only be found for the nondiagonal stress elements where r approaches zero and the solutions for temperature and normal stresses become divergent ( Figure 2).

CONCLUSION
We revisited the analytical solution presented by Booker and Savvidou 9 in their original publication on a consolidation problem around a point heat source. In addition to the resolution of several minor inaccuracies, notation errors, and conventional inconsistencies, a correction to the solution of the stress tensor was presented. The corrected solution provided here was verified numerically using finite element analysis, and it appears to be in agreement with the numerical results. This comparison was introduced into the automatic test suite (benchmarks) of the OpenGeoSys project and is available with the source code.