## Abstract

The Green’s function of a bimaterial infinite domain with a plane interface is applied to thermal analysis of a spherical underground heat storage tank. The heat transfer from a spherical source is derived from the integral of the Green’s function over the spherical domain. Because the thermal conductivity of the tank is generally different from soil, the Eshelby’s equivalent inclusion method (EIM) is used to simulate the thermal conductivity mismatch of the tank from the soil. For simplicity, the ground with an approximately uniform temperature on the surface is simulated by a bimaterial infinite domain, which is perfectly conductive above the ground. The heat conduction in the ground is investigated for two scenarios: First, a steady-state uniform heat flux from surface into the ground is considered, and the heat flux is disturbed by the existence of the tank due to the conductivity mismatch. A prescribed temperature gradient, or an eigen-temperature gradient, is introduced to investigate the local temperature field in the neighborhood of the tank. Second, when a temperature difference exists between the water in the tank and soil, the heat transfer between the tank and soil depends on the tank size, conductivity, and temperature difference, which provide a guideline for heat exchange design for the tank size. The modeling framework can be extended to two-dimensional cases, periodic, or transient heat transfer problems for geothermal well operations. The corresponding Green’s functions are provided for those applications.

## 1 Introduction

Thermal management of building environment through heat exchange with the ground has attracted significant attention and applications for energy efficiency improvement and environmental benefits [1,2]. A water tank is generally employed to store energy from building roof [3] and to exchange heat with the ground, which keeps at a stable temperature in general. By using the thermal mass of the earth, we can manage the indoor temperature of buildings as if they were embedded in the earth, like a basement. Figure 1 shows a schematic illustration of a bidirectional geothermal system [3]. The water tank serves as both a thermal energy storage (TES) and a heat exchanger with the surrounding earth. A heat pump can regulate the heating or cooling supplies to the building with a high coefficient of performance (COP) for the ground source heat pump (GSHP) system.

The performance of the bidirectional geothermal system relies on the heat transfer capacity through the surface of the water tank. It is of great significance to understand the heat conduction and temperature distribution of an underground water tank under the thermal loading for the system design and performance analysis. Particularly, given the temperature of heat source, the heat transfer capacity from the storage tank to the earth should be understood in the design phase, so that the size and the location of the storage tank can be determined to maximize the energy harvesting and energy efficiency.

A single storage tank underground has been modeled by an inhomogeneity embedded in an infinite domain [3,4], and the effect of the ground surface has not been taken into account. Particularly, the ground temperature profile exhibits variation with the depth and time due to the surface temperature change [5]. It is necessary to consider the effect of the ground surface on the temperature distribution and heat conduction in the neighborhood of the tank, particularly when the diameter of the tank is not too small in comparison with the depth.

Green’s function has been utilized as a powerful tool in solving various partial/ordinary differential equations and is also the foundation of the boundary integral equation method [6]. Serving as the fundamental solution to these problems with a generalized mathematical form, Green’s functions have been applied to many engineering fields including civil, mechanical [7], electrical, materials engineering [8], and applied physics and science [9]. Particularly, Green’s functions have been used extensively in the thermal analysis considering different types of heat sources, material configurations, and boundary conditions [10].

Green’s function is typically applicable to a homogeneous material with heat sources. For inhomogeneities embedded in a matrix, the material mismatch may lead to disturbance of the local field. Eshelby’s equivalent inclusion method (EIM) provides a unique way to simulate the material mismatch by an eigen-field [11–13]. Although it was originally established for elastic problems, the EIM has been applied to many other physical problems and the cases of multiple inhomogeneities with particle interactions. Some representative applications include the prediction of elastic moduli [12], thermal expansion coefficients [14–16], and steady-state thermal conductivity [17,18].

In the literature, inhomogeneity problems in steady-state heat conduction have been successfully solved by standard analytical methods [19], semi-analytical methods [20,21], and numerical methods. Later Hatta and Taya [17] unified these methods to solve the inhomogeneity problems of steady-state heat conduction by analogy with the EIM in elasticity. He also extended the EIM to the case of multiple ellipsoidal inhomogeneities in an infinite isotropic matrix. Yin et al. [18] used the EIM to consider the interfacial thermal resistance of spherical inhomogeneities. Because the integral of the Green’s function over a spherical or ellipsoidal inclusion can be analytically derived, the explicit solution can be obtained. By using numerical integrals, EIM can be extended to transient heat conduction problems as well [22].

