Explicit Multipole Formula for the Local Thermal Resistance in an Energy Pile—The Line-Source Approximation

: This paper presents a closed-form quite handy formula for the local thermal resistance R b between the temperature of the bulk heat-carrier ﬂuid in the pipes, equally spaced on a concentric circle inside a circular energy pile, and the mean temperature at the periphery of the pile. The so-called multipole method is used to calculate the temperature ﬁeld. An important improvement of the multipole method is presented, where Cauchy’s mean value theorem of analytical functions is used. The formula for thermal resistance R b0 for the zero-order approximation ( J = 0), where only line heat sources at the pipes are used, is presented. The errors using zeroth-order approximation ( J = 0) are shown to be quite small by comparisons with eight-order approximation ( J = 8) with its accuracy of more than eight digits. The relative error for the local thermal resistance R b0 for the zero-order approximation ( J = 0) lies below 5% for a wide range of input parameter values. These ranges are judged to cover most practical cases of application. The smallest local thermal resistance R bmin is, with some exceptions, obtained when the pipes lie directly in contact with the pile periphery. A neat formula for this minimum is presented.


Introduction
Ground-source heat pump (GSHP) systems are among the most energy-efficient and environmentally clean heating and cooling technologies available today [1]. Due to their higher efficiency, they offer substantial potential for energy savings and carbon emissions reduction in buildings [2]. The most common application of GSHP systems is with vertical borehole heat exchangers [3]. Boreholes are routinely drilled to 100-to 400-m depths and require a specific drilling area as well as specialized drilling equipment. The high capital cost of installing borehole heat exchangers is one of the underlying reasons that hamper the widespread deployment of GSHP systems in many places around the world. For new constructions, an alternative to using borehole heat exchangers is to use piled foundations with embedded heat exchange elements as ground heat exchangers. The thermally active foundation piles, commonly referred to as energy or thermal piles, are dual-function structural elements that transmit mechanical loads through unstable soil layers to lower ground levels with high bearing capacity as well as exchange heat with the surrounding ground to provide heating and cooling to the building. They do not require the capital cost and land requirements of borehole heat exchangers and have the potential to lower the environmental footprint of GSHP systems by using fewer materials with low embodied energy.
Energy piles are constructed in many different configurations, as shown in Figure 1, and have typical diameters varying between 0.1 m and up to 3.0 m and lengths ranging from 10 to 60 m [4].
The number and arrangement of heat transfer pipes in the pile depend on the size and construction method of the energy piles. Driven piles of steel or concrete usually have one to two U-pipes arranged in a series or parallel configuration. The heat transfer pipes are placed towards the center of the pile in the hollow driven piles and are positioned around the reinforcement bars in the solid driven piles. Continuous flight auger (CFA) piles have one to three U-pipes installed in the center of the pile. Bored piles have one to several U-pipes mounted on the reinforcement frame close to the boundary of the pile. The differences in construction methods and geometrical arrangements lead to different thermal characteristics of energy piles.
Energies 2020, 13, x FOR PEER REVIEW 2 of 21 The number and arrangement of heat transfer pipes in the pile depend on the size and construction method of the energy piles. Driven piles of steel or concrete usually have one to two U-pipes arranged in a series or parallel configuration. The heat transfer pipes are placed towards the center of the pile in the hollow driven piles and are positioned around the reinforcement bars in the solid driven piles. Continuous flight auger (CFA) piles have one to three U-pipes installed in the center of the pile. Bored piles have one to several U-pipes mounted on the reinforcement frame close to the boundary of the pile. The differences in construction methods and geometrical arrangements lead to different thermal characteristics of energy piles. The current design and analysis methods of energy piles are primarily based on the methods developed for borehole heat exchangers [5]. This is done under the assumption that the thermal characteristics of energy piles and borehole heat exchangers are similar, which is not entirely accurate but simplifies the analysis. Energy piles usually have a greater number of heat transfer pipes and a far smaller aspect (length to diameter) ratio, which differentiates their thermal response from borehole heat exchangers. A key characteristic in this regard is the thermal resistance, which characterizes the heat transfer inside the ground heat exchangers. The thermal resistance Rb of a ground heat exchanger is the ratio of the temperature difference between the circulating fluid in the heat transfer pipes and the surrounding ground to the heat flux per unit length of the ground heat exchanger. The larger number of pipes and their many possible arrangements in energy piles complicate the heat transfer and make the calculation of thermal resistance calculation more complex than for borehole heat exchangers.
There are two main approaches for calculating the thermal resistance of energy piles. In the first approach, thermal resistance is estimated using in situ thermal response testing of energy piles [6,7]. This approach is not generally useful for design purposes as it can only estimate thermal resistance of an existing energy pile. In the second approach, thermal resistance is computed using mathematical models based on conductive heat transfer. Several mathematical models have been proposed and used in this regard. Readers are referred to studies [8][9][10] for a detailed survey and comparison of these models. Some of the commonly used models for calculating the thermal resistance of energy piles include those of Loveridge and Powrie [11], Paul [12], Sharqawy et al. [13], Bauer et al. [14], and Diao  The current design and analysis methods of energy piles are primarily based on the methods developed for borehole heat exchangers [5]. This is done under the assumption that the thermal characteristics of energy piles and borehole heat exchangers are similar, which is not entirely accurate but simplifies the analysis. Energy piles usually have a greater number of heat transfer pipes and a far smaller aspect (length to diameter) ratio, which differentiates their thermal response from borehole heat exchangers. A key characteristic in this regard is the thermal resistance, which characterizes the heat transfer inside the ground heat exchangers. The thermal resistance R b of a ground heat exchanger is the ratio of the temperature difference between the circulating fluid in the heat transfer pipes and the surrounding ground to the heat flux per unit length of the ground heat exchanger. The larger number of pipes and their many possible arrangements in energy piles complicate the heat transfer and make the calculation of thermal resistance calculation more complex than for borehole heat exchangers.
There are two main approaches for calculating the thermal resistance of energy piles. In the first approach, thermal resistance is estimated using in situ thermal response testing of energy piles [6,7]. This approach is not generally useful for design purposes as it can only estimate thermal resistance of an existing energy pile. In the second approach, thermal resistance is computed using mathematical models based on conductive heat transfer. Several mathematical models have been proposed and used in this regard. Readers are referred to studies [8][9][10] for a detailed survey and comparison of these models. Some of the commonly used models for calculating the thermal resistance of energy piles include those of Loveridge and Powrie [11], Paul [12], Sharqawy et al. [13], Bauer et al. [14], The radius of the energy pile is rb. The index b for borehole is kept from other multipole papers. The thermal conductivity in the ground outside the energy pile, r > rb, is λ (W/(m·K)), and it is λb inside the circular energy pile region outside the pipes. The outer radius of the pipes is rp, and Rp denotes the thermal resistance from bulk fluid in the pipes to the outer pipe wall (per meter pipe, (K·m)/W). The pipe resistance Rp is the sum of the thermal resistances of the pipe wall (an annulus) and the fluid boundary layer. The heat injection from each pipe to the ground is q0 (W/m). The value of q0 is negative for heat extraction from the ground. The required fluid temperature Tf for the heat carrier fluid to achieve the prescribed heat flux q0 is to be calculated.
The (steady-state) thermal problem requires that an outer ground temperature is prescribed. As will be discussed in Section 6, it is sufficient to prescribe the average temperature Tbav around the periphery r = rb of the energy pile. This way is close to the thermal problem that was first presented in [19].
The thermal problem has a high degree of symmetry. The temperature field in the wedge indicated in the left-hand Figure 2 around the x-axis with the opening angle φ = 2π/N is repeated around each of the N pipes. The temperature field is symmetric in y, so the upper and lower halves of the wedge mirror each other.
The primary input data are: The thermal resistance Rp ((K·m)/W) may be represented by a nondimensional parameter β, and the ratio between the two conductivities by the parameter σ: The positions of the N pipes on the circle r = rc become using the complex exponential ze: The positions of consecutive pipes are obtained by the complex rotation ze or increase of the polar angle by 2π/N. The last pipe N lies on the positive x-axis at z = rc.
The circles around the pipes and the pile must not intersect. This gives the following restrictions on rp, rc, and rb for any N:  The radius of the energy pile is r b . The index b for borehole is kept from other multipole papers. The thermal conductivity in the ground outside the energy pile, r > r b , is λ (W/(m·K)), and it is λ b inside the circular energy pile region outside the pipes. The outer radius of the pipes is r p , and R p denotes the thermal resistance from bulk fluid in the pipes to the outer pipe wall (per meter pipe, (K·m)/W). The pipe resistance R p is the sum of the thermal resistances of the pipe wall (an annulus) and the fluid boundary layer.
The heat injection from each pipe to the ground is q 0 (W/m). The value of q 0 is negative for heat extraction from the ground. The required fluid temperature T f for the heat carrier fluid to achieve the prescribed heat flux q 0 is to be calculated.
The (steady-state) thermal problem requires that an outer ground temperature is prescribed. As will be discussed in Section 6, it is sufficient to prescribe the average temperature T bav around the periphery r = r b of the energy pile. This way is close to the thermal problem that was first presented in [19].
The thermal problem has a high degree of symmetry. The temperature field in the wedge indicated in the left-hand Figure 2 around the x-axis with the opening angle ϕ = 2π/N is repeated around each of the N pipes. The temperature field is symmetric in y, so the upper and lower halves of the wedge mirror each other.
The primary input data are: The thermal resistance R p ((K·m)/W) may be represented by a nondimensional parameter β, and the ratio between the two conductivities by the parameter σ: The positions of the N pipes on the circle r = r c become using the complex exponential z e : z n = r c · z n e , z e = e 2·π·i N = cos(2 π / N) + i · sin(2 π / N), n = 1, 2, . . . , N; z 1 = r c · z 1 e , z 2 = r c · z 2 e , . . . , z N = r c · z N e = r c ; |z e | = 1, arg(z e ) = 2 π N . (3) The positions of consecutive pipes are obtained by the complex rotation z e or increase of the polar angle by 2π/N. The last pipe N lies on the positive x-axis at z = r c . The circles around the pipes and the pile must not intersect. This gives the following restrictions on r p , r c , and r b for any N: The lower limit is the case when the pipes are squeezed in the center so that consecutive pipes touch each other. In the upper limit, the pipes touch the wall of the energy pile. The left-hand inequalities give an upper limit for N depending on r p /(r b − r p ) as indicated to the right. The pipes may be the two shanks of U-pipes and then N is an even number. However, all formulas in the paper are valid for any N including N = 1.
The thermal resistance R b ((K·m)/W) between the heat carrier fluid in the pipes and the periphery of the energy pile with the average temperature T bav is defined by the relation: The inverse of the thermal resistance is the corresponding thermal conductance K b (W/(K·m)) of the energy pile.

