Explicit Multipole Formulas for Calculating Thermal Resistance of Single U-Tube Ground Heat Exchangers ^†

Abstract

Borehole thermal resistance is both an important design parameter and a key performance characteristic of a ground heat exchanger. Another quantity that is particularly important for ground heat exchangers is the internal thermal resistance between the heat exchanger pipes. Both these resistances can be calculated to a high degree of accuracy by means of the well-known multipole method. However, the multipole method has a fairly intricate mathematical algorithm and is thus not trivial to implement. Consequently, there is considerable interest in developing explicit formulas for calculating borehole resistances. This paper presents derivation and solutions of newly derived second-order and higher-order multipole formulas for calculating borehole thermal resistance and total internal thermal resistance of single U-tube ground heat exchangers. A new and simple form of the first-order multipole formula is also presented. The accuracy of the presented formulas is established by comparing them to the original multipole method. The superiority of the new higher-order multipole formulas over the existing formulas is also demonstrated.

1. Introduction

The use of ground source heat pump (GSHP) systems to provide heating and cooling in buildings has increased at a rapid rate in the last two decades or so. Stimulated by energy prices, technology advances and environmental concerns, the thermal energy utilized by these systems has increased from approximately 500 GWh in 1995 to over 91,000 GWh in 2015. During the same period, the worldwide installed capacities of these systems have increased from under 2000 MWt in 1995 to over 50,000 MWt in 2015 [1]. Comparable growth levels are expected for the future.

A typical GSHP system consists of a heat pump, a ground heat exchanger, and auxiliary systems for storage and distribution of thermal energy. The sizing and performance of a heat pump system depends upon the coefficient of performance of the heat pump, which in turn is a function of heat source and heat sink temperatures. Compared to ambient air, ground, in general, is a far superior source or sink of thermal energy because of its relatively stable temperature levels over the year. In a GSHP system, a ground heat exchanger is used to transfer heat to or from the ground. The ground heat exchanger can be of open or closed type. In an open system groundwater is directly used as the heat carrier fluid, whereas in a closed system the heat carrier fluid is circulated in a closed loop, which can be horizontal or vertical. Various heat exchanger configurations can be used in closed-loop vertical systems, including single or double U-tubes, and simple or complex coaxial pipes. Among all types, a borehole heat exchanger with a single U-tube is by far the most commonly used ground heat exchanger in practice because of its low cost, small space requirements, and ease of installation. The scope of this paper is also limited to the application of single U-tubes in borehole heat exchangers.

Ground thermal conductivity (λ) and borehole thermal resistance (R_b) are the two principal parameters that govern the heat transfer mechanism of a borehole heat exchanger, and thus influence the sizing and performance of the overall GSHP system [2]. The heat transfer outside the borehole boundary is dictated by the thermal conductivity of the ground, whereas the heat transfer inside the borehole is characterized by the borehole thermal resistance between the heat carrier fluid and the borehole wall. A high ground thermal conductivity is beneficial for the ground heat transfer. However, being an intrinsic property of the ground, ground thermal conductivity cannot be controlled in practice. On the other hand, a low borehole thermal resistance is desirable for better heat transfer inside the borehole heat exchanger. The borehole thermal resistance depends upon the physical arrangement and the thermal properties of borehole components including grouting, heat exchanger pipes, and the heat carrier fluid. Its value can be engineered to a certain extent by optimizing the geometry and layout of the borehole heat exchanger and by choosing appropriate materials for the borehole components.

For any given borehole configuration, the borehole thermal resistance can either be estimated theoretically [3] or measured experimentally [4]. The theoretical methods to estimate borehole thermal resistance include analytical [5,6,7] or empirical [8,9,10,11] formulas based on one-dimensional or two-dimensional steady-state conductive heat transfer in the borehole. As demonstrated by [12,13], there exists a great disparity between borehole thermal resistance values calculated from different theoretical methods. Many of the simplified methods work well in certain situations and poorly in others.

The multipole method is an analytical method to determine thermal resistances for any number of arbitrarily placed pipes in a composite region. It is based on two-dimensional steady-state conductive heat transfer in a borehole, and uses a combination of line heat sources and so-called multipoles to determine thermal resistances for any number of arbitrarily placed pipes in a composite region. The method was originally developed by [14,15], and was later revised by [16]. The main difference between the versions of [14] and [16] was in regards to prescribing a constant temperature condition at a suitable radial distance outside the borehole. In the later version, this condition was replaced by a condition on the mean temperature around the borehole wall, which simplified the algorithm.

The multipole method has been compared to and tested against several numerical methods. Young [17] compared the tenth-order multipole method with a boundary-fitted coordinates finite volume methods for 18 combinations of shank spacing and grout conductivities, and observed a maximum difference of 0.26% between the two methods. Lamarche et al. [12] calculated a root-mean-square difference of 0.003 between the first-order multipole method and a finite element numerical model for 72 studied cases with varying grout and soil conductivities, and pipe sizes and positing. He [18] studied four cases of different borehole diameters and grout conductivities using a three-dimensional numerical model and the multipole method, and reported a maximum difference of 0.19%. Liao et al. [11] compared the first-order multipole method with a two-dimensional numerical method and noted a maximum relative difference of 0.2% for 18 combinations of shank spacing and grout conductivities. Al-Chalabi [19] studied 12 different cases of pipe positioning, and grout and soil conductivities, using the first-order multipole method and a two-dimensional finite element method, and observed a maximum difference of 1.3% between the two methods.

The accuracy and the complexity of the multipole method increases with the number of multipoles used for the calculation. When implemented in a computer program, the order of the multipoles to be used for a calculation is typically prescribed to ten, which was the maximum possible order in the original implementation [14] of the multipole method. Popular ground heat exchanger programs Earth Energy Designer (EED, v3.2. BLOCON: Lund, Sweden) [20] and Ground Loop Heat Exchanger Professional (GLHEPro) [21] also use tenth-order multipoles when calculating the borehole thermal resistance. The tenth-order multipole calculations have an accuracy of over eight decimal digits [22].