For inhomogeneities in a semi-infinite domain or a bimaterial infinite domain, the thermal analysis has not been well studied because the boundary effect will play a role when the inhomogeneities are close to the surface or the interface of the bimaterial. Actually, the semi-infinite case can be considered as a bimaterial infinite domain with one material being perfectly conductive or insulative for uniform temperature and zero heat flux boundary conditions, respectively [23]. The material mismatch between inhomogeneity and matrix as well as the bimaterial will induce the discontinuity and singularity of thermal fields. The elastic problems for inhomogeneities in a semi-infinite domain have been solved for three-dimensional (3D) [24] and two-dimensional (2D) cases [25], respectively.

This article aims to solve the boundary value problem [26] for an inhomogeneity in a semi-infinite domain and apply it to the geothermal system with a spherical heat storage tank. The Green’s function derived in this article is applied to formulate one inhomogeneity inside a semi-infinite matrix with nonuniform “eigen-temperature gradient” (ETG) [18]. For an ellipsoidal inhomogeneity in an infinite domain, the ETG should be uniform [11,12], and classic micromechanics-based models were built upon it [27–30]. When the boundary effect or particle interactions are considered, the eigen-strain for elastic problems or ETG for thermal problems will not be uniform any more and can be represented by the Taylor expansion on each inhomogeneity [28]. Initially, the Green’s function for a bimaterial infinite domain will be constructed. The inhomogeneity problem is then solved with the EIM and verified with the finite element method. Although the present analysis focuses on the steady-state heat transfer of a spherical underground tank in a semi-infinite domain, it can be extended to cases with other shape of tanks [31,32], periodic temperature change condition, and transient heat transfer problems [23]. The Green’s functions of a bimaterial for 2D infinite domain, sinusoidal temperature field, and transient heat transfer problems are derived for future applications of transient heat conduction.

## 2 Formulation of Green’s Function

*D*composed of upper

*D*

^{+}and lower

*D*

^{−}parts, divided by the interface $S$ of

*x*

_{3}= 0, where

*D*

^{+}exhibits homogeneous thermal conductivity

*K*′ and the

*D*

^{−}, homogeneous thermal conductivity

*K*″. At the steady state, the governing equation is written as follows:

*q*

_{V}is the volumetric heat source and the thermal conductivity $K(x)$ can be either

*K*′ or

*K*″ depending on the location $x$. When

*q*

_{V}= 0, the heat equation reduces to Laplace’s equation. When $qV=\delta (x\u2032)$ with $\delta (x\u2032)$ being the Dirac Delta function for a unit point heat source

*Q*, the solution is the Green’s function caused by the unit point source

*Q*at $x\u2032$. Following the imaging displacement setting [33], and modifying Walpole’s formulation [34], the temperature field in the infinite domain can be written in terms of Green’s function as follows [19]:

*Q*(

*x*

_{1}′,

*x*

_{2}′,

*x*

_{3}′), an imaginary point referring to the interface $S$(

*x*′

_{3}= 0) is constructed as the image point of

*Q*as $Q\xaf(x1\u2032,x2\u2032,\u2212x3\u2032)$ to help derive the Green’s function for an infinite bimaterial. For an arbitrary field point of interest

*P*(

*x*

_{1},

*x*

_{2},

*x*

_{3}), the position vectors of

*P*referred to

*Q*and $Q\xaf$ are defined as $r$ and $s$, with their distance given

*P*is located on different sides of the interface due to the superposition of image and infinite parts. Without the loss of generality, we assume the point heat source is above the interface (

*x*′

_{3}> 0). Following Walpole’s formulation [34] and using the interface continuity, we construct Green’s function as follows:

*K*″ =

*K*′, the bimaterial infinite domain recovers a single-material infinite 3D domain, and the Green’s function is the Newtonian potential [23]. When

*K*″ = 0, the lower half exhibits zero heat flux. The bimaterial infinite domain is equivalent to the upper half, or a semi-infinite domain, with zero heat flux on the boundary, which is corresponding to the Mindlin’s solution for the elastic problem [35]. When