The Multipole Method Applied to the Energy Pile
This section presents the application of the multipole method to an energy pile with multiple U-pipes equally spaced on a circle.

Temperature Field from the N Line Heat Sources
Complex notations and Cartesian and polar coordinates are used: The complex-valued logarithm of any complex function f (z) is defined by: Explicit expressions for the temperature field are derived in Appendices A.1-A.3. The temperature field has the following form: where W c (z, r c ) is obtained from the sum of the N complex-valued line sources in Figure 2, left.
The following compact expressions are derived in detail in Appendices A.1-A.3, Equation (A12): The temperature field involves the absolute value of two polynomials in z: Energies 2020, 13, 5445 6 of 24 The formulas for the temperature in the pile and ground regions become: This is the line-source solution without using multipoles. The above temperature field satisfies the following conditions:

•
Laplace equation in the pile and ground regions, • The heat flux of q 0 (W/m) from each pipe, • The internal boundary conditions of continuous temperature at the pile periphery, • The internal boundary conditions of continuous radial heat flux at the pile periphery (the thermal conductivity is λ b in the pile region and λ in the outside ground), • The average temperature of T bav at the pile periphery.
The temperature field in polar coordinates becomes from (11) and (12): The two denominators in the logarithms become: The temperature around the wall of the energy pile becomes: The wall temperature varies around the average T bav roughly following cos(N·ϕ).

