The Green’s Function-Based Thermal Analysis of a Spherical Geothermal Tank in a Semi-Infinite Domain

Wang, Tengxiang; Wu, Chunlin; Zhang, Liangliang; Yin, Huiming

doi:10.1115/1.4054568

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.

Fig. 1

View large Download slide

Schematic illustration of a spherical thermal energy storage tank embedded in the earth for a bidirectional geothermal system [3]

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

Figure 2 shows an infinite bimaterial D composed of upper D⁺ and lower D⁻ parts, divided by the interface

S

of x₃ = 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:

- \nabla (K (x) \nabla T) = q_{V} (x)

(1)

where 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

q_{V} = δ (x^{'})

with

δ (x^{'})

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^{'}

⁠. 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]:

T (x) = G (x, x^{'}) δ (x^{'})

(2)

where the Green’s function

G (x, x^{'})

needs to be derived with bimaterials in virtue of the image method. In general, the form of

G (x, x^{'})

depends on the location of point of interest

P (x)

and the location of point heat source

Q (x^{'})

⁠. For the point heat source Q(x₁′, x₂′, x₃′), an imaginary point referring to the interface

S

(x′₃ = 0) is constructed as the image point of Q as

\bar{Q} (x_{1}^{'}, x_{2}^{'}, - x_{3}^{'})

to help derive the Green’s function for an infinite bimaterial. For an arbitrary field point of interest P(x₁, x₂, x₃), the position vectors of P referred to Q and

\bar{Q}

are defined as

r

and

s

⁠, with their distance given

\begin{matrix} r = [(x_{1} - x_{1}^{'})^{2} + (x_{2} - x_{2}^{'})^{2} + (x_{3} - x_{3}^{'})^{2}]^{1 / 2} \\ s = [(x_{1} - x_{1}^{'})^{2} + (x_{2} - x_{2}^{'})^{2} + (x_{3} + x_{3}^{'})^{2}]^{1 / 2} \end{matrix}

(3)

Green’s function exhibits different forms when 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′₃ > 0). Following Walpole’s formulation [34] and using the interface continuity, we construct Green’s function as follows:

G (x, x^{'}) = {\begin{matrix} \frac{1}{4 π K^{'} r} + \frac{K^{'} - K^{″}}{K^{'} + K^{″}} \frac{1}{4 π K^{'} s} & for x_{3} \geq 0 \\ \frac{1}{2 π (K^{'} + K^{″}) r} & for x_{3} < 0 \end{matrix}

(4)

Therefore, for a distributed source field

q_{V} (x^{'})

with

x^{'} = (x_{1}^{'}, x_{2}^{'}, x_{3}^{'})

⁠, the temperature at a field point

x = (x_{1}, x_{2}, x_{3})

can be written as follows:

T (x) = \int_{D} G (x, x^{'}) q_{V} (x^{'}) d x^{'}

(5)

When 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.

Fig. 2

View large Download slide

A point heat source Q inside the bimaterial matrix

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

Figure 3 shows an infinite matrix domain D with a subdomain Ω centered at

x^{c}

with an interface

S (x_{3} = 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 - x^{c}

with the reference point at the particle’s center, such as

T_{i}^{*} (x) = T_{i}^{0} + T_{k i}^{1} (x_{k} - x_{k}^{c}) + T_{k l i}^{2} (x_{k} - x_{k}^{c}) (x_{l} - x_{l}^{c}) + \dots, for x \in Ω

(6)

The constitutive law over the upper semi-infinite matrix with thermal conductivity K′ is rewritten as follows:

q_{i} = - K^{'} δ_{i j} (T_{, j} - T_{j}^{*})

(7)

where T*_i = 0 for

x \in (D - Ω)

⁠. Substituting the aforementioned heat flux density into the equilibrium equation q_i,i = 0 in the absence of a heat source provides:

K^{'} T_{, i i} = K^{'} T_{i, i}^{*}

(8)

Fig. 3

View large Download slide

One particle embedded in the bimaterial matrix

By applying Green’s function that satisfies the boundary conditions, the temperature can be written in terms of the eigen-temperature gradient as follows:

T (x) = \int_{D} G (x, x^{'}) (- K^{'}) \frac{\partial T_{i}^{*}}{\partial x_{i}^{'}} d x^{'}

(9)

Applying the Gauss theorem and considering zero eigen-temperature gradient at the boundary of the semi-infinite domain, the aforementioned equation can be rewritten as follows:

T (x) = \int_{Ω} \frac{\partial G (x, x^{'})}{\partial x_{i}^{'}} K^{'} T_{i}^{*} d x^{'} = \int_{Ω} Γ_{i} K^{'} T_{i}^{*} d x^{'}

(10)

where the modified Green’s function

Γ_{i} = - \frac{1}{4 π K^{'}} [{(\frac{1}{r})}_{, i} + \frac{K^{'} - K^{″}}{K^{'} + K^{″}} M_{i j} {(\frac{1}{s})}_{, j}]

(11)

where

M = [\begin{matrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & - 1 \end{matrix}]

provides the mirror projection of a vector. Consider the ETG as a continuous tensor function over the inhomogeneity and thus write it in a polynomial form in Eq. . For simplicity, we take into account of the constant, linear, and quadratic terms only as follows:

T_{i}^{*} (x^{'}) = {\begin{matrix} T_{i}^{0} + T_{i k}^{1} (x_{k}^{'} - x_{k}^{c}) + T_{i k l}^{2} (x_{k}^{'} - x_{k}^{c}) (x_{l}^{'} - x_{l}^{c}), & x^{'} \in Ω \\ 0, & x^{'} \in (D - Ω) \end{matrix}

(12)

Thus, by taking partial differentiation of Eq. and inserting it with Eq. , we can obtain

\begin{matrix} T_{, i} (x) & = \int_{Ω} Γ_{j, i} K^{'} [T_{j}^{0} + T_{j k}^{1} (x_{k}^{'} - x_{k}^{c}) + T_{j k l}^{2} (x_{k}^{'} - x_{k}^{c}) (x_{l}^{'} - x_{l}^{c})] d x^{'} = - \frac{1}{4 π} (D_{i j} T_{j}^{0} + D_{i j k} T_{j k}^{1} + D_{i j k l} T_{j k l}^{2}) \end{matrix}

(13)

where D terms can be sorted in the general formula,

\begin{aligned} D_{i j k l \dots, p q \dots} & = \int_{Ω} [{(\frac{(x_{k}^{'} - x_{k}^{c}) (x_{l}^{'} - x_{l}^{c}) \dots}{r})}_{, i j p q \dots} + \frac{K^{'} - K^{″}}{K^{'} + K^{″}} M_{j m} {(\frac{(x_{k}^{'} - x_{k}^{c}) (x_{l}^{'} - x_{l}^{c}) \dots}{s})}_{, i m p q \dots}] d x^{'} \\ = Φ_{k l \dots, i j p q \dots} + \frac{K^{'} - K^{″}}{K^{'} + K^{″}} M_{j m} M_{k s} M_{l w} \dots {\bar{Φ}}_{s w \dots, i m p q \dots} \end{aligned}

(14)

Here, D_ij, D_ijk, and D_ijkl are written in the terms of $\bar{Φ}$ ⁠, 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 ${\bar{Φ}}_{k, i j}, {\bar{Φ}}_{k l, i j}$ ⁠, 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 $T_{i}^{*} (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.

Considering a uniform far-field heat flux applied on D, the equivalent heat flux condition is written as follows:

K^{Ω} (T_{, i}^{0} + T_{, i}^{'}) = K^{'} (T_{, i}^{0} + T_{, i}^{'} - T_{i}^{*})

(15)

Substituting Eq. into Eq. , the equivalent heat flux conditions with Taylor series expansion can be expressed as follows:

{\begin{matrix} (K^{Ω} - K^{'}) (- 4 π T_{, i}^{0} + [D_{i j} (x^{c}) T_{j}^{0} + D_{i j k} (x^{c}) T_{j k}^{1} + D_{i j k l} (x^{c}) T_{j k l}^{2}]) \\ = 4 π K^{'} T_{i}^{0} \\ (K^{Ω} - K^{'}) [D_{i j, r} (x^{c}) T_{j}^{0} + D_{i j k, r} (x^{c}) T_{j k}^{1} + D_{i j k l, r} (x^{c}) T_{j k l}^{2}] = 4 π K^{'} T_{i r}^{1} \\ (K^{Ω} - K^{'}) [D_{i j, p q} (x^{c}) T_{j}^{0} + D_{i j k, p q} (x^{c}) T_{j k}^{1} + D_{i j k l, p q} (x^{c}) T_{j k l}^{2}] = 8 π K^{'} T_{i p q}^{2} \end{matrix}

(16)

where the first set of equations include free index 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

(T_{j}^{0}, T_{j k}^{1}, T_{j k l}^{2})

to be solved. After the ETG field is obtained, the temperature and its gradient can be calculated with Eqs. and .

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^{+} \cup D^{-}$ 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₃ = 0 and the domain D is subjected to a far-field load $q_{3}^{0} = 10 W / m^{2}$ ⁠. Consider a spherical inhomogeneity with radius a = 1 m and thermal conductivity K^Ω = 2 W/m · K embedded in the D⁺ (⁠ $x_{3}^{c} > 0$ ⁠). In Fig. 4, the symmetric properties are utilized for the finite element method (FEM) simulation, and the domain with the length of 20a is applied to reduce boundary effects. In the following, heat flux along x₃ 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, $x_{3}^{c} = 2 a, 1.5 a$ ⁠, 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 $x_{3} - x_{3}^{c} = - 1$ or 1 m. Considering the domain integral of Green’s function with uniform eigen-field, composed of Φ and $\bar{Φ}$ ⁠, where Φ is constant within the inhomogeneity domain and $\bar{Φ}$ 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.

Fig. 4

View large Download slide

Schematic plot of mesh in FEM analysis of heat flux on one inhomogeneity embedded in infinite bimaterial domain

Fig. 5

Comparison of the heat flux distribution along x3 with distance [ − 3a, 3a] between FEM and EIM for x3c= (a) 1a, (b) 1.5a, and (c) 2a

View large Download slide

Comparison of the heat flux distribution along x₃ with distance [ − 3a, 3a] between FEM and EIM for $x_{3}^{c} =$ (a) 1a, (b) 1.5a, and (c) 2a

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₃ 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 $q_{3}^{0} = - 1 w / m^{2}$ ⁠. 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 $\bar{T} - T_{3}^{0} d$ ⁠, where $\bar{T}$ 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₁ and x₂, the variation is comparatively small compared with that of loading direction (x₃). 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.

Fig. 6

Variation of T¯−T30d subjected to an uniform heat flux q30=−1W/m2 along x3 with distance [−2, 2] m to the center of water tank x3c = −2, −4, −6, and −8 m

View large Download slide

Variation of $\bar{T} - T_{3}^{0} d$ subjected to an uniform heat flux $q_{3}^{0} = - 1 W / m^{2}$ along x₃ with distance [−2, 2] m to the center of water tank $x_{3}^{c}$ = −2, −4, −6, and −8 m

Fig. 7

Variation of T¯−T30d subjected to an uniform heat flux q30=−1W/m2 along x3 with distance [−2, 2] m to the center of water tank xc = −4 m when KΩ = 10, 20, 50, and 100 W/m · K

View large Download slide

Variation of $\bar{T} - T_{3}^{0} d$ subjected to an uniform heat flux $q_{3}^{0} = - 1 W / m^{2}$ along x₃ with distance [−2, 2] m to the center of water tank $x^{c}$ = −4 m when K^Ω = 10, 20, 50, and 100 W/m · K

5.2 Heat Exchange Between the Water Tank and Surrounding Soil.

Shown in Fig. 1, a spherical water tank with radius 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³. Notice that when the point heat source is located underground, Eq. changes as (1/4π K′)(1/r − 1/s) due to K″ → ∞ when x₃ < 0. Let the earth exhibit a uniform far-field temperature T⁰ = 20°C, the temperature field at

x

is expressed as follows:

T (x) = \bar{T} (x) + T^{0} = \int_{D} [G (x, x^{'}) q_{V} (x^{'}) + Γ_{i} (x, x^{'}) T_{i}^{*} (x^{'})] d x^{'} + T^{0}

(17)

where

T_{i}^{*} (x^{'})

is described in Eq. . As demonstrated in Sec. 4, the quadratic term provides the most accurate results compared to uniform and linear ones. Hence, in the following case studies, only quadratic term is applied to illustrate heat fields. Considering the water tank with K^Ω = 10 W/m · K located at different depths, shown in Fig. 8, the temperature fields along x₃ and x₁ 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.

Fig. 8

Variation of temperature subjected to uniform heat source Q = 500 W along (a) x3 and (b) x1 direction with distance [−2, 2] m to the center of water tank x3c = −2, −4, −6, and −8 m

View large Download slide

Variation of temperature subjected to uniform heat source Q = 500 W along (a) x₃ and (b) x₁ direction with distance [−2, 2] m to the center of water tank $x_{3}^{c}$ = −2, −4, −6, and −8 m

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.

Fig. 9

Variation of temperature subjected to uniform heat source Q = 500 W along (a) x3 and (b) x1 direction with distance [−2, 2] m to the center of water tank xc=−4m when KΩ = 10, 20, 50, and 100 W/m · K

View large Download slide

Variation of temperature subjected to uniform heat source Q = 500 W along (a) x₃ and (b) x₁ direction with distance [−2, 2] m to the center of water tank $x^{c} = - 4 m$ when K^Ω = 10, 20, 50, and 100 W/m · K

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⁻³; in the extreme case of 0 radius, strong singular effects of heat source exist, which explains significant changes in smaller radius cases.

Fig. 10

Variation of temperature subjected to uniform heat source Q = 500 W along (a) x3 and (b) x1 directions with distance [−3, 3] m to the center of water tank xc = −8 m with radius a = 0.5, 1, and 2 m

View large Download slide

Variation of temperature subjected to uniform heat source Q = 500 W along (a) x₃ and (b) x₁ directions with distance [−3, 3] m to the center of water tank $x^{c}$ = −8 m with radius a = 0.5, 1, and 2 m

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² [23], which leads to the same analysis in the infinite domain.

Fig. 11

View large Download slide

Variation of temperature subjected to uniform heat source Q = 500 W at the center of water tank $x_{3}^{c} = - 3, - 5$ ⁠, −8 m and infinity with radius ranging in [0.2, 2.5] m

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.

The method can be extended to two-dimensional cases, periodic and transient heat transfer problems for geothermal well operations. The corresponding Green’s function has been provided in Appendix A.

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

To get the two-dimensional (2D) Green’s function for a bimaterial infinite domain, we extend the point heat source Q along x₂ axis for both directions from −∞ to +∞. Therefore, the 2D Green’s function can be derived by a direct integral along x₂ axis from the 3D alternatives as follows:

G^{2 D} (x, x^{'}) = \int_{- \infty}^{\infty} G (x, x^{'}) d x_{2}^{'}

(A1)

where

G (x, x^{'})

is given in Eq. . By using the Hadamard regularization with the finite integral part [37], we can obtain

G^{2 D} (x, x^{'}) = {\begin{matrix} - \frac{1}{2 π K^{'}} \ln r - \frac{K^{'} - K^{″}}{K^{'} + K^{″}} \frac{1}{2 π K^{'}} \ln s & for x_{3} \geq 0 \\ - \frac{1}{π (K^{'} + K^{″})} \ln r & for x_{3} < 0 \end{matrix}

(A2)

where

r = \sqrt{(x_{1} - x_{1}^{'})^{2} + (x_{3} - x_{3}^{'})^{2}}

and

s = \sqrt{(x_{1} - x_{1}^{'})^{2} + (x_{3} + x_{3}^{'})^{2}}

⁠. The Green’s function derived from three-dimensional case can be verified with two-dimensional equilibrium equation of heat conduction as we know lnr is a harmonic function for a two-dimensional case.

Green’s Function for Transient Heat Transfer for Bimaterial Domain

To obtain the fundamental solution of a transient heat transfer problem, the Fourier and Laplace transforms are generally applied to change the governing equation of a partial differential equation into an algebraic equation in the transformed domain. Using the corresponding inverse transforms, the fundamental solution in spatial and temporal domain can be determined. For bimaterials, Zhou and Han [38] applied the Fourier transform of horizontal plane and Laplace transform in temporal domain to derive the Green’s function for anisotropic bimaterial domain, and the Green’s function in the spatial and time domain exists in forms of numerical inverse transforms. The closed-form Green’s function exhibits significant value for the robustness and simplicity. Here, we follow their procedure and derive the explicit form of the Green’s function for practical applications with the EIM. Consider the heat equation in a full-space isotropic domain,

\dot{T} (x, t) - α (x) T_{, i i} (x, t) = q_{V} (x, t)

(A3)

where

α (x) \equiv K / ρ C_{p}

is the thermal diffusivity with

ρ

the density and C_p the specific heat capacity. For simplicity, the heat source is normalized by

ρ C_{p}

in comparison with Eq. [39]. When the heat source is a unit impulse source at

(x^{'}, t^{'})

⁠, let

q_{V} (x, t) = δ (t - t^{'}) δ (x - x^{'})

⁠, the solution of the temperature field can be written in terms of the Green’s function as follows:

T ({x - x}^{'}, t - t^{'}) = G ({x, x}^{'}, t, t^{'})

(A4)

which can be used for general distributed heat source by the spatial and time integrals. Considering the initial temperature distribution is 0, to derive the explicit form of the Green’s function, the Fourier transform is applied to Eq. in the x₁−x₂ plane for the unit impulse source,

\dot{\tilde{T}} - α (- \tilde{T} ξ_{m} ξ_{m} + {\tilde{T}}_{, 33}) = e^{- i ξ_{m} x_{m}^{'}} δ (x_{3} - x_{3}^{'}) δ (τ)

(A5)

where

\tilde{T} = \tilde{T} (ξ_{1}, ξ_{2}, x_{3}, t)

⁠, dummy index m = 1, 2, and

τ = t - t^{'}

⁠. The solution of

\tilde{T}

can be written as follows:

\tilde{T} (ξ_{1}, ξ_{2}, x_{3}, t) = C e^{- ξ_{m} ξ_{m} t α} \frac{e^{- x_{3}^{2} / 4 α t}}{\sqrt{4 π α t}} + e^{i ξ_{m} x_{m}^{'} - ξ_{m} ξ_{m} τ α} \frac{e^{- (x_{3} - x_{3}^{'})^{2} / 4 α τ}}{\sqrt{4 π α τ}}

(A6)

When

τ \to 0

the last term reduces to

δ (x_{3} - x_{3}^{'})

⁠. Applying the inverse Fourier transform in x₁ and x₂ directions, the Green’s function of a full-space with an uniform material is

e^{- r^{2} / 4 α τ} / {\sqrt{4 π α τ}}^{3}

⁠. Now, considering the interface S is located x₃ = 0, one can write the following continuity equations,

T (x, t) |_{x_{3} = 0^{+}} = T (x, t) |_{x_{3} = 0^{-}} and q_{3} (x, t) |_{x_{3} = 0^{+}} = q_{3} (x, t) |_{x_{3} = 0^{-}}

(A7)

From Zhou and Han [38] and Walpole [34], the solution is the superposition of a full-space solution and an image part. In such case, the solution can be written as follows:

\tilde{T} (ξ_{1}, ξ_{2}, x_{3}, t) = {\begin{matrix} {\tilde{T}}^{\infty} (ξ_{1}, ξ_{2}, x_{3}, t) + m_{1} \bar{\tilde{T}} (ξ_{1}, ξ_{2}, x_{3}, t) & for x_{3} \geq 0 \\ m_{2} \tilde{T} (ξ_{1}, ξ_{2}, x_{3}, t) & for x_{3} < 0 \end{matrix}

(A8)

where

\bar{\tilde{T}}

represents the function with image source; when x₃ < 0, the coefficients merge as the expressions are similar. Applying inverse Fourier transform in Eq. and substituting it into Eq. , in spatial domain, the fundamental solution can be obtained in terms of the Green’s function,

\begin{aligned} G ({x, x}^{'}, t, t^{'}) & = \\ {\begin{matrix} \frac{1}{(4 π α^{'} τ)^{3 / 2}} {e^{- r^{2} / 4 α^{'} τ} + \frac{K^{'} α^{″} - K^{″} α^{'}}{K^{'} α^{″} + K^{″} α^{'}} e^{- s^{2} / 4 α^{'} τ}} & for x_{3} \geq 0 \\ \frac{1}{(4 π α^{″} τ)^{3 / 2}} \frac{2 K^{'} (α^{″})^{5 / 2} (α^{'})^{- 3 / 2}}{K^{'} α^{″} + K^{″} α^{'}} e^{- r_{s}^{2} / 4 α^{'} τ} e^{(- r^{2} + r_{s}^{2}) / 4 α^{″} τ} & for x_{3} < 0 \end{matrix} \end{aligned}

(A9)

where

r_{s} = \sqrt{(x_{1} - x_{1}^{'})^{2} + (x_{2} - x_{2}^{'})^{2} + (x_{3}^{'})^{2}}

⁠; superscripts ′ and ″ of K and α denote material properties in D⁺ and D⁻, respectively.

When the heat load is harmonic, i.e.,

q_{V} = δ ({x - x}^{'}) e^{i ω t}

⁠, the fundamental solution can be derived using same continuity conditions in Eq. and the solution can be similarly expressed in terms of the Green’s function as follows:

G (x, x^{'}, ω) = {\begin{matrix} {\tilde{T}}^{\infty} (x_{1}, x_{2}, x_{3}, ω) + m_{3} \bar{\tilde{T}} (x_{1}, x_{2}, x_{3}, ω) & for x_{3} \geq 0 \\ m_{4} \tilde{T} (x_{1}, x_{2}, x_{3}, ω), & for x_{3} < 0 \end{matrix}

(A10)

where

m_{3} = \frac{(1 + i) (K^{'} - K^{″}) + 2 (c^{'} K^{'} - c^{″} K^{″}) r_{s}}{(1 + i) (K^{'} + K^{″}) + 2 (c^{'} K^{'} + c^{″} K^{″}) r_{s}}

(A11a)

m_{4} = \frac{2 K^{'} ((1 + i) + 2 c^{'} α^{″} r_{s}) e^{(i - 1) r_{s} (c^{'} - c^{″})}}{α^{'} ((1 + i) (K^{'} + K^{″}) + 2 (c^{'} K^{'} + c^{″} K^{″}) r_{s})}

(A11b)

where

c = \sqrt{ω / 2 α}

⁠.

{\tilde{T}}^{\infty}

and

\tilde{T}

are the Green’s function for heat equation with harmonic heat load with material properties D⁺ and D⁻, respectively.

{\tilde{T}}^{\infty} (x, x^{'}, ω) = \frac{1}{4 π α^{'} r} e^{i ω t} e^{(i - 1) c^{'} r}

(A12)

and

\bar{\tilde{T}}

is for image source,

\bar{\tilde{T}} (x, x^{'}; ω) = \frac{1}{4 π α^{'} s} e^{i ω t} e^{(i - 1) c^{'} s}

(A13)

Notice that for harmonic heat load, we have assumed

e^{i ω t}

as time-related function. Yet the heat field is real, hence the imaginary part should be dropped.

Appendix B: The Integral of the Harmonic Potential and Its Derivatives

The second-order partial derivatives of domain integrals are as follows,

\begin{aligned} Φ_{, i j} & = {\begin{matrix} \frac{4 π a^{3}}{3 r^{3}} (δ_{i j} - 3 n_{i} n_{j}) & for x \in (D - Ω) \\ - \frac{4 π}{3} δ_{i j} & for x \in Ω \end{matrix} \end{aligned}

(B1)

\begin{aligned} Φ_{k, i j} & = {\begin{matrix} \frac{4 π a^{5}}{5 r^{4}} (- δ_{i j} n_{k} - δ_{i k} n_{j} - δ_{j k} n_{i} + 5 n_{i} n_{j} n_{k}), & for x \in (D - Ω) \\ \frac{4 π}{5} [- δ_{i j} (x_{k} - x_{k}^{c}) - δ_{i k} (x_{j} - x_{j}^{c}) - δ_{j k} (x_{i} - x_{i}^{c})], & for x \in Ω \end{matrix} \end{aligned}

(B2)

\begin{aligned} Φ_{k l, i j} & = {\begin{matrix} \frac{16 π a^{7}}{5 r^{5}} n_{i} n_{j} n_{k} n_{l} - \frac{16 π a^{7}}{35 r^{5}} (δ_{i j} n_{k} n_{l} + δ_{k j} n_{i} n_{l} + δ_{l j} n_{i} n_{k}) \\ - \frac{20 π a^{7}}{35 r^{5}} (δ_{k i} n_{j} n_{l} + δ_{l i} n_{j} n_{k}) + \frac{4 π a^{7}}{35 r^{5}} (δ_{l j} δ_{k i} + δ_{l i} δ_{k j}) \\ + \frac{π a^{5}}{105 r^{5}} [(12 a^{2} - 28 r^{2}) δ_{i j} δ_{k l} + (84 r^{2} - 60 a^{2}) δ_{k l} n_{i} n_{j}] & for x \in (D - Ω) \\ - \frac{4 π r^{2}}{7} (δ_{i k} n_{j} n_{l} + δ_{l i} n_{j} n_{k} + δ_{i j} n_{l} n_{k} + δ_{k j} n_{i} n_{l} + δ_{l j} n_{l} n_{k}) \\ + (\frac{2 π a^{2}}{5} - \frac{2 π r^{2}}{7}) (δ_{i k} δ_{j l} + δ_{i l} δ_{j k}) \\ + (\frac{4 π r^{2}}{35} - \frac{4 π a^{2}}{15}) δ_{i j} δ_{k l} + \frac{8 π r^{2}}{35} n_{i} n_{j} δ_{k l} & for x \in Ω \end{matrix} \end{aligned}

(B3)

Notice that for the image potentials, i.e.,

\bar{Φ}

can be obtained by using the three aforementioned equations but switching the source point as

\bar{x}

⁠.

The third-order partial derivatives of domain integrals are as follows:

\begin{aligned} Φ_{, i j r} & = {\begin{matrix} \frac{4 π a^{3}}{r^{4}} (5 n_{i} n_{j} n_{r} - δ_{i j} n_{r} - δ_{i r} n_{j} - δ_{j r} n_{i}) & for r > a \\ 0 & for r \leq a \end{matrix} \end{aligned}

(B4)

\begin{aligned} (B 5) & Φ_{, i j p q} & = {\begin{matrix} \frac{4 π a^{3}}{r^{5}} [- 35 n_{i} n_{j} n_{p} n_{q} + 5 (δ_{i j} n_{p} n_{q} + δ_{i p} n_{j} n_{q} + δ_{i q} n_{j} n_{p} + δ_{i p} n_{i} n_{q} + δ_{j q} n_{i} n_{p} + δ_{p q} n_{i} n_{j}) - δ_{i j} δ_{p q} - δ_{i p} δ_{j q} - δ_{j p} δ_{i q}] for r > a \\ 0 for r \leq a \end{matrix} \end{aligned}

(B5)

\begin{aligned} Φ_{k, i j r} & = {\begin{matrix} \frac{4 π a^{5}}{5 r^{5}} [- 35 n_{i} n_{j} n_{r} n_{k} + 5 (δ_{i j} n_{k} n_{r} + δ_{i k} n_{j} n_{r} + δ_{i r} n_{j} n_{k} + δ_{j k} n_{i} n_{r} + δ_{j r} n_{i} n_{k} + δ_{k r} n_{i} n_{j}) - δ_{i j} δ_{k r} - δ_{i k} δ_{j r} - δ_{j k} δ_{i r}] for r > a \\ \frac{4 π}{5} (- δ_{i j} δ_{k r} - δ_{i k} δ_{j r} - δ_{j k} δ_{i r}) for r \leq a \end{matrix} \end{aligned}

(B6)

\begin{aligned} Φ_{k, i j p q} & = {\begin{matrix} \frac{4 π a^{5}}{r^{6}} [63 n_{i} n_{j} n_{p} n_{q} n_{k} - 7 (δ_{i j} n_{p} n_{q} n_{k} + δ_{i k} n_{j} n_{p} n_{q} + δ_{i p} n_{j} n_{q} n_{k} + δ_{i q} n_{j} n_{p} n_{k} \\ + δ_{j k} n_{i} n_{p} n_{q} + δ_{j p} n_{i} n_{q} n_{k} + δ_{j q} n_{i} n_{p} n_{k} + δ_{k p} n_{i} n_{j} n_{q} + δ_{k q} n_{i} n_{j} n_{p} + δ_{p q} n_{i} n_{j} n_{k}) \\ + (δ_{j k} δ_{p q} + δ_{j p} δ_{k q} + δ_{k p} δ_{j q}) n_{i} + (δ_{i k} δ_{p q} + δ_{i p} δ_{k q} + δ_{k p} δ_{i q}) n_{j} + (δ_{i j} δ_{p q} + δ_{i p} δ_{j q} + δ_{j p} δ_{i q}) n_{k} \\ + (δ_{i j} δ_{k q} + δ_{i k} δ_{j q} + δ_{j k} δ_{i q}) n_{p} + (δ_{i j} δ_{k p} + δ_{i k} δ_{j p} + δ_{j k} δ_{i p}) n_{q}] for r > a \\ 0 for r \leq a \end{matrix} \end{aligned}

(B7)

Φ_{k l, i j r} = {\begin{matrix} \frac{- 144 π a^{7}}{5 r^{6}} n_{i} n_{j} n_{k} n_{l} n_{r} + \frac{16 π a^{7}}{5 r^{6}} (δ_{i j} n_{k} n_{l} n_{r} + δ_{l j} n_{i} n_{k} n_{r} + δ_{k j} n_{i} n_{l} n_{r} + δ_{i r} n_{j} n_{l} n_{k} + δ_{j r} n_{i} n_{l} n_{k} + δ_{k r} n_{i} n_{j} n_{l}) \\ + \frac{4 π a^{7}}{r^{6}} (δ_{k i} n_{j} n_{l} n_{r} + δ_{l i} n_{j} n_{k} n_{r}) + (\frac{4 π a^{7} - 4 π a^{5} r^{2}}{r^{6}}) δ_{k l} n_{i} n_{j} n_{r} \\ + [\frac{- 16 π a^{7}}{35 r^{6}} (δ_{i j} δ_{l r} + δ_{l j} δ_{i r}) - \frac{20 π a^{7}}{35 r^{6}} δ_{l i} δ_{j r}] n_{k} + [\frac{- 16 π a^{7}}{35 r^{6}} (δ_{i j} δ_{k r} + δ_{k j} δ_{i r}) - \frac{20 π a^{7}}{35 r^{6}} δ_{k i} δ_{j r}] n_{l} \\ + [\frac{- 16 π a^{7}}{35 r^{6}} (δ_{j k} δ_{l r} + δ_{k r} δ_{j l}) + \frac{π a^{5} (28 r^{2} - 20 a^{2})}{35 r^{6}} δ_{k l} δ_{j r}] n_{i} + [\frac{- 20 π a^{7}}{35 r^{6}} (δ_{k i} δ_{l r} + δ_{l i} δ_{k r}) + \frac{π a^{5} (28 r^{2} - 20 a^{2})}{35 r^{6}} δ_{k l} δ_{i r}] n_{j} \\ + [\frac{- 4 π a^{7}}{35 r^{6}} (δ_{l j} δ_{k i} + δ_{l i} δ_{k j}) + \frac{π a^{5} (28 r^{2} - 20 a^{2})}{35 r^{6}} δ_{k l} δ_{i j}] n_{r} for r > a \\ - \frac{4 π r}{7} (δ_{i k} δ_{j r} + δ_{i j} δ_{k r} + δ_{k j} δ_{i r}) n_{l} - \frac{4 π r}{7} (δ_{l i} δ_{j r} + δ_{i j} δ_{l r} + δ_{l j} δ_{i r}) n_{k} \\ + [\frac{- 4 π r}{7} (δ_{j k} δ_{l r} + δ_{k r} δ_{j l}) + \frac{8 π r}{35} δ_{k l} δ_{j r}] n_{i} + [\frac{- 4 π r}{7} (δ_{i k} δ_{l r} + δ_{k r} δ_{i l}) + \frac{8 π r}{35} δ_{k l} δ_{i r}] n_{j} \\ + [\frac{- 4 π r}{7} (δ_{i k} δ_{l j} + δ_{j k} δ_{i l}) + \frac{8 π r}{35} δ_{k l} δ_{i j}] n_{r} for r \leq a \end{matrix}

(B8)

(B9)

For 2D cases, the integrals of the harmonic potential and its derivatives can be written as follows:

\begin{matrix} Φ = {\begin{matrix} 2 π a^{2} \ln r & for x \in (D - Ω) \\ π r^{2} + π a^{2} [2 \ln (a) - 1] & for x \in Ω \end{matrix} \\ Φ_{, i} = {\begin{matrix} - 2 π a^{2} x_{i} / r^{2} = \frac{- 2 π a^{2}}{r} n_{i} & for x \in (D - Ω) \\ - 2 π x_{i} = - 2 π r n_{i} & for x \in Ω \end{matrix} \\ Φ_{p, i} = {\begin{matrix} \frac{π}{4} δ_{p i} (2 a^{2} - r^{2}) - \frac{π}{2} r^{2} n_{p} n_{i} & for x \in (D - Ω) \\ \frac{π a^{4}}{4 r^{2}} δ_{p i} - \frac{π a^{4}}{2 r^{2}} n_{p} n_{i} & for x \in Ω \end{matrix} \\ Φ_{m n, i} = {\begin{matrix} \frac{π}{2} (δ_{m i} x_{n} + δ_{n i} x_{m}) (a^{2} - \frac{2 r^{2}}{3}) - \frac{2 π x_{m} x_{n} x_{i}}{3} - δ_{m n} x_{i} (\frac{π a^{2}}{2} - \frac{π r^{2}}{6}) & for x \in (D - Ω) \\ (δ_{m i} x_{n} + δ_{n i} x_{m}) \frac{π a^{6}}{6 r^{4}} - x_{i} x_{m} x_{n} \frac{2 π a^{6}}{3 r^{6}} - δ_{m n} x_{i} (\frac{π a^{4}}{2 r^{2}} - \frac{π a^{6}}{6 r^{4}}) & for x \in Ω \end{matrix} \end{matrix}

(B10)

References

1.

Yin

,

H.

,

Yang

,

D.

,

Kelly

,

G.

, and

Garant

,

J.

,

2013

, “

Design and Performance of a Novel Building Integrated Pv/thermal System for Energy Efficiency of Buildings

,”

Sol. Energy

,

87

, pp.

184

–

195

.

Google Scholar

Crossref

2.

Yin

,

H.

,

Zadshir

,

M.

, and

Pao

,

F.

,

2021

,

Building Integrated Photovoltaic Thermal Systems: Fundamentals, Designs and Applications

,

Academic Press

,

Cambridge, MA

.

3.

Yumrutaş

,

R.

, and

Ünsal

,

M.

,

2012

, “

Energy Analysis and Modeling of a Solar Assisted House Heating System With a Heat Pump and an Underground Energy Storage Tank

,”

Sol. Energy

,

86

(

3

), pp.

983

–

993

.

Google Scholar

Crossref

4.

Ingersoll

,

L. R.

,

1948

, “

Theory of the Ground Pipe Heat Source for the Heat Pump

,”

Heat. Pip. Air Condition.

,

20

, pp.

119

–

122

.

5.

Badache

,

M.

,

Eslami-Nejad

,

P.

,

Ouzzane

,

M.

,

Aidoun

,

Z.

, and

Lamarche

,

L.

,

2016

, “

A New Modeling Approach for Improved Ground Temperature Profile Determination

,”

Renewable Energy

,

85

(

10

), pp.

436

–

444

.

Google Scholar

Crossref

6.

Pan

,

E.

, and

Chen

,

W.

,

2015

,

Static Green’s Functions in Anisotropic Media

,

Cambridge University Press

,

Cambridge, UK

.

7.

Klein

,

B.

,

Register

,

L. F.

,

Hess

,

K.

,

Deppe

,

D.

, and

Deng

,

Q.

,

1998

, “

Self-Consistent Green’s Function Approach to the Analysis of Dielectrically Apertured Vertical-Cavity Surface-Emitting Lasers

,”

Appl. Phys. Lett.

,

73

(

23

), pp.

3324

–

3326

.

Google Scholar

Crossref

8.

Hanson

,

G. W.

,

2008

, “

Dyadic Green’s Functions and Guided Surface Waves for a Surface Conductivity Model of Graphene

,”

J. Appl. Phys.

,

103

(

6

), p.

064302

.

Google Scholar

Crossref

9.

Greenberg

,

M. D.

,

2015

,

Applications of Green’s Functions in Science and Engineering

,

Dover Publications

,

Mineola, NY

.

10.

Cole

,

K.

,

Beck

,

J.

,

Haji-Sheikh

,

A.

, and

Litkouhi

,

B.

,

2010

,

Heat Conduction Using Greens Functions

,

CRC Press

,

Boca Raton, FL

.

Google Scholar

Crossref

11.

Eshelby

,

J. D.

, and

Peierls

,

R. E.

,

1957

, “

The Determination of the Elastic Field of an Ellipsoidal Inclusion, and Related Problems

,”

Proc. R. Soc. London, A

,

241

(

1226

), pp.

376

–

396

.

Google Scholar

Crossref

12.

Eshelby

,

J. D.

, and

Peierls

,

R. E.

,

1959

, “

The Elastic Field Outside an Ellipsoidal Inclusion

,”

Proc. R. Soc. London, A

,

252

(

1271

), pp.

561

–

569

.

Google Scholar

Crossref

13.

Zhang

,

L.

,

Lin

,

Q.

,

Chen

,

F.

,

Zhang

,

Y.

, and

Yin

,

H.

,

2020

, “

Micromechanical Modeling and Experimental Characterization for the Elastoplastic Behavior of a Functionally Graded Material

,”

Int. J. Solids Struct.

,

206

, pp.

370

–

382

.

Google Scholar

Crossref

14.

Sakata

,

S.

,

Ashida

,

F.

, and

Kojima

,

T.

,

2010

, “

Stochastic Homogenization Analysis for Thermal Expansion Coefficients of Fiber Reinforced Composites Using the Equivalent Inclusion Method With Perturbation-Based Approach

,”

Comput. Struct.

,

88

(

7–8

), pp.

458

–

466

.

Google Scholar

Crossref

15.

Takei

,

T.

,

Hatta

,

H.

, and

Taya

,

M.

,

1991

, “

Thermal Expansion Behavior of Particulate-Filled Composites II: Multi-Reinforcing Phases (Hybrid Composites)

,”

Mater. Sci. Eng. A

,

131

(

1

), pp.

145

–

152

.

Google Scholar

Crossref

16.

Yin

,

H.

,

Paulino

,

G.

,

Buttlar

,

W.

, and

Sun

,

L.

,

2007

, “

Micromechanics-Based Thermoelastic Model for Functionally Graded Particulate Materials With Particle Interactions

,”

J. Mech. Phys. Solids

,

55

(

1

), pp.

132

–

160

.

Google Scholar

Crossref

17.

Hatta

,

H.

, and

Taya

,

M.

,

1986

, “

Equivalent Inclusion Method for Steady State Heat Conduction in Composites

,”

Int. J. Eng. Sci.

,

24

, pp.

1159

–

1170

.

Google Scholar

Crossref

18.

Yin

,

H. M.

,

Paulino

,

G. H.

,

Buttlar

,

W. G.

, and

Sun

,

L. Z.

,

2008

, “

Effective Thermal Conductivity of Functionally Graded Particulate Nanocomposites With Interfacial Thermal Resistance

,”

ASME J. Appl. Mech.

,

75

(

5

), p.

051113

.

Google Scholar

Crossref

19.

Carslaw

,

H. S.

, and

Jaeger

,

J. C.

,

1986

,

Conduction of Heat in Solids

, 2nd ed.,

Oxford University Press

,

Oxford, New York

.

20.

Koroteeva

,

O.

,

Mogilevskaya

,

S.

,

Crouch

,

S.

, and

Gordeliy

,

E.

,

2010

, “

A Computational Technique for Evaluating the Effective Thermal Conductivity of Isotropic Porous Materials

,”

Eng. Anal. Bound. Elem.

,

34

(

9

), pp.

793

–

801

.

Google Scholar

Crossref

21.

Mogilevskaya

,

S.

,

Kushch

,

V.

,

Koroteeva

,

O.

, and

Crouch

,

S.

,

2012

, “

Equivalent Inhomogeneity Method for Evaluating the Effective Conductivities of Isotropic Particulate Composites

,”

J. Mech. Mater. Struct.

,

7

(

1

), pp.

103

–

117

.

Google Scholar

Crossref

22.

Wu

,

C.

,

Wei

,

Z.

, and

Yin

,

H.

,

2021

, “

Virtual and Physical Experiments of Encapsulated Phase Change Material Embedded in Building Envelopes

,”

Int. J. Heat Mass Transfer

,

172

(

6

), p.

121083

.

Google Scholar

Crossref

23.

Yin

,

H.

,

Song

,

G.

,

Zhang

,

L.

, and

Wu

,

C.

,

2022

,

The Inclusion-Based Boundary Element Method

,

Academic Press

,

Cambridge, MA

.

24.

Liu

,

Y.

,

Song

,

G.

, and

Yin

,

H.

,

2015

, “

Boundary Effect on the Elastic Field of a Semi-infinite Solid Containing Inhomogeneities

,”

Proc. R. Soc. A

,

471

(

2179

), p.

20150174

.

Google Scholar

Crossref

25.

Dang

,

X.

,

Liu

,

Y.

,

Wang

,

L.

, and

Wang

,

J.

,

2019

, “

Solutions of the Elastic Fields in a Half-Plane Region Containing Multiple Inhomogeneities With the Equivalent Inclusion Method and the Applications to Properties of Composites

,”

Acta Mech.

,

230

(

5

), pp.

1529

–

1547

.

Google Scholar

Crossref

26.

Ostoja-Starzewski

,

M.

, and

Schulte

,

J.

,

1996

, “

Bounding of Effective Thermal Conductivities of Multiscale Materials by Essential and Natural Boundary Conditions

,”

Phys. Rev. B

,

54

(

1

), pp.

278

–

285

.

Google Scholar

Crossref

27.

Hill

,

R.

,

1963

, “

Elastic Properties of Reinforced Solids: Some Theoretical Principles

,”

J. Mech. Phys. Solids

,

11

(

5

), pp.

357

–

372

.

Google Scholar

Crossref

28.

Mura

,

T.

,

1987

,

Micromechanics of Defects in Solids

, 2nd ed. (

Mechanics of Elastic and Inelastic Solids

),

Springer, Netherlands

.

Google Scholar

Crossref

29.

Nemat-Nasser

,

S.

, and

Hori

,

M.

,

2013

,

Micromechanics: Overall Properties of Heterogeneous Materials

,

Elsevier

,

Amsterdam, The Netherlands

.

30.

Willis

,

J. R.

,

1981

,

Advances in Applied Mechanics

, Vol.

21

,

C.-S.

Yih

, ed.,

Elsevier

,

Amsterdam, The Netherlands

, pp.

1

–

78

.

31.

Wu

,

C.

, and

Yin

,

H.

,

2021

, “

Elastic Solution of a Polygon-Shaped Inclusion With a Polynomial Eigenstrain

,”

ASME J. Appl. Mech.

,

88

(

6

), p.

061002

.

Google Scholar

Crossref

32.

Wu

,

C.

,

Zhang

,

L.

, and

Yin

,

H.

,

2021

, “

Elastic Solution of a Polyhedral Particle With a Polynomial Eigenstrain and Particle Discretization

,”

ASME J. Appl. Mech.

,

88

(

12

), p.

121001

.

Google Scholar

Crossref

33.

Liu

,

Y. J.

,

Song

,

G.

, and

Yin

,

H. M.

,

2015

, “

Boundary Effect on the Elastic Field of a Semi-infinite Solid Containing Inhomogeneities

,”

Proc. R. Soc. A

,

471

(

2179

), p.

20150174

.

Google Scholar

Crossref

34.

Walpole

,

L.

,

1997

, “

An Inclusion in One of Two Joined Isotropic Elastic Half-Spaces

,”

IMA J. Appl. Math.

,

59

(

2

), pp.

193

–

209

.

Google Scholar

Crossref

35.

Mindlin

,

R. D.

,

1936

, “

Force at a Point in the Interior of a Semi-Infinite Solid

,”

Physics

,

7

(

5

), pp.

195

–

202

.

Google Scholar

Crossref

36.

Rongved

,

L.

,

1955

, “

Force Interior to One of Two Joined Semi-Infinite Solid

,”

Proceedings of the Second Midwestern Conference on Solid Mechanics

, pp.

1

–

13

.

37.

Hadamard

,

J.

,

1932

,

Le Probléme De Cauchy Et Les équations Aux Dérivées Partielles Linéaires Hyperboliques

,

Hermann & Cie

,

Paris

.

38.

Zhou

,

J.

, and

Han

,

X.

,

2020

, “

Three-Dimensional Green’s Functions for Transient Heat Conduction Problems in Anisotropic Bimaterial

,”

Int. J. Heat Mass Transfer

,

146

(

4

), p.

118805

.

Google Scholar

Crossref

39.

Haberman

,

R.

,

2012

,

Applied Partial Differential Equations With Fourier Series and Boundary Value Problems

,

Addison-Wesley

,

Boston, MA

.

The Green’s Function-Based Thermal Analysis of a Spherical Geothermal Tank in a Semi-Infinite Domain

Abstract

1 Introduction

2 Formulation of Green’s Function

3 One Inhomogeneity Embedded in an Infinite Bimaterial Matrix

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

5 Application to an Underground Heat Storage Tank

5.1 Heat Flux Disturbance by the Water Tank.

5.2 Heat Exchange Between the Water Tank and Surrounding Soil.

6 Conclusions

Acknowledgment

Conflict of Interest

Data Availability Statement

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

Two-Dimensional Bimaterial Under the Steady-State Heat Conduction

Green’s Function for Transient Heat Transfer for Bimaterial Domain

Appendix B: The Integral of the Harmonic Potential and Its Derivatives

References

Get Email Alerts

Cited By

ASME Journals

ASME Conference Proceedings

ASME eBooks

Resources

Opportunities

The Green’s Function-Based Thermal Analysis of a Spherical Geothermal Tank in a Semi-Infinite Domain

Abstract

1 Introduction

2 Formulation of Green’s Function

3 One Inhomogeneity Embedded in an Infinite Bimaterial Matrix

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

5 Application to an Underground Heat Storage Tank

5.1 Heat Flux Disturbance by the Water Tank.

5.2 Heat Exchange Between the Water Tank and Surrounding Soil.

6 Conclusions

Acknowledgment

Conflict of Interest

Data Availability Statement

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

Two-Dimensional Bimaterial Under the Steady-State Heat Conduction

Green’s Function for Transient Heat Transfer for Bimaterial Domain

Appendix B: The Integral of the Harmonic Potential and Its Derivatives

References

Get Email Alerts

Cited By

Related Articles

Related Proceedings Papers

Related Chapters

ASME Journals

ASME Conference Proceedings

ASME eBooks

Resources

Opportunities

This Feature Is Available To Subscribers Only