On the adverse side, the actual multipole method has a quite rigorous mathematical formulation and a fairly complex algorithm. Its implementation in computer programs requires a considerable amount of coding (the original implementation [14] in FORTRAN was nearly 600 lines in length). However, it is possible to simplify the multipole method to explicit closed-form formulas assuming that heat exchanger pipes are placed symmetrically about the center of the borehole. In reality, heat exchanger pipes may be located anywhere in the borehole as long as they do not overlap each other. The position of pipes also varies along the depth of the borehole. However, in the absence of any a priori knowledge of the pipe positions, it is reasonable to assume that the heat exchanger pipes are symmetrically placed about the center of the borehole. So far, closed-form multipole formulas for zeroth-order and first-order have been developed for the case of a single U-tube with symmetrical pipes [23].

In this paper, we present derivation and solutions of new second-order and higher-order multipole formulas for calculating borehole thermal resistance of single U-tube ground heat exchangers with symmetric pipes. The presented formulas are tested against the original multipole method (i.e., the tenth-order multipole calculation) as well as previously-derived zeroth-order and first-order multipole formulas.

2. Borehole Thermal Resistance

The borehole thermal resistance (R_b) is the ratio of the temperature difference between the heat carrier fluid and the borehole wall to the heat transfer rate per unit length of the borehole. It is defined locally at a specific depth in the borehole as represented in Equation (1), where T_f,loc is the local mean fluid temperature, T_b is the average borehole wall temperature and q_b is the heat transfer rate per unit length of the borehole.

$$R_{\mathrm{b}}=\frac{T_{\mathrm{f},\mathrm{loc}}-T_{\mathrm{b}}}{q_{\mathrm{b}}}$$

(1)

Borehole thermal resistance is sometimes treated as a sum of resistances in series. The total borehole resistance (R_b) i.e., fluid-to-ground resistance is considered to be made up of two major parts: pipe resistance (R_p) and grout resistance (R_g). The pipe resistance includes both conductive resistance of the pipe (R_pc) and convective resistance of the fluid (R_pic). The grout resistance constitutes the thermal resistance between the outer pipe wall of the U-tube and the borehole wall. Its value depends on the pipe resistance and the ground thermal conductivity [13]. As the individual pipe resistances are in parallel with each other, the pipe resistance is generally calculated for a single pipe and is then divided by the total number (N) of pipes.

$$R_{\mathrm{b}}=R_{\mathrm{g}}+\frac{R_{\mathrm{p}}}{N}=R_{\mathrm{g}}+\frac{R_{\mathrm{pc}}+R_{\mathrm{pic}}}{N}$$

(2)

Nevertheless, in order to fully understand the concept of borehole thermal resistance, it is helpful to use a resistance network. Figure 1 shows an example of such a network for a single U-tube borehole. Other networks forms can also be formulated, e.g., see [11,24] but any such representation is an approximation to reality under network-specific assumptions and restrictions.

The three resistance Δ network of Figure 1 is based on the following relations between heat flows q₁ and q₂ and fluid temperatures T_f1 and T_f2.

$$q_1=\frac{T_{\mathrm{f}1}-T_{\mathrm{b},\mathrm{av}}}{R_{1\mathrm{b}}}+\frac{T_{\mathrm{f}1}-T_{\mathrm{f}2}}{R_{12}}$$

(3)

$$q_2=\frac{T_{\mathrm{f}2}-T_{\mathrm{b},\mathrm{av}}}{R_{2\mathrm{b}}}+\frac{T_{\mathrm{f}2}-T_{\mathrm{f}1}}{R_{12}}$$

(4)

In the above equations R_1b and R_2b are thermal resistances between pipe 1 or 2 and the borehole wall, respectively, and R₁₂ is the fluid-to-fluid thermal resistance between pipes 1 and 2. These resistances include the thermal resistances of the fluid and the pipe. Using the Δ network of Figure 1, the total borehole resistance, R_b, can be determined by setting the fluid temperatures T_f1 = T_f2 and solving for two parallel resistances R_1b and R_2b:

$$R_{\mathrm{b}}=\frac{R_{1\mathrm{b}}{R}_{2\mathrm{b}}}{R_{1\mathrm{b}}+R_{2\mathrm{b}}}$$

(5)

The Δ network can also be used to determine resistance R_a by setting the heat flows q₁ = −q₂. The resistance, R_a, is the total internal resistance between the upward and downward flowing legs of the U-tube when there is no net heat flow from the borehole. It is related to the fluid-to-fluid internal resistance R₁₂ by the thermal network of Figure 1. The resistance R₁₂, which is sometimes also referred to as the direct coupling resistance, is the effect of the network representation and is not a directly measurable physical resistance. Resistance R₁₂ is in parallel with the series resistances R_1b and R_2b.

$$R_{\mathrm{a}}=\frac{R_{12}(R_{1\mathrm{b}}+R_{2\mathrm{b}})}{R_{12}+R_{1\mathrm{b}}+R_{2\mathrm{b}}}$$

(6)

The concept of total internal resistance R_a is critical to the understanding of the effective borehole thermal resistance $R_{\mathrm{b}}^{*}$, which is the thermal resistance between the heat carrier fluid, characterized by a simple mean T_f of the inlet and outlet temperatures, and the average borehole wall temperature T_b. The effective borehole thermal resistance is mathematically expressed as Equation (7).

$$R_{\mathrm{b}}^{*}=\frac{T_{\mathrm{f}}-T_{\mathrm{b}}}{q_{\mathrm{b}}}$$

(7)

The difference between borehole thermal resistance R_b, defined by Equation (1), and effective borehole thermal resistance $R_{\mathrm{b}}^{*}$, defined by Equation (7), is that the borehole thermal resistance applies locally at a specific depth in the borehole, whereas effective borehole thermal resistance applies to the entire borehole. The effective borehole thermal resistance is what is measured in an experiment, for example an in-situ thermal response test. The experimental measurements of mean fluid temperature taken at the top of the borehole include the effects of thermal short-circuiting between the upward and downward flows in the U-tube. Consequently, the effective borehole thermal resistance is often higher than the borehole thermal resistance due to the higher fluid temperature caused by the thermal short-circuiting between the U-tube pipes.