Boundary Condition at the Pipes
The boundary condition at the periphery of any pipe remains to be fulfilled. It is sufficient to consider pipe N at z N = r c due to the rotational symmetry, where the temperature field is repeated for every increase of the polar angle ϕ by 2π/N. See Figure 2 and the discussion at Equation (A13).
Let ρ be the local radial distance from the center of pipe N and ψ the local polar angle from this center z = r c as indicated in Figure 2, right: The circular region around pipe N inside the energy pile area extends at least to the pipe radius r p . Let h p (W/(m 2 ·K)) denote the heat transfer coefficient (per unit pipe area) between the heat carrier Energies 2020, 13, 5445 7 of 24 fluid in the pipes and the outer side of the pipe wall. The pipe thermal conductance K p (W/(m·K)) is obtained by multiplication by the length of the circumference: In the right-hand formula, (2), left, is used. The boundary condition involves a heat balance between the radial heat flux at the outside pipe wall and, the difference between the bulk fluid temperature T f and the temperature outside the pipe multiplied by the heat transfer coefficient h p . The temperature outside the pipe and the heat flux vary around the pipe with the local polar angle ψ: It should be noted that the above two heat fluxes have the dimension of W/m 2 . The above exact boundary condition at each point of the periphery may be written in the following way: The sum of the temperature outside the pipe and the radial heat flux over the pipe divided by h p shall be equal to (a constant) T f for all ψ. The subscript BC is used to denote this left-hand boundary-condition function of ψ. In a more precise solution, the boundary condition may be fulfilled with increasing accuracy by adding multipole components to the temperature field with suitably adjusted strength factors. This is, however, left to a sequel paper. Here, only line heat sources are used, which means that the boundary condition is fulfilled approximately as a mean heat balance for the considered pipe.
The integral of q pN (ψ) around the pipe is equal to the heat flux q 0 . Integration of (20) gives: Using (18), center, the following relation between the fluid temperature T f0 and the mean temperature outside the pipe is obtained in the present approximation without multipoles (J = 0): The subscript 0 is added for the fluid temperature to underline that it is the value for J = 0. The remaining task is now to calculate this mean temperature around the pipe from the temperature field (8) and (9), upper line.

Mean Temperature at the Outer Periphery of the Pipes
The temperature T c (z) consists of the line source at z = r c and a remaining part that is an analytic function of z near z = r c . The direct contributions from the N pipes in (9), upper left, is rewritten in the following way: Energies 2020, 13, 5445 The full solution W c (z, r c ) consists of the line source with its singularity at z = r c and an analytic part W a (z, r c ): The integral of the line source becomes: There is a nice mean-value theorem for analytic functions by Cauchy, which is directly applicable to the integral of W a (z, r c ) around the pipe at z = r c . Let f (z) be a function that is analytic on and inside a circle with the radius r p around a point z 0 in the complex z-plane. Then, the mean value of the function on the circle is equal to the value of f (z) at the center of the circle: This formula may be found in any standard book on analytical functions, for example [22]. An important idea in the improved new approach to the multipole technique presented here is to use this mean-value theorem for the analytic function W a (z, r c ): The integral of T c around the pipe in (22) gets a logarithmic contribution (25) from the line source and another one (27) from the analytical part: The value of W a (z, r c ) for z = r c is from (23) and (24): The fluid temperature now becomes from (22), (27), and (29):

Thermal Resistance of the Energy Pile
The main goal of this paper is to derive a formula for the (local, steady-state) thermal resistance between the heat carrier fluid in the pipes and the periphery of the energy pile with its prescribed average temperature T bav . This resistance becomes from (5) and (30): Energies 2020, 13, 5445 The formula is based on the line-source approximation for the temperature field. An important task is then to determine the accuracy of the formula for various input parameter values; see Section 8.

The Thermal Resistance R b0 as a Function of r c
The thermal resistance R b0 of the energy pile depends on the input parameters (1) and in particular on the radial distance r c from the center z = 0 to the pipes. Equation (31) may be written as: The dependence on r c is expressed by the function R'(r c ): The range of allowed values for r c is limited by the restrictions that the circles of pipes and energy pile must not intersect each other, (4).
An interesting question is now what is the best choice for r c ? The thermal resistance shall be as small as possible. The resistance should normally decrease as the distance to the pile periphery decreases. The derivative of R'(r c ) is readily obtained from (33): The derivative is negative for all r c if σ is negative, i.e., λ b ≤ λ. Then, the minimum occurs at the right end r c = r b − r p . The derivative becomes positive if the right-hand expression within the square brackets is negative. This may happen for larger ratios of λ b /λ and r c close to r b − r p , but for N greater than, say 5, the minimum lies very close to the right-hand limit, and this limit can be used quite safely in the present applications to determine the minimum of R b0 : R bmin = R b0 | r c =r b −r p .
For λ b larger than λ, it may be better to distance the pipes slightly from the energy pile periphery to avoid too much heat flow in the ground region with its lower conductivity. The effect on the minimum is, however, negligible in the present applications. For the cases N = 1, 2, 3, 4, the effect is somewhat greater and worth considering.

Minimum Thermal Resistance for r c = r b − r p
The minimum value for the thermal resistance of the energy pile occurs exactly or with good accuracy when the pipes touch the periphery of the energy pile. The formula for the minimum is obtained from (31): The minimum resistance depends on the ratio r p /r b in the following way: Energies 2020, 13, 5445 10 of 24 The ratio r p /r b is rather small. The following approximation (derived from suitable Taylor expansions of the logarithms) shows the behavior for small ratios: The two quadratic terms in r p may be neglected for say r p ≤ 0.1. The minimum resistance of the energy pile when the pipes touch the pile wall, r c = r b − r p , becomes, with good accuracy:

Formulas and Graphs for the Temperature Field
The temperature field in Cartesian and polar coordinates is given in Section 3.1, Equations (11)- (14). Here, the temperature field is plotted and discussed for a reference case.