*K*″ → ∞, the lower half is perfectly conductive and cannot exhibit any temperature gradient. Therefore, the bimaterial infinite domain is equivalent to the upper half, or a semi-infinite domain, with a constant temperature on the boundary, which is corresponding to the Rongved’s solution [36] for the elastic problem.

When the heat source is periodic, such as daily repetitive temperature change, the Green’s function is a harmonic function with time. When an impulse heat source is considered, the transient Green’s function will be attenuated with time. Both harmonic Green’s function and transient Green’s function for a bimaterial infinite domain can be constructed in the same fashion for both 3D and 2D problems and are provided in Appendix A. For simplicity, this article focuses on the static thermal analysis only, it can be extended to the harmonic and transient heat transfer problems in the same fashion with the Green’s functions in Appendix A.

## 3 One Inhomogeneity Embedded in an Infinite Bimaterial Matrix

*D*with a subdomain Ω centered at $xc$ with an interface $S(x3=0)$ of zero thermal resistance to separate the bimaterial matrix. The particle’s thermal conductivity is

*K*

^{Ω}. According to Mura [28], an inhomogeneity is defined as a subdomain Ω in a matrix with different material properties from the rest (

*D*− Ω), which can be simulated by an eigenstrain. The heat conduction problem can be treated similarly with an ETG on the particle, which is continuous over Ω and can be written in terms of the Taylor expansion of $x\u2212xc$ with the reference point at the particle’s center, such as

*K*′ is rewritten as follows:

*T**

_{i}= 0 for $x\u2208(D\u2212\Omega )$. Substituting the aforementioned heat flux density into the equilibrium equation

*q*

_{i,i}= 0 in the absence of a heat source provides:

*D*terms can be sorted in the general formula,

Here, *D*_{ij}, *D*_{ijk}, and *D*_{ijkl} are written in the terms of $\Phi \xaf$, which represent the components of the image part, and thus, the integral is always outside the particle. Following a straightforward but lengthy procedure [23,28], one can obtain the explicit expressions of Φ_{k,ij}, Φ_{kl,ij} and $\Phi \xafk,ij,\Phi \xafkl,ij$, which are listed in Appendix B. The temperature field is solved in Eq. (13) for the case of semi-infinite domain *D* with thermal conductivity *K*′ containing an inclusion Ω with ETG $Ti*(x)$.

Extending from the inclusion problem to the inhomogeneity problem [23,28] with different thermal conductivity *K*^{Ω}, the heat field can be solved by replacing the inhomogeneity with the same matrix material and applying an ETG field yet to be determined [17]. For one inhomogeneity in an infinite homogeneous domain with uniform far-field heat flux density, the eigen-temperature gradient is uniform by analogy with the case in the elastic field [11,12]. However, in most cases, the uniformity of eigen-fields does not hold due to several factors, such as size effect [31,32] and interactions between inhomogeneities [23]. Even for the case of one inhomogeneity, the existence of boundary effects [24] may disturb eigen-fields as well, which will be illustrated and discussed in Sec. 5 specifically. In such cases, the assumption of uniform eigen-field may result in the poor accuracy of solutions. Hence, in this article, a polynomial-form ETG in Eq. (6) is applied in following case studies of Secs. 4 and 5.

*D*, the equivalent heat flux condition is written as follows:

*i*= 1, 2, 3 with three equations, the second set of equations include free index

*i*and

*r*with 3 × 3 = 9 equations, and the third set of equations with free index

*i*,

*p*, and

*q*include 3 × 6 = 18 equations considering the symmetry of

*p*and

*q*. Overall, the system of linear equations contains 30 unknowns $(Tj0,Tjk1,Tjkl2)$ to be solved. After the ETG field is obtained, the temperature and its gradient can be calculated with Eqs. (10) and (13).

## 4 Numerical Verification of Equivalent Inclusion Method Results With Finite Element Method Simulation