It is possible to calculate the effective borehole thermal resistance from the borehole thermal resistance [23,25,26]. Equations (8) and (9), developed by [23], are commonly used for this purpose. These equations correspond to two distinct boundary conditions of uniform borehole wall temperature and uniform heat flux along the borehole, respectively. Both these are limiting conditions and the real situation falls somewhere in between. Therefore, the effective borehole thermal resistance is commonly determined as the mean value between these two equations. Calculation of the effective borehole thermal resistance from Equations (8) and (9) requires knowledge of total internal thermal resistance R_a and direct coupling resistance R₁₂ between the two U-tube legs, respectively, in addition to the borehole thermal resistance R_b.

$$R_{\mathrm{b}}^{*}=R_{\mathrm{b}}+\frac{1}{3R_{\mathrm{a}}}{\left(\frac{H}{{\rho }_{\mathrm{f}}{c}_{\mathrm{f}}{V}_{\mathrm{f}}}\right)}^2$$

(8)

$$R_{\mathrm{b}}^{*}=R_{\mathrm{b}}\eta \cdot \coth \eta ,\eta =\frac{H}{{\rho }_{\mathrm{f}}{c}_{\mathrm{f}}{V}_{\mathrm{f}}}\cdot \frac{1}{2R_{\mathrm{b}}}\sqrt{1+\frac{4R_{\mathrm{b}}}{R_{12}}}$$

(9)

3. Multipole Method for N Pipes in a Borehole

The general multipole method may be applied for any number, N, of pipes in the borehole as shown in Figure 2. The center of pipe n lies at (x_n,y_n). The thermal conductivity of the grout region outside the pipes is λ_b (W/m-K), and of the ground outside the boreholes is λ. The borehole radius is r_b and the outer radius of the pipes is r_p. The thermal resistance from the fluid in the pipe to the grout adjacent to the pipe is R_p (m-K/W).

The (steady-state) temperature field T(x,y) is to be determined for any set of prescribed heat fluxes q_n, n = 1, …, N, and any prescribed average temperature T_b,av around the borehole wall. The fluid temperature in pipe n is T_fn, n = 1, …, N. From this the borehole thermal resistances are determined. The temperature field consists of a linear combination of line heat source solutions and so-called multipole solutions at the center of each pipe:

$$T(x,y)=T_{\mathrm{b},\mathrm{av}}+{\displaystyle \sum _{n=1}^N\frac{q_n}{2\pi {\lambda }_{\mathrm{b}}}\cdot \mathrm{Re}\left[W_0(z,z_n)\right]}+\mathrm{Re}\left[{\displaystyle \sum _{n=1}^N{\displaystyle \sum _{j=1}^J{P}_{n,j}\cdot W_j(z,z_n)}}\right]$$

(10)

The following complex notations are used:

$$z=x+i\cdot y,\overline{z}=x-i\cdot y,z_n=x_n+i\cdot y_n,n=1,\dots ,N.$$

(11)

The real part of the function W₀ (z,z_n) gives the temperature field for a line heat source at z = z_n with the strength $q_n=2\pi {\lambda }_{\mathrm{b}}$(W/m):

$$W_0(z,z_n)=\{\begin{array}{cc}\ln \left(\frac{r_{\mathrm{b}}}{z-z_n}\right)+\sigma \cdot \ln \left(\frac{r_{\mathrm{b}}^2}{r_{\mathrm{b}}^2-\overline{z}\cdot z_n}\right)& \left|z\right|\le r_{\mathrm{b}}\\ \left(1+\sigma \right)\cdot \ln \left(\frac{r_{\mathrm{b}}}{z-z_n}\right)+\sigma \cdot \frac{1+\sigma }{1-\sigma }\cdot \ln \left(\frac{r_{\mathrm{b}}}{z}\right)& \left|z\right|\ge r_{\mathrm{b}}\end{array}$$

(12)

The upper expression of Equation (12) is valid in the borehole outside the pipes and the lower one in the ground outside the borehole. The first top logarithm represents an ordinary line heat source, while the second term represents a line source at the mirror point $z=r_{\mathrm{b}}^2/\overline{z_{\mathrm{b}}}$. The two expressions, and the ensuing radial heat flux at the grout and ground side, are constructed to become equal at the borehole radius $\left|z\right|=r_{\mathrm{b}}$. This means that the solution satisfies the internal boundary conditions at the borehole wall. The average temperature T_b,av around the borehole wall is zero.

The last term in Equation (10) for the temperature field involves a multipole solution of order j at z = z_n:

$$W_j(z,z_n)=\{\begin{array}{cc}{\left(\frac{r_{\mathrm{p}}}{z-z_n}\right)}^j+\sigma \cdot {\left(\frac{r_{\mathrm{p}}\cdot \overline{z}}{r_{\mathrm{b}}^2-\overline{z}\cdot z_n}\right)}^j& \left|z\right|\le r_{\mathrm{b}}\\ \left(1+\sigma \right){\left(\frac{r_{\mathrm{p}}}{z-z_n}\right)}^j& \left|z\right|\ge r_{\mathrm{b}}\end{array}j=1,2,\dots$$

(13)

The above solutions satisfy the internal boundary conditions of continuous temperature and radial heat flux at the borehole wall. The average temperature at the borehole wall is zero. Multipoles up to order j = J are used. The first top term in the above expressions, in local polar coordinates around the outer periphery of pipe n, becomes:

$${\left(\frac{r_{\mathrm{p}}}{z-z_n}\right)}^j=\left[z-z_n=\rho \cdot e^{i\cdot \psi }\right]=\frac{r_{\mathrm{p}}^j}{{\rho }^j}\cdot e^{-i\cdot \psi \cdot j}=\begin{array}{c}\\ \left[\rho =r_{\mathrm{p}}\right]=\cos (j\cdot \psi )-i\cdot \sin (j\cdot \psi ),-\pi <\psi \le \pi .\\ \end{array}$$

(14)

The multipole strength factors P_n,j are complex-valued numbers. This means the expression (10) can account for a Fourier expansion of any temperatures around the pipes with Cosine and Sine terms up to order j = J.

The boundary condition at the periphery of pipe n involves a balance between the radial heat flux and the temperature difference between fluid temperature T_fn and the temperature outside the pipe, which varies around the pipe. Let as in Equation (14), $\rho$ be the radial distance from the center of pipe n and ψ the polar angle. The temperature in the grout at the pipe wall and the radial heat flux multiplied by the pipe resistance per unit pipe area are given by the two left-hand terms below:

$${T(x,y)|}_{\rho =r_{\mathrm{p}},\psi }+\underset{-\beta \cdot r_{\mathrm{p}}}{\underbrace{R_{\mathrm{p}}\cdot 2\pi r_{\mathrm{p}}\cdot (-{\lambda }_{\mathrm{b}})}}\cdot {\frac{\partial }{\partial \rho }\left[{T(x,y)|}_{\rho ,\psi }\right]|}_{\rho =r_{\mathrm{p}}}={\tilde{T}}_{\mathrm{f}n}(\psi )\simeq T_{\mathrm{f}n}.$$

(15)

The left-hand expression will depend on the polar angle ψ. The exact boundary condition is that the expression becomes constant and equal to the required fluid temperature T_fn to obtain the prescribed heat fluxes q_n.

The temperature field of Equation (10), taken for any set of multipole factors P_n,j, is inserted in Equation (15). It is possible by quite lengthy calculation involving complex-valued formulations and Taylor expansions to separate the expressions in Fourier exponentials of any order j. The Fourier coefficient of order j for any pipe n depends on the heat fluxes and multipole factors:

$${\tilde{T}}_{\mathrm{f}n}(\psi )=\mathrm{Re}\left[F_{n,0}(q,P)+{\displaystyle \sum _{j=1}^{\infty }{F}_{n,j}(q,P)\cdot e^{i\cdot \psi \cdot j}}\right].$$

(16)

Here, the vector q represents the N prescribed heat fluxes q_n, and the matrix P represents the N·J complex-valued multipole factors P_n,j. The multipole factors are now chosen so that the Fourier coefficients up to order J become zero. This gives N·J equations for the multipole factors:

$$F_{n,j}(q,P)=0,n=1,\dots ,N,j=1,\dots ,J\to P_{n,j}$$

(17)

The equations are complex-valued with a real and an imaginary part. This means that the number of unknowns and equations are 2·N·J for the general case. The zero-order term in Equation (16) determines the required fluid temperature T_fn, and the terms above j = J represent the error for the heat balance over the pipe wall as a function of the polar angle ψ:

$$T_{\mathrm{f}n}=\mathrm{Re}\left[F_{0,n}(q,P)\right],{\tilde{T}}_{\mathrm{f}n}(\psi )-T_{\mathrm{f}n}=\mathrm{Re}\left[{\displaystyle \sum _{j=J+1}^{\infty }{F}_{k,n}(q,P)\cdot e^{i\cdot \psi \cdot j}}\right].$$

(18)

4. Two Symmetrically Placed Pipes in a Borehole

Figure 3 shows the case that is considered in this study. There are two symmetrically placed pipes in a grouted borehole. The center of the two pipes lie at (x_p,0) and (−x_p,0), respectively. The prescribed heat fluxes from the pipes are q₁ and q₂ (W/m). The corresponding fluid temperatures are T_f1 and T_f2. The main objective is to present a very precise method to calculate the relations between heat flows and fluid temperatures.

The primary input data are:

$${\lambda }_{\mathrm{b}},\lambda ,r_{\mathrm{b}},x_{\mathrm{p}},r_{\mathrm{p}},R_{\mathrm{p}},q_1,q_2,J,T_{\mathrm{b},\mathrm{av}}.$$

(19)

The thermal resistance may be represented by a dimensionless parameter β and the ratio between the two conductivities by the parameter σ. The circles of the pipes and the borehole should not intersect. This gives:

$$r_{\mathrm{p}}\le x_{\mathrm{p}}\le r_{\mathrm{b}}-r_{\mathrm{p}};\beta =2\pi {\lambda }_{\mathrm{b}}{R}_{\mathrm{p}},\sigma =\frac{{\lambda }_{\mathrm{b}}-\lambda }{{\lambda }_{\mathrm{b}}+\lambda }.$$

(20)

The number of equations and unknowns for multipole factors for two pipes are in general 2·2·J. But in the considered symmetrical case, the multipole factors P_n,j turn out to be real-valued numbers. The number of unknowns and equations are reduced to 2·J for the two sets of multipoles P_1,j and P_2,j.

4.1. Multipole Relations for Even and Odd Solutions

The multipoles for the two pipes depend on the prescribed heat fluxes. The temperature field is symmetric in x for equal heat fluxes:

$$Evencase:s=+1,q_2=s\cdot q_1,T_{\mathrm{b},\mathrm{av}}=0\Rightarrow T(x,y)=T(-x,y).$$

(21)

The temperature field is odd in x for antisymmetric heat fluxes:

$$Oddcase:s=-1,q_2=s\cdot q_1,T_{\mathrm{b},\mathrm{av}}=0\Rightarrow T(x,y)=-T(-x,y).$$

(22)

The average temperature at the borehole wall is chosen to be zero in both cases: T_b,av = 0.

There is a set of multipole factors for the even case and another one for the odd case. It is found for these two cases of high symmetry that the multipoles of pipe 2 are closely related to the multipoles of pipe 1 in the even and odd cases. The following relations are valid for the even case s = +1 and for the odd case s = −1:

$$P_{2,j}^{\pm }=s\cdot {(-1)}^j\cdot P_{1,j}^{\pm },\{\begin{array}{l}+\leftrightarrow s=1\\ -\leftrightarrow s=-1,\end{array}j=1,\dots ,J.$$

(23)

The upper index + represents the even temperature solution and s = +1, while the upper index − represents the odd solution and s = −1. With these notational conventions, the formulas below for the even and odd cases can be presented in the same formulas, where s = ±1 becomes a parameter.

The temperature solutions for the even and odd cases determine the required fluid temperatures, which define an even and an odd thermal resistance $R_J^{+}$ and $R_J^{-}$, respectively, in the following way:

$$Evencase{q}_2=q_1,T_{\mathrm{b},\mathrm{av}}=0:T_{\mathrm{f}2}^{+}=T_{\mathrm{f}1}^{+},T_{\mathrm{f}1}^{+}=R_J^{+}\cdot q_1.$$

(24)

$$Oddcase{q}_2=-q_1,T_{\mathrm{b},\mathrm{av}}=0:T_{\mathrm{f}2}^{-}=-T_{\mathrm{f}1}^{-},T_{\mathrm{f}1}^{-}=R_J^{-}\cdot q_1.$$

(25)

These even and odd thermal resistances are obtained from the multipole solutions. They are the required fluid temperatures for unit heat flux. These thermal resistances depend on the choice of the number of multipoles J. It will be shown below that the accuracy increases rapidly with increasing J.

4.2. Formulas for Even and Odd Thermal Resistances

The temperature field of Equation (10) may be used without the right-hand sum of multipoles. In this case J = 0, and only the sum of line sources are used. The thermal resistances for the two symmetrically placed pipes then become:

$$R_0^{+}=R_{\mathrm{p}}+\frac{1}{2\pi {\lambda }_{\mathrm{b}}}\cdot \left[\ln \left(\frac{r_{\mathrm{b}}^2}{2r_{\mathrm{p}}{x}_{\mathrm{p}}}\right)+\sigma \cdot \ln \left(\frac{r_{\mathrm{b}}^4}{r_{\mathrm{b}}^4-x_{\mathrm{b}}^4}\right)\right],$$

(26)

$$R_0^{-}=R_{\mathrm{p}}+\frac{1}{2\pi {\lambda }_{\mathrm{b}}}\cdot \left[\ln \left(\frac{2x_{\mathrm{p}}}{r_{\mathrm{p}}}\right)+\sigma \cdot \ln \left(\frac{r_{\mathrm{b}}^2+x_{\mathrm{p}}^2}{r_{\mathrm{b}}^2-x_{\mathrm{p}}^2}\right)\right].$$

(27)

The correction to the zero-order resistance for any positive J is given by a formula of the type:

$$R_J^{\pm }=R_0^{\pm }-B_J^{\pm },s=\pm 1,J=1,2,3,\dots .$$

(28)

The quantities $B_J^{+}$ and $B_J^{-}$ are obtained from the even and odd solutions for the considered multipole order J. In the formulas below the following auxiliary (dimensionless) parameters are used:

$$p_0=\frac{r_{\mathrm{p}}}{2x_{\mathrm{p}}},p_1=\frac{r_{\mathrm{p}}{x}_{\mathrm{p}}}{r_{\mathrm{b}}^2-x_{\mathrm{p}}^2},p_2=\frac{r_{\mathrm{p}}{x}_{\mathrm{p}}}{r_{\mathrm{b}}^2+x_{\mathrm{p}}^2},b_k=\frac{1-k\beta }{1+k\beta },k=1,2,\dots .$$

(29)

The B-values for any positive J are obtained from the formula:

$$B_J^{\pm }=\frac{1}{2\pi {\lambda }_b}\cdot {\left(V_J^{\pm }\right)}^{\mathrm{tr}}·\left[{\left(M_J^{\pm }\right)}^{-1}\right]·V{b}_J^{\pm },s=\pm 1,J=1,2,3,\dots .$$

(30)

The larger dot indicates a matrix/vector multiplication. The formula involves two vectors V_J and Vb_J, and a matrix M_J for the even case s = +1 and the odd case s = −1:

$${\left(V_J^{\pm }\right)}^{tr}=\left(\begin{array}{ccc}{V}_1^{\pm }& \cdots & V_J^{\pm }\end{array}\right),M_J^{\pm }=\left(\begin{array}{ccc}{M}_{1,1}^{\pm }& \cdots & M_{1,J}^{\pm }\\ ⋮& \ddots & ⋮\\ M_{J,1}^{\pm }& \cdots & M_{J,J}^{\pm }\end{array}\right),V{b}_J^{\pm }=\left(\begin{array}{c}{b}_1\cdot V_1^{\pm }\\ ⋮\\ b_J\cdot V_J^{\pm }\end{array}\right).$$

(31)

Note that inverses of the M matrices appear in Equation (30). The components of the vectors V_J and Vb_J are given by:

$$V_k^{\pm }={(-p_0)}^k\cdot s+\sigma \cdot p_1^k+\sigma \cdot {(-p_2)}^k\cdot s,V{b}_k^{\pm }=b_k\cdot V_k^{\pm },k=1,\dots ,J.$$

(32)

The M matrices are given by:

$$M_{k,j}^{\pm }=k\cdot {\delta }_{k,j}+b_k\cdot A_{k,j}^{\pm },k=1,\dots ,J,j=1,\dots ,J;{\delta }_{k,j}=\{\begin{array}{cc}1& j=k\\ 0& j\ne k\end{array}.$$

(33)

Here, the two matrices A⁺ and A⁻, which account for various interactions between line sources and multipoles in the boundary conditions at the two pipes, are fairly complicated:

$$\begin{array}{l}{A}_{k,j}^{\pm }=\frac{(k+j-1)!}{(k-1)!\cdot (j-1)!}\cdot {(-p_0)}^{k+j}\cdot s+\sigma kj\cdot {\displaystyle \sum _{i=0}^{\min (k,j)}\frac{(k+j-1-i)!}{i!\cdot (k-i)!\cdot (j-i)!}\cdot {\left(2p_0\right)}^i}\left[p_1^{k+j-i}+{(-p_2)}^{j+k-i}\cdot s\right],\\ \\ A_{k,j}^{\pm }=A_{j,k}^{\pm },s=\pm 1,k=1,\dots ,J,j=1,\dots ,J.\end{array}$$

(34)

It may be noted that the A matrices are symmetric.

The matrix formula, Equation (30), becomes a sum over all components of the inverse of the M matrix:

$$B_J^{\pm }=\frac{1}{2\pi {\lambda }_{\mathrm{b}}}\cdot {\displaystyle \sum _{k=1}^J{\displaystyle \sum _{j=1}^J{V}_k^{\pm }\cdot {\left[{\left(M_J^{\pm }\right)}^{-1}\right]}_{k,j}}}\cdot b_j\cdot V_j^{\pm },J=1,2,3,\dots ,s=\pm 1.$$

(35)

4.3. Explicit Multipole Formulas for First, Second and Third-Order

The zero-order resistances $R_0^{+}$ and $R_0^{-}$ are given by the explicit formulas of Equations (26) and (27), respectively. Explicit formulas may also be derived for J = 1 and J = 2.

The first-order $B_1^{\pm }$ corrections are given by Equation (35) for J = 1, in which case the inverse of the matrix is equal the inverse of the single matrix component:

$$B_1^{\pm }=\frac{1}{2\pi {\lambda }_{\mathrm{b}}}\cdot \frac{V_1^{\pm }\cdot b_1\cdot V_1^{\pm }}{M_{1,1}^{\pm }}.$$

(36)

The quantities in this formula are defined in Equations (32)–(34):

$$\begin{array}{l}{V}_1^{\pm }=-p_0\cdot s+\sigma p_1-\sigma p_2\cdot s,V{b}_1^{\pm }=b_1\cdot V_1^{\pm },s=\pm 1,\\ \\ M_{1,1}^{\pm }=1+b_1\cdot A_{1,1}^{\pm },A_{1,1}^{\pm }=p_0^2\cdot s+\sigma \cdot \left[p_1\cdot \left(p_1+2p_0\right)+p_2\cdot \left(p_2-2p_0\right)\cdot s\right].\end{array}$$

(37)

The final formula is then:

$$B_1^{\pm }=\frac{1}{2\pi {\lambda }_{\mathrm{b}}}\cdot \frac{b_1\cdot {\left(-p_0\cdot s+\sigma p_1-\sigma p_2\cdot s\right)}^2}{1+b_1\cdot \left\{p_0^2\cdot s+\sigma \cdot \left[p_1\cdot \left(p_1+2p_0\right)+p_2\cdot \left(p_2-2p_0\right)\cdot s\right]\right\}}.$$

(38)

Here, the five dimensionless parameters are defined by Equations (20) and (29).

The second-order $B_2^{\pm }$ corrections are given by Equation (35) for J = 2:

$$B_2^{\pm }=\frac{1}{2\pi {\lambda }_{\mathrm{b}}}\cdot {\left(V_2^{\pm }\right)}^{tr}·\left[{\left(M_2^{\pm }\right)}^{-1}\right]·V{b}_2^{\pm }.$$

(39)

The M matrices and their inverses are:

$$M_2^{\pm }=\left(\begin{array}{cc}1+b_1{A}_{1,1}^{\pm }& b_1{A}_{1,2}^{\pm }\\ b_2{A}_{1,2}^{\pm }& 2+b_2{A}_{2,2}^{\pm }\end{array}\right),{\left(M_2^{\pm }\right)}^{-1}=\frac{1}{\det (M_2^{\pm })}\cdot \left(\begin{array}{cc}2+b_2{A}_{2,2}^{\pm }& -b_1{A}_{1,2}^{\pm }\\ -b_2{A}_{1,2}^{\pm }& 1+b_1{A}_{1,1}^{\pm }\end{array}\right).$$

(40)

The first-order components in Equation (37) for J = 1 are still valid. The second-order components of the symmetric A matrices and V vectors are from Equations (32) and (34):

$$\begin{array}{l}{A}_{1,2}^{\pm }=A_{2,1}^{\pm }=-2p_0^3\cdot s+2\sigma \cdot \left[p_1^2\cdot \left(p_1+2p_0\right)-p_2^2\cdot \left(p_2-2p_0\right)\cdot s\right],\\ \\ A_{2,2}^{\pm }=6p_0^4\cdot s+2\sigma \cdot \left[p_1^2\left(3p_1^2+8p_1^{}{p}_0+4p_0^2\right)+p_2^2\left(3p_2^2-8p_2^{}{p}_0+4p_0^2\right)\cdot s\right],\\ \\ V_2^{\pm }=p_0^2\cdot s+\sigma \cdot p_1^2+\sigma \cdot p_2^2\cdot s.\end{array}$$

(41)

A more explicit formula can now be obtained from Equations (35), (39) and (40):

$$B_2^{\pm }=\frac{1}{2\pi {\lambda }_{\mathrm{b}}}\cdot \frac{b_1{\left(V_1^{\pm }\right)}^2\cdot \left(2+b_2{A}_{2,2}^{\pm }\right)-2b_1{b}_2{V}_1^{\pm }{V}_2^{\pm }\cdot A_{1,2}^{\pm }+b_2{\left(V_2^{\pm }\right)}^2\cdot \left(1+b_1{A}_{1,1}^{\pm }\right)}{\left(1+b_1{A}_{1,1}^{\pm }\right)\cdot \left(2+b_2{A}_{2,2}^{\pm }\right)-b_1{b}_2{\left(A_{1,2}^{\pm }\right)}^2}.$$

(42)

The third-order $B_J^{\pm }$ corrections for J = 3 involve the inversion of the 3 × 3 M matrices. A more explicit formula cannot be given. The third-order $B_3^{\pm }$ corrections are given by Equation (30) for J = 3:

$$B_3^{\pm }=\frac{1}{2\pi {\lambda }_{\mathrm{b}}}\cdot {\left(V_3^{\pm }\right)}^{tr}·\left[{\left(M_3^{\pm }\right)}^{-1}\right]·V{b}_3^{\pm }.$$

(43)

Here, $V_k^{\pm }$ and $V{b}_k^{\pm }$ are given by Equation (31) for J = 3. The M matrices have the form:

$$M_3^{\pm }=\left(\begin{array}{ccc}1+b_1{A}_{1,1}^{\pm }& b_1{A}_{1,2}^{\pm }& b_1{A}_{1,3}^{\pm }\\ b_2{A}_{1,2}^{\pm }& 2+b_2{A}_{2,2}^{\pm }& b_2{A}_{2,3}^{\pm }\\ b_3{A}_{1,3}^{\pm }& b_3{A}_{2,3}^{\pm }& 3+b_3{A}_{3,3}^{\pm }\end{array}\right).$$

(44)

Here, $A_{k,j}^{\pm }$ is from Equation (34) for indices k or j equal to three. The third-order $B_3^{\pm }$ corrections of Equation (43) are given by Equation (45) for J = 3:

$$B_3^{\pm }=\frac{1}{2\pi {\lambda }_{\mathrm{b}}}\cdot {\displaystyle \sum _{k=1}^3{\displaystyle \sum _{j=1}^3{V}_k^{\pm }\cdot {\left[{\left(M_3^{\pm }\right)}^{-1}\right]}_{k,j}}}\cdot b_j\cdot V_j^{\pm }.$$

(45)

The formula contains all 9 components of the inverse matrix to $M_3^{\pm }$.

4.4. Thermal Resistances R_1b, R₁₂, R_b and R_a

Final formulas for thermal resistances R_1b, R₁₂, R_b and R_a can now be obtained by again considering the Δ network of Figure 1. For the case of equal pipes in symmetric positions, the thermal resistances R_1b is equal to R_2b, and thus the heat fluxes of Equations (3) and (4) become:

$$q_1=\frac{T_{\mathrm{f}1}-0}{R_{1\mathrm{b}}}+\frac{T_{\mathrm{f}1}-T_{\mathrm{f}2}}{R_{12}},q_2=\frac{T_{\mathrm{f}2}-0}{R_{1\mathrm{b}}}+\frac{T_{\mathrm{f}2}-T_{\mathrm{f}1}}{R_{12}}.$$

(46)

Using Equations (24), (25) and (46), the thermal networks for even (q₂ = q₁) and odd (q₂ = −q₁) cases can be drawn as shown in Figure 4.

As seen from the even case of Figure 4, the borehole thermal resistance R_b between the fluid in pipes and the borehole wall consists of two parallel equal resistances, each of value R_1b. On the other hand, from the odd case of Figure 4, it can be noticed that the total internal resistance R_a between the two pipes consists of a pair of equal series resistances, each of value 0.5 R₁₂. Hence, Equations (5) and (6) can now be simplified to:

$$R_{\mathrm{b}}(J)=\frac{R_{1\mathrm{b}}(J)}{2},R_{\mathrm{a}}(J)=\frac{2R_{1\mathrm{b}}(J)R_{12}(J)}{2R_{1\mathrm{b}}(J)+R_{12}(J)},J=0,1,2,3,\dots .$$

(47)

Final formulas for network resistances R_1b and R₁₂ are obtained in the form of even and odd thermal resistance $R_J^{+}$ and $R_J^{-}$ using Figure 4 and Equations (24) and (25).

$$R_{1\mathrm{b}}(J)=R_J^{+},\frac{1}{R_J^{-}}=\frac{1}{R_{1\mathrm{b}}(J)}+\frac{2}{R_{12}(J)},R_{12}(J)=\frac{2R_J^{+}{R}_J^{-}}{R_J^{+}-R_J^{-}},J=0,1,2,\dots ,$$

(48)

Similarly, final formulas for thermal resistances R_b and R_a are obtained in the form of even and odd thermal resistance $R_J^{+}$ and $R_J^{-}$ using Equations (47) and (48).

$$R_{\mathrm{b}}(J)=\frac{R_J^{+}}{2},R_{\mathrm{a}}(J)=2R_J^{-},J=0,1,2,3,\dots .$$

(49)

For Equations (48) and (49), the values of thermal resistance $R_J^{+}$ and $R_J^{-}$ are obtained from Equations (28) and (38) for the first-order multipoles, from Equations (28) and (42) for the second-order multipoles, and from Equations (28) and (45) for the third-order multipoles.

5. Comparison with Existing Multipole Solutions

In this section, the new higher-order multipole formulas for borehole thermal resistance and total internal thermal resistance are compared with previously published results. The comparison is made using a reference dataset provided by [13], who compared and benchmarked the borehole thermal resistance estimations from ten different analytical methods against the tenth-order multipole method for 216 different cases. The authors showed that compared to other methods, the zeroth-order and first-order multipole formulas provide better accuracies over the entire range of studied parameters. In this paper, we will also benchmark the newly developed second-order and third-order multipole formulas against the tenth-order multipole method. The second-order and third-order multipole formulas will also be compared to each other and to the zeroth-order and first-order multipole formulas to demonstrate improvements in the accuracy of the calculated results.

Table 1. Summary of comparison cases provided by [13].

Parameters	Levels	No. of Levels
Ratio of borehole diameter to outer pipe diameter 1, 2rb/2rp	3, 6, 9	3
Shank spacing configuration, 2xp (mm)	Close, Moderate, Wide 2 For rb/rp = 3, 2xp = 32, 43, 64 For rb/rp = 6, 2xp = 32, 75, 160 For rb/rp = 9, 2xp = 32, 107, 256	3
Ground thermal conductivity, λ (W/m-K)	1, 2, 3, 4	4

¹ Pipe outer diameter (2r_p) is fixed at 32 mm, borehole diameters (2r_b) are 96 mm, 192 mm, and 288 mm. ² Corresponding to Paul’s [8] A–C configurations.

provides the detailed summary of the 216 comparison cases provided by the reference dataset used in this paper for the comparison of the second-order and third-order multipole formulas. The cases cover three different borehole diameters of 96 mm, 192 mm, and 288 mm. The U-tube outer pipe diameter value is held fixed at 32 mm for all cases. The total pipe resistance R_p also remains constant at 0.05 m-K/W. For each borehole diameter, three shank spacing configurations, i.e., close, moderate and wide, corresponding, respectively, to Paul’s [8] Configuration A, Configuration B and Configuration C, are considered. Four levels of ground thermal conductivity ranging from 1–4 W/m-K, and six levels of grout thermal conductivity ranging from 0.6–3.6 W/m-K are used. Given the existing and reasonably foreseeable values of design parameters, the 216 cases used for the comparison bracket almost all real-world single U-tube borehole heat exchangers.

Results of the second-order and third-order multipole formulas for calculating the grout thermal resistance and the total internal thermal resistance are summarized in

Table 2. Mean and maximum absolute percentage errors in calculation of the grout thermal resistance (R_g) for the 216 cases provided by [13].

Method	Shank Spacing Configuration	Ground Conductivity
		Low (0.6–1.2 W/m-K)		Moderate (1.2–2.4 W/m-K)		High (2.4–3.6 W/m-K)
		Mean	Max	Mean	Max	Mean	Max
Zeroth-order multipole	Close	6.1	12.4	2.9	8.0	1.0	2.5
	Moderate	3.3	10.8	1.5	6.5	0.4	1.6
	Wide	8.9	30.4	1.9	11.2	0.3	1.8
First-order multipole	Close	0.2	0.4	0.2	0.6	0.7	1.5
	Moderate	0.0	0.2	0.0	0.2	0.1	0.5
	Wide	0.5	2.2	0.0	0.2	0.1	0.6
Second-order multipole	Close	0.0	0.0	0.1	0.2	0.2	0.5
	Moderate	0.0	0.0	0.0	0.1	0.0	0.1
	Wide	0.1	0.3	0.0	0.0	0.0	0.1
Third-order multipole	Close	0.0	0.0	0.0	0.1	0.1	0.2
	Moderate	0.0	0.0	0.0	0.0	0.0	0.0
	Wide	0.0	0.0	0.0	0.0	0.0	0.0

and

Table 3. Mean and maximum absolute percentage errors in calculation of the total internal thermal resistance (R_a) for the 216 cases provided by [13].

Method	Shank Spacing Configuration	Ground Conductivity
		Low (0.6–1.2 W/m-K)		Moderate (1.2–2.4 W/m-K)		High (2.4–3.6 W/m-K)
		Mean	Max	Mean	Max	Mean	Max
Zeroth-order multipole	Close	23.3	37.6	6.9	15.0	1.2	2.7
	Moderate	1.8	7.8	1.1	6.0	0.3	1.8
	Wide	1.9	8.5	0.4	2.3	0.2	1.3
First-order multipole	Close	3.2	5.9	0.6	0.7	0.7	0.7
	Moderate	0.2	0.8	0.0	0.2	0.1	0.2
	Wide	0.3	1.2	0.0	0.1	0.0	0.1
Second-order multipole	Close	0.4	1.0	0.2	0.2	0.1	0.2
	Moderate	0.0	0.0	0.0	0.0	0.0	0.0
	Wide	0.0	0.2	0.0	0.0	0.0	0.0
Third-order multipole	Close	0.0	0.1	0.1	0.1	0.0	0.0
	Moderate	0.0	0.0	0.0	0.0	0.0	0.0
	Wide	0.0	0.0	0.0	0.0	0.0	0.0

, respectively. Each entry in these two tables represents the mean or maximum error in percentage for a sample containing two three values of grout thermal conductivity, three values of borehole diameter, and four values of ground thermal conductivity. The errors have been determined by comparing the results of multipole formulas to the tenth-order multipole method.

Table 2 shows that the grout thermal resistance (R_g) values obtained from the second-order and third-order multipole formulas are, respectively, within 0.5% and 0.2% of the tenth-order multipole method for all 216 cases. Also, the mean absolute percentage error of the results obtained from the second-order and third-order multipole formula are, respectively, smaller than 0.2% and 0.1%. In comparison, the mean and maximum absolute percentage errors for the zeroth-order multipole formula are as high as 9% and 30%, respectively. The first-order multipole formula has smaller errors than the zeroth-order formula. The mean and maximum absolute percentage errors for the first-order multipole formula are 0.7% and 2.2%, respectively.

Table 3 shows that the total internal thermal resistance (R_a) values calculated from the second order and third-order multipole formulas are, respectively, within 1% and 0.1% of the tenth-order multipole method for all 216 cases. The mean absolute percentage error of the results obtained from the second-order and third-order multipole formulas never exceed 0.4% and 0.1%, respectively. In comparison, the zeroth-order and the first-order multipole formulas give maximum absolute percentage errors of approximately 38% and 6%, respectively. The mean absolute percentage errors of the zeroth-order and the first-order multipole expressions are as high as 23% and 3%, respectively.

As shown and discussed above, compared to the zeroth-order and first-order multipole formulas, the accuracy of the newly developed second-order and third-order multipole formulas to calculate borehole thermal resistance and total internal thermal resistance is higher by several orders of magnitude. Even though the second-order and third-order multipole formulas presented in this paper are more complicated than other analytical expressions including the zeroth-order and first-order formulas, they are still simple enough to apply for computation purposes. The implementation of the second-order multipole formulas involves approximately ten lines of coding of rather compact and simple algebraic expressions. Whereas, the implementation of third-order multipole formulas requires programing of a more complex set of equations involving vector and matrix algebra.

Although, both second-order and third-order multipole formulas offer a significant improvement over the original implementation of the multipole method, which required about 600 lines of FORTRAN coding, but the second-order multipole formulas are intrinsically more straightforward to implement, while providing almost similar levels of accuracy as third-order multipole formulas. Figure 5, Figure 6 and Figure 7 present a selection of representative comparison results of second-order multiple formulas against zeroth-order, first-order and tenth-order multipoles. The results shown in these figures are for a single ground thermal conductivity of 4.0 W/m-K. The left-side figures show the grout thermal resistance R_g values, and the right-side ones show the total internal thermal resistance R_a values, plotted against the grout thermal conductivity. Each figure presents three curves corresponding to close, moderate and wide shank spacings, shown in black, blue and red colors, respectively. The exact value of the shank spacing for each case is provided in Table 1. It must be pointed out that multipole formulas presented in the previous section, calculate the borehole thermal resistance and not the grout thermal resistance. However, in order to be consistent with the dataset provided by [13], the values of grout thermal resistance have been calculated and presented in Figure 5, Figure 6 and Figure 7. The grout thermal resistance values have been determined from Equation (2), by subtracting half of the fixed pipe resistance of 0.05 m-K/W from the corresponding borehole thermal resistance values obtained from the multipole formulas. Computing the grout thermal resistance directly by disregarding the pipe resistance gives erroneous results for all but zeroth-order multipole calculations.

Figure 5, Figure 6 and Figure 7 further demonstrate the ability of second-order multiple formulas to calculate the grout thermal resistance R_g and the total internal thermal resistance R_a. The thermal resistance values calculated from the second-order multiple formulas are always within 1% of the original tenth-order multipole method over the entire range of parameters. Hence, it can be concluded that due to their excellent accuracy and relative ease of implementation, the second-order multipole formulas are recommended for calculation of borehole thermal resistance and total internal thermal resistance for all cases where the two legs of the U-tube are placed symmetrically in the borehole.

6. Conclusions

This paper has presented new second-order and third-order multipole formulas for calculating borehole thermal resistance and total internal thermal resistance. The newly derived formulas can be used for all single U-tube applications where the two legs of the U-tube are symmetrically placed in the borehole. The presented formulas have been compared with the original multipole method, as well as the previously derived zeroth-order and first-order explicit multipole formulas. Both the second-order and third-order multipole formulas provide significant accuracy improvements over the zeroth-order and the first-order multipole formulations. The thermal resistances calculated from the second-order and third-order multipole formulas are always within 1% and 0.2% of the original tenth-order multipole method, respectively. Second-order multipole formulas are, however, more straightforward to use and are, thus, recommended for use in practice.