Reference Case
The following data are chosen as the reference case: The pipe resistance R p and the conductivity parameter σ become: The average temperature T bav at the pile wall is put to zero in all formulas and graphs below for simplicity. The resistance of the energy pile and the fluid temperature become (31):

Plots of the Temperature Field
The overall temperature field T (x, y) = T pol (r, ϕ) for the reference case (39) is presented in Figure 3, left. Six isotherms, T = 1, 0, −1, −2, −3, −4, are shown. The black circle r = r b shows the periphery of the energy pipe. The isotherm T = 0 winds around the pile periphery along which the average temperature is zero. The pipes with the fluid temperature T f0 = 1.916 in heat carrier fluid lie inside the N = 8 isothermal circles for T = 1. The temperatures at a few particular points become: The isotherm T = −1 has a wavy form influenced by the position of the pipes. The deviation from a circle decreases outwardly. The isotherm T = −3 is very close to a perfect circle. Figure 3, right, shows a surface plot of the temperature field for the reference case. The pipes are clearly seen as temperature spikes. The very flat region around the center z = 0 with values close to T(0,0) = 0.465 is to be noted.
The isotherm T = −1 has a wavy form influenced by the position of the pipes. The deviation from a circle decreases outwardly. The isotherm T = −3 is very close to a perfect circle.   Figure 4 shows radial profiles Tpol (r, φ) from (13) and (14) in the interval 0.2 < r < 0.5 for φ = 0, π/(5N), (3π)/(5N), and π/N. The red top curve for φ = 0 tends to infinity at the position r = rc of the line source. The bottom curve shows the radial profile for the radial line φ = π/N, which lies in the middle between pipe n = N and pipe n = 1 in Figure 2. The black dotted line shows the radial average temperature Tpolav (r), which is defined below in (46). The four curves converge to the radial average temperature for r greater than 0.4. The difference between the top and bottom curve is 0.091, 0.015, 0.004, and 0.00006 at r = 0.4, 0.5, 0.6, and 1, respectively.   (13) and (14) in the interval 0.2 < r < 0.5 for ϕ = 0, π/(5N), (3π)/(5N), and π/N. The red top curve for ϕ = 0 tends to infinity at the position r = r c of the line source. The bottom curve shows the radial profile for the radial line ϕ = π/N, which lies in the middle between pipe n = N and pipe n = 1 in Figure 2. The black dotted line shows the radial average temperature T polav (r), which is defined below in (46). The four curves converge to the radial average temperature for r greater than 0.4. The difference between the top and bottom curve is 0.091, 0.015, 0.004, and 0.00006 at r = 0.4, 0.5, 0.6, and 1, respectively.  The four curves lie very close to the value T (0, 0) = 0.465 at the center for 0 < r < 0.2. The deviation from temperature at r = 0 is less than 0.0002 for r < 0.1 and less than 0.06 for r < 0.2. The temperature field is very flat in the center region of Figure 3 (left) and 3 (right). This means that changes in the thermal conductivity or introduction of a return pipe with the same fluid temperature in the central region r < 0.2 will have a negligible, or perhaps marginal, effect on the value of the thermal resistance.

Total Radial Heat Flux: Thermal Resistance to Any Radius r > rb
The Laplace equation in polar coordinates for a thermal conductivity that depends on r reads: The equation is valid everywhere except at the N heat sources on the circle r = rc. Let Tpolav (r) denote  The four curves lie very close to the value T (0, 0) = 0.465 at the center for 0 < r < 0.2. The deviation from temperature at r = 0 is less than 0.0002 for r < 0.1 and less than 0.06 for r < 0.2. The temperature field is very flat in the center region of Figure 3 (left) and 3 (right). This means that changes in the thermal conductivity or introduction of a return pipe with the same fluid temperature in the central region r < 0.2 will have a negligible, or perhaps marginal, effect on the value of the thermal resistance.

Total Radial Heat Flux: Thermal Resistance to Any Radius r > r b
The Laplace equation in polar coordinates for a thermal conductivity that depends on r reads: The equation is valid everywhere except at the N heat sources on the circle r = r c . Let T polav (r) denote the average temperature obtained from integration over −π < ϕ < π. Integration of (43) in ϕ gives: The heat flux in the ϕ-direction is the same at ϕ = π and ϕ = −π. The right-hand term above vanishes and the average temperature T polav (r) satisfies the one-dimensional radial heat conduction equation except at r = r c , where the total radial heat flux q totrad (r) increases from zero to N·q 0 : The temperature T polav (r) becomes constant for r < r c . It is zero at r = r b (for T bav = 0) and it varies logarithmically as ln (r b /r) for r > r c , inversely proportional to the thermal conductivities λ b and λ: This average temperature is shown by the dotted line in Figure 4. It should be noted that the derivative of the average temperature is discontinuous at r = r b following λ (r), (43) right, while the heat flux is continuous. It should also be noted that T polav (0) from (46), right, gives the same value T (0) = 0.465 as in (42).
The thermal resistance R b refers to the resistance between the fluid in the pipes and the average temperature at the periphery of the energy pile, (5). The corresponding thermal resistance between the fluid in the pipes and the average temperature at any radius r 0 > r b is obtained in the following way: The resistance from fluid to r = r 0 becomes equal to the sum of R b and the thermal resistance of the annular region r b < r < r 0 with its thermal conductivity λ: This simple thermal network is illustrated in Figure 5.
Energies 2020, 13, x FOR PEER REVIEW 12 of 21 The thermal resistance Rb refers to the resistance between the fluid in the pipes and the average temperature at the periphery of the energy pile, (5). The corresponding thermal resistance between the fluid in the pipes and the average temperature at any radius r0 > rb is obtained in the following way: The resistance from fluid to r = r0 becomes equal to the sum of Rb and the thermal resistance of the annular region rb < r < r0 with its thermal conductivity λ: This simple thermal network is illustrated in Figure 5.

Deviation from the Exact Boundary Condition at the Pipes
The boundary condition at the pipes is discussed in Section 3.2. Equation (20) gives the exact boundary condition at the periphery of the pipes:

Deviation from the Exact Boundary Condition at the Pipes
The boundary condition at the pipes is discussed in Section 3.2. Equation (20) gives the exact boundary condition at the periphery of the pipes: The temperature at the outside of the pipe as a function of z is given by (8) and (9), top line. The radial heat flux q ρN (ψ) is defined in (19), right. It is fairly straight-forward to determine this radial derivative with respect to ρ from the formula for the temperature: This gives an explicit formula for T BC (ψ). Figure 6 shows the function (49) for the exact boundary condition around the periphery of the pipes for the reference case (39). The chosen fluid temperature T f0 is the average or mean value of the boundary function over ψ as derived above in Section 3.2, Equations (28) and (30). The difference between the two curves represents the deviation from the exact boundary condition. The largest difference 0.092 between T BC (ψ) and the chosen T f0 occurs at ψ = ±π. The relative error lies below 5%. The temperature at the outside of the pipe and the radial heat flux q pN , (50), divided by the heat transfer coefficient h p are also shown. These curves vary more than T BC (ψ), but they counteract each other so that the deviation in the boundary condition is much smaller.

Comparison with Exact Results from Multipoles Solutions of Order J = 8
Higher-order multipole solutions (J = 1, 2, …) give the thermal resistance with increasing accuracy. The general solution for any order of multipoles J will be reported in a sequel paper. Multipoles up to J = 8 are used in the comparisons of accuracy for the thermal resistance of the energy pile.

Error for Rb (J) as a Function of J
The error for Rb (J) compared to Rb (8) is:

Comparison with Exact Results from Multipoles Solutions of Order J = 8
Higher-order multipole solutions (J = 1, 2, . . . ) give the thermal resistance with increasing accuracy. The general solution for any order of multipoles J will be reported in a sequel paper. Multipoles up to J = 8 are used in the comparisons of accuracy for the thermal resistance of the energy pile.

Error for R b (J) as a Function of J
The error for R b (J) compared to R b (8) is: Energies 2020, 13, 5445 14 of 24 Figure 7 shows the absolute relative error |∆R b (J)| for J = 0, 1, . . . 7 for the reference case.
to J = 8 are used in the comparisons of accuracy for the thermal resistance of the energy pile.

Error for Rb (J) as a Function of J
The error for Rb (J) compared to Rb (8) is: Figure 7 shows the absolute relative error |∆Rb (J)| for J = 0, 1, … 7 for the reference case. The relative error for J = 0 is 0.0069 or 0.69%. There is a dramatic increase in accuracy with increasing J. This is typical for the multipole method. The multipole method gives the solution and in particular the thermal resistance Rb (8) with 8 digits of accuracy. The error for thermal resistance is always much smaller than the relative deviations for the boundary conditions at the pipes; see Figure 6. ∆Rb0 (N, rb, rc, rp, σ, β)

for a Suitable Set of Parameter Values
The input parameters (1) determine the temperature field and the thermal resistance of the energy pile. The temperature field (11) and (12) is proportional to q0/λb. The non-dimensional relative The relative error for J = 0 is 0.0069 or 0.69%. There is a dramatic increase in accuracy with increasing J. This is typical for the multipole method. The multipole method gives the solution and in particular the thermal resistance R b (8) with 8 digits of accuracy. The error for thermal resistance is always much smaller than the relative deviations for the boundary conditions at the pipes; see Figure 6.

Relative Error ∆R b0 (N, r b , r c , r p , σ, β) for a Suitable Set of Parameter Values
The input parameters (1) determine the temperature field and the thermal resistance of the energy pile. The temperature field (11) and (12) is proportional to q 0 /λ b . The non-dimensional relative error ∆R b0 = ∆R b (0) with J = 0 in (51) will depend on the number of pipes N, two dimensionless parameters, and the three radii: It may be noted that the relative error is unchanged when the three radii are scaled with the same factor k > 0: ∆R b0 (N, σ, β, r p , r c , r b ) = ∆R b0 (N, σ, β, k · r p , k · r c , k · r b ).
To investigate the performance of the proposed formula, a parametric study is made, varying the number of pipes N, the pile and ground thermal conductivities, the dimensionless pipe thermal resistance, and the pile geometry. The parameters and levels for the full parametric study, representing 6 × 5 × 3 × 3 × 4 = 1080 cases, are summarized in (54): N = 2, 4, 6, 8, 10, 12; 2 · r b = 0.16, 0.3, 0.6, 1.2, 2.4; r c = r c,min , 2 r b 3 , r b − r p ; 2 · r p = 0.032; λ b / λ = 0.5, 1, 2, or σ = − 1 3 , 0, 1 3 ; β = 0.25, 0.5, 1, 2. (54) These first 3 parameters represent 90 different geometries; in all cases, the outer diameter of the U-pipe is 0.032 m. The pile diameters are between 0.16 and 2.4 m. The smallest diameter is smaller than a typical precast driven thermal pile and the largest diameter is closer to what might be found in a large CFA or bored thermal pile. Taken together, the considered diameters bracket almost all energy pile installations. The three r c radii correspond to close, moderate, and wide pipe spacings. The pipes touch each other and the wall of the energy pile for r c equal to r c,min (close spacing), and r c equal to r b − r p (wide spacing), respectively. These cases are typical of CFA and hollow driven piles, respectively. The case r c equal to 2r b /3 (moderate spacing) is representative of a bored energy pile where heat exchanger pipes are mounted on the reinforcement ring. Again, the three spacings virtually bracket all possible spacings. Similarly, the three levels of sigma bracket the most likely ratios of the borehole and ground thermal conductivities found in practice. Table 1 shows the relative errors ∆R b0 (N, r b , r c , r p , σ, β) in percent. The grey cells in the table indicate relative errors of 10% or more. More than 50% of the errors are marked as zero, which means that the error is below 0.5%. For r b greater than or equal to 0.3 m, the relative error ∆R b0 (N, r b , r c , r p , σ, β) is less than or equal to 5% for 643 of the 648 cases and less than 10% for all cases. For r b = 0.15 m, the relative error ∆R b0 (N, r b , r c , r p , σ, β) is less than or equal to 5% for 156 of the 216 cases and less than 10% for all except a few cases. The relative error exceeds 10% only for 8 cases when N is greater than or equal to 8 and β is equal to 0.25 or 0.50 (i.e., for very low R p and R b values). These cases are very unlikely in practice and are possible only under the implausible conditions of extremely low thermal conductivity of energy piles and laminar flow in the heat exchanger pipes. Moreover, although the relative error in these cases is large, the absolute error is small and of the order of 0.002 (K·m)/W) or less. For r b equal to 0.08 m, the relative error ∆R b0 (N, r b , r c , r p , σ, β) is greater than 10% for over one-third of the considered cases. With N equal to 12 and r p equal to 0.016 m, the r c,min becomes larger than 2r b /3, and hence the moderate spacing condition does not hold. The relative error for these cases with r b equal to 0.08 m is classified as not applicable (n/a) in Table 1. An interesting observation from Table 1 is that the relative error ∆R b0 (N, r b , r c , r p , σ, β) is positive for β less than 1 and is negative for β larger than 1. This suggests that the zero-order formula overestimates the thermal resistance of energy piles for β less than 1 and underestimates the thermal resistance for β greater than 1. For all cases with β equal to 1, the relative error is fairly close to zero. This information is valuable when setting margins of safety for design.
One parameter that can affect the accuracy of the multipole Formula (39) is the pipe radius. To demonstrate the impact of different pipe radii on the relative errors of the parametric study (54), Tables 2 and 3 present relative errors ∆R b0 (N, r b , r c , r p , σ, β) for r b equal to 150 mm and r p equal to 12.5 mm, and r b equal to 300 mm and r p equal to 20 mm, respectively. Comparing the results of Tables 1 and 2 indicates that for energy piles with r b equal to 0.150 m, the relative errors ∆R b0 (N, r b , r c , r p , σ, β) become considerably lower for r p equal to 0.0125 m than for r p equal to 0.016 mm. For r p equal to 0.0125 m, the relative error exceeds 10% only when N is greater than or equal to 10, and β equal to 0.25. Comparison of the results of Tables 1 and 3 shows that for energy piles with r b equal to 0.3 m, the relative errors ∆R b0 (N, r b , r c , r p , σ, β) are marginally higher for r p equal to 0.02 m than for r p equal to 0.016 m. Nevertheless, the relative error for r p equal to 0.02 m exceeds 10% only for the single case of N equal to 12, σ equal to −0.33, and β equal to 0.25. Table 1. Relative error ∆R b0 (N, r b , r c , r p , σ, β) in percent for N = 2 to 12, and r c = r c,min (close spacing), 2r b /3 (moderate spacing), and r b − r p (wide spacing) a . The pipe radius r p is 0.016 mm.
a Relative error larger than 10% shown with a grey background. Table 3. Relative error ∆R b0 (N, r b , r c , r p , σ, β) in percent for N = 2 to 12, and r c = r c,min (close spacing), 2r b /3 (moderate spacing), and r b − r p (wide spacing) a . The pipe radius r p is 0.02 m.

Further Discussion and Concluding Summary
The newly developed zeroth-order multipole formula provides a fast, flexible, and accurate method for calculating the thermal resistance of energy piles. The newly presented zero-order multipole Formula (31) for R b0 may be used in virtually all engineering applications for energy piles with r b greater than or equal to 0.3 m and r p less than or equal to 0.02 m with pipes equally spaced on a circle. The relative error ∆R b0 (N, r b , r c , r p , σ, β) may exceed 10% in a few odd cases, but the absolute error in such cases is very small. Furthermore, the zero-order multipole Formula (31) may also be used for all energy piles with r b between 0.15 and 0.3 m and r p less than or equal to 0.0125 m if N is less than or equal to 10, and β is greater than or equal to 0.25. For applications with r p equal to 16 mm and r b equal to 0.15 to 0.3 m, the zero-order multipole Formula (31) may be used if N is less than or equal to 8. The zero-order multipole formula may also be used for all energy pile applications if β is equal to 1 or N is equal to 2.
To the best of the authors' knowledge, zeroth-order multipole Formula (31) is the first, and the only closed-form formula published to date, that can be used to calculate the thermal resistance of energy piles with any number of heat exchanger pipes. The proposed formula marks a substantial improvement on the current practice of sizing energy piles with three or more U-pipes based on a guessed value of thermal resistance, which can lead to improperly sized systems. The explicit formula will facilitate designers to step beyond the limits of using a conservative thermal resistance estimate, and will enable them to optimize the energy pile design by making informed decisions about design choices, such as the position and the number of heat exchanger pipes in the pile. The proper sizing of energy piles will, in turn, allow for the overall system size to be reduced and the system performance to be improved.
A possible direction for future research is the development of higher-order multipole formulas for calculating the local thermal resistance of energy piles. The higher-order formulas will further improve the accuracy of the thermal resistance calculation, especially for small-diameter energy piles (less than 300 mm). The derivation of such formulas is expected in the near future and they will be presented soon in a sequel paper.
To conclude, in this paper, a closed-form Formula (31) for the thermal resistance of energy piles with any number of heat exchanger pipes evenly spaced on a concentric circle inside the pile periphery was derived. The proposed formula is based on the zero-order multipole approximation (J = 0) and considers only line heat sources at the pipes. It was demonstrated that the newly derived formula can calculate the thermal resistance of large-diameter energy piles (r b ≥ 0.3 m) with an accuracy better than 10% essentially under all conditions. The relative errors are less than or equal to 5% for over 99% of the 648 cases analyzed in this study. It was also demonstrated that the zero-order formula can also be used for smaller-diameter energy piles (0.15 ≤ r b < 0.3 m) for most conditions, except when the dimensionless pipe thermal resistance β is very small and the number N of heat exchanger pipes in the energy pile tops 8. For these specific cases and for energy piles with particularly smaller diameters (r b < 0.15 m), the use of higher-order multipoles is recommended-closed-form formulas for which will be presented in our forthcoming papers. This paper also presents a handy Formula (35) for the minimum thermal resistance of energy piles when the pipes touch the periphery of the pile.