To verify the aforementioned algorithm in Sec. 3, consider an infinite bimaterial domain $D=D+\u222aD\u2212$ composed of two dissimilar isotropic material as shown in Fig. 2 that *K*′ = 1 W/m · K and *K*″ = 10 W/m · K. Without the loss of any generality, the interface *S* is constructed as *x*_{3} = 0 and the domain *D* is subjected to a far-field load $q30=10W/m2$. Consider a spherical inhomogeneity with radius *a* = 1 m and thermal conductivity *K*^{Ω} = 2 W/m · K embedded in the *D*^{+} ($x3c>0$). In Fig. 4, the symmetric properties are utilized for the finite element method (FEM) simulation, and the domain with the length of 20*a* is applied to reduce boundary effects. In the following, heat flux along *x*_{3} direction is compared between FEM and EIM with three cases using uniform, linear, and quadratic ETG, respectively. Note that as an inhomogeneity approaches the interface, boundary effects will disturb the heat fields more significantly. As illustrated in Ref. [24], the boundary effects are generally evaluated as the ratio of distance-to-boundary and radius of the inhomogeneity; therefore, three cases of locations are considered, $x3c=2a,1.5a$, and *a* (inhomogeneity touching interface). Shown in Fig. 5, heat flux of EIM results of three cases are compared with FEM results. When the inhomogeneity is not close to the interface, say 2 *a*, though “uniform” curve exhibits some discrepancies comparing with “linear” and “quadratic” ones, all of them agree well with FEM with difference less than $0.15%$. As the inhomogeneity moves closer toward the interface, deviations are observed in Fig. 5(b) that only “linear” and “quadratic” terms could provide good analysis, especially in neighborhood of $x3\u2212x3c=\u22121$ or 1 m. Considering the domain integral of Green’s function with uniform eigen-field, composed of Φ and $\Phi \xaf$, where Φ is constant within the inhomogeneity domain and $\Phi \xaf$ is the combination of linear and inverse proportional function, which is the reason why using the uniform ETG is barely possible to produce a complicated heat field. Eventually, when the inhomogeneity touches the interface, in which singularity and discontinuity issues arise as the source approaches the boundary. The “quadratic” provides the best predictions among EIM results; however, discrepancies are still observed in the entrance region of the inhomogeneity. Beside the comparison, it is noticed that, at the interface, although the heat flux in the third direction is continuous, its slope changes due to the mismatch of material properties.

## 5 Application to an Underground Heat Storage Tank

To simulate the heat transfer of a spherical underground heat storage tank to the ground, for simplicity, the ground is assumed with a uniform temperature on the surface and in the air and thus the whole system approximated by a bimaterial infinite domain with an infinitely large conductivity above the ground, so that we can focus on the heat conduction within the ground or the other semi-infinite domain.

First, a steady-state heat transfer from the surface into the ground is considered, and the temperature field is disturbed by the existence of the tank due to the conductivity mismatch. A prescribed temperature gradient, or ETG, is introduced to calculate the local temperature field in the neighborhood of the tank. Second, in the heat exchange mode, given a heat exchange demand, the steady-state temperature field is calculated. The temperature difference exists between the water in the tank and soil’s far-field temperature can be derived. The relation between the heat flux, temperature gap, tank size, and depth can be set up as a guideline for GSHP tank design. Parametric studies are conducted for the demonstration of the design features.

Although the earth surface temperature is affected by heat convection with air, heat conduction with ground, and radiation to the sky and from the sun, as well as surface moisture evaporation [5], it generally remains uniform in a planar area with a relatively homogeneous soil. For simplicity, we use a perfectly conductive material with *K*″ → ∞ in the air to simulate the heat transfer in the ground with a uniform surface temperature *T*^{s} and *K*′ = 0.519 W/m · K [3]. Although this assumption simplifies the heat and mass transfer features on the earth surface, it does not change the heat transfer in the neighborhood of the tank very much. In this case, the bimaterial infinite domain is reduced to a semi-infinite domain with a constant temperature on the surface, and the bimaterial Green’s function is reduced to the Green’s function for a semi-infinite domain accordingly.

To promote the heat exchange in the tank, thermal conductive pipe is filled in the tank, which exhibits an effective heat conductivity, significantly higher than the soil, say *K*^{Ω} = 10 − 100 W/m · K. This article only considers the steady-state heat transfer and transient heat conduction problem will be investigated in the future. Two scenarios are considered in the following: (1) disturbance to a uniform heat flux *h*_{3} from the ground surface to the deep earth and (2) heat exchange from the water tank to the surrounding soil.

### 5.1 Heat Flux Disturbance by the Water Tank.

Consider a water tank with radius *a* = 1 m and thermal conductivity *K*^{Ω} embedded underground. Without the loss of any generality, for the steady-state heat conduction in a homogeneous material with a uniform heat flux, a linearly distributed temperature, hence a constant heat flux, say $q30=\u22121w/m2$. In Fig. 6, the water tank is placed at different depths as 2, 4, 6, and 8 m and its thermal conductivity *K*^{Ω} = 10 W/m · K. To compare the disturbance of different thermal conductivity, consider $T\xaf\u2212T30d$, where $T\xaf$ is defined in Eq. (17) as the reference temperature to the ground surface and *d* is the depth of the center of the water tank. Comparing among the four temperature curves in Fig. 6, water tank at the depth of 2 m exhibits larger difference when distance to the center of water tank is approaching 2 m, while the other three curves almost overlap. Such phenomenon can be interpreted as boundary effects of constant temperature and the rapid vanish of second components of kernel function (1/*s*). Shown in Fig. 7, the disturbed effects brought by thermal conductivity is investigated, i.e., *K*^{Ω} = 10, 20, 50 and 100 W/m · K when the depth of the water tank is 4 m. It is noticed that although different thermal conductivity change the heat fields in *x*_{1} and *x*_{2}, the variation is comparatively small compared with that of loading direction (*x*_{3}). The four temperature curves in Fig. 7 mainly differ in [−1, 1] m within the inhomogeneity itself that the smaller *K*^{Ω}, the larger the difference of temperatures, which agrees well with the definition of thermal conductivity. However, the disturbance by *K*^{Ω} is relatively smaller compared with that by depth, where the boundary effects dominate.

### 5.2 Heat Exchange Between the Water Tank and Surrounding Soil.

*a*= 1 m and thermal conductivity

*K*

^{Ω}is constructed underground with distance

*L*to the ground. In summer, besides the power collected by the solar panels, the temperature of the roof also increases due to the solar radiation. To avoid decrease of working efficiency and extend the life span of the panels, the system aims to control the temperature of the roof. In the sense of steady-state heat conduction, say the water tank absorbs heat at the rate of

*Q*= 500 W; hence, the volume heat source

*q*

_{V}≡ 1500/4π W/m

^{3}. Notice that when the point heat source is located underground, Eq. (4) changes as (1/4π

*K*′)(1/

*r*− 1/

*s*) due to

*K*″ → ∞ when

*x*

_{3}< 0. Let the earth exhibit a uniform far-field temperature

*T*

^{0}= 20°C, the temperature field at $x$ is expressed as follows:

*K*

^{Ω}= 10 W/m · K located at different depths, shown in Fig. 8, the temperature fields along

*x*

_{3}and

*x*

_{1}direction are plotted. In such case, the closer the water tank is to the ground surface, the lower the temperature field, which can be interpreted as the confinement of the boundary conditions of constant temperature. When the water tank moves deeper, the boundary effects reduces significantly, so that the temperature differences are much smaller comparing between depths at 2, 4, and 6, 8 m. In Fig. 8(b), except for the case “Depth 2 m,” the other three curves exhibit a comparatively symmetric distribution, which demonstrates that the boundary effect vanishes rapidly with the distance.

In Fig. 9, four cases of different thermal conductivity *K*^{Ω} are investigated when the depth of water tank is 4 m. With higher *K*^{Ω}, the temperature variation within the inhomogeneity domain is smaller; however, the temperature fields in the neighborhood of inhomogeneity are similar.

In Figs. 10(a) and 10(b), the effect brought by radius *a* is presented while keeping the same heat source *Q*. As the water tank shrinks, the volume heat source *q*_{V} increases as proportional to *r*^{−3}; in the extreme case of 0 radius, strong singular effects of heat source exist, which explains significant changes in smaller radius cases.

In Fig. 11, the water tank with radius ranging in [0.2, 2.5] m are placed at a depth of 3, 5, and 8 m and infinity. Compared with boundary effects, the singular effects dominate, which results in similar temperature curves when radii are less than 0.5 m. As discussed earlier, the boundary effects reduces rapidly with distance; hence, the cases with depth equaling 5 and 8 m exhibit small differences. When the depth increases to infinity, the boundary effects by the image part 1/*s* vanish in the order of *r*^{2} [23], which leads to the same analysis in the infinite domain.

In actual geothermal applications, the temperature profile of the ground changes with time and season. The periodic and transient temperature variations should be considered. The present formulation can be extended to harmonic or transient heat conduction problems by the introduction of the Green’s functions for those problems, which have been provided in Appendix A.

## 6 Conclusions

The Green’s function of the steady-state heat conduction problem in a 3D infinite bimaterial is used to solve the inhomogeneity problem that a bimaterial contains a particle with a different thermal conductivity. Eshelby’s equivalent inclusion method is used to obtain the analytical solution, which is verified with the finite element method. Tailorable accuracy can be obtained with the polynomial eigen-temperature gradient. When one semi-infinite domain exhibits a zero or infinitely large thermal conductivity, the solution can be used to study the steady-state heat conduction in the other semi-infinite domain with an isothermal or adiabatic surface. The solution is applied to study the heat transfer of a geothermal tank with the earth with the following conclusions:

When heat is transferred from the surface to the deep earth, the water tank exhibits a higher average temperature than the earth at the same depth, and vice versa.

Boundary effects at the interface dominate the disturbance of heat field when the water tank is close to the ground surface; when the water tank moves gradually deeper, the image part of the Green’s function vanishes rapidly, and the heat field behaves similarly as an infinite homogeneous domain.

Subject to a certain heat source

*Q*, a larger water tank could lower the disturbance to the heat field, while the heat field in a smaller water tank exhibits comparatively significant temperature rises.The thermal conductivity of the inhomogeneity mainly changes the distribution of heat fields within itself; larger thermal conductivity results in “flatter” temperature curves and vice versa.

## Acknowledgment

This work is sponsored by the National Science Foundation IIP #1738802, IIP #1941244, CMMI #1762891, U.S. Department of Agriculture NIFA #2021-67021-34201, and National Natural Science Foundation of China (Grant No. 12102458), whose support is gratefully acknowledged.

## Conflict of Interest

There are no conflicts of interest.

## Data Availability Statement

The datasets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request.

### Appendix A: Other Green’s Functions for a Bimaterial Infinite Domain

Section 2 provides the Green’s function of bimaterials for the 3D infinite domain in the steady state. The method is applicable to the transient heat conduction problem as well, which is more relevant in the actual applications. In the following, we derive the Green’s functions for periodic heat source and transient heat transfer. Moreover, the 3D Green’s functions can be extended to the 2D bimateiral infinite domain. We will provide all Green’s functions. The present model with EIM can be extended to all cases for wide applications.

#### Two-Dimensional Bimaterial Under the Steady-State Heat Conduction

*x*

_{2}axis for both directions from −∞ to +∞. Therefore, the 2D Green’s function can be derived by a direct integral along

*x*

_{2}axis from the 3D alternatives as follows:

*r*is a harmonic function for a two-dimensional case.

#### Green’s Function for Transient Heat Transfer for Bimaterial Domain

*C*

_{p}the specific heat capacity. For simplicity, the heat source is normalized by $\rho Cp$ in comparison with Eq. (1) [39]. When the heat source is a unit impulse source at $(x\u2032,t\u2032)$, let $qV(x,t)=\delta (t\u2212t\u2032)\delta (x\u2212x\u2032)$, the solution of the temperature field can be written in terms of the Green’s function as follows:

*x*

_{1}−

*x*

_{2}plane for the unit impulse source,

*m*= 1, 2, and $\tau =t\u2212t\u2032$. The solution of $T~$ can be written as follows:

*x*

_{1}and

*x*

_{2}directions, the Green’s function of a full-space with an uniform material is $e\u2212r2/4\alpha \tau /4\pi \alpha \tau 3$. Now, considering the interface

*S*is located

*x*

_{3}= 0, one can write the following continuity equations,

*x*

_{3}< 0, the coefficients merge as the expressions are similar. Applying inverse Fourier transform in Eq. (A8) and substituting it into Eq. (A7), in spatial domain, the fundamental solution can be obtained in terms of the Green’s function,

*K*and α denote material properties in

*D*

^{+}and

*D*

^{−}, respectively.

*a*)

*b*)

*D*

^{+}and

*D*

^{−}, respectively.