## Abstract

Slender-body theory (SBT) for transient heat transfer from bodies whose lengths are much larger than their radius into a conductive medium is derived. SBT uses matched asymptotic expansions of inner and outer solutions. An analytical inner solution for heat transfer from a circular cross section is matched to an outer solution obtained using Green’s functions. An efficient numerical implementation is obtained based on a judicious choice of the discrete elements used to represent the body and implementation of the fast multipole method (FMM). The SBT requires a one-dimensional spatial discretization only along the axis of the body in contrast to the three-dimensional discretization for finite-element models. Two case studies, heat transfer from two parallel cylinders and heat transfer from a slinky-coil heat exchanger, are used to show the speed and accuracy of the SBT model and its ability to model interacting slender bodies of finite length and bodies with centreline curvature and internal advective heat flow.

## 1. Introduction

Simulating and predicting the heat exchange of slender bodies undergoing transient heat transfer in conductive media is of great interest and importance in the design and analysis of several applications. Examples are (i) wellbore heat transmission in oil, gas or geothermal wells, (ii) reservoir heat transfer in ground-source heat pumps or underground thermal energy storage systems using slinky-coil, horizontal pipe or borehole heat exchangers and (iii) cooling or heating elements in concrete foundations and industrial processes.

Various heat transfer solution methods are available and in some cases sufficient, depending on the complexity and compactness of the slender body heat transfer problem. Simple geometries, e.g. a straight finite cylinder, can be handled with analytical methods such as the Laplace transform, separation of variables and Bessel functions, Green’s functions and perturbation methods [1–6]. More complex but compact geometries, meaning the slender body is, e.g. irregularly curved but closely packed together, can be simulated with finite-element (FEM) or finite-volume methods [7,8]. Nevertheless, using these numerical methods for non-compact geometries typically requires long computational times in order to achieve accurate results. In this paper, we present a hybrid method that efficiently handles these type of geometries. The method uses the *slender-body theory* (SBT), an asymptotic analysis to allow for thin non-compact geometries, in combination with a numerical implementation to allow these geometries to have complex shapes. Hereby, one can reduce from discretizing a three-dimensional domain to discretizing over a one-dimensional axis of the slender body. Especially for longer time scales, when the domain needed to be studied becomes *L*_{t} the total slender body length, much shorter computational times can be obtained and much less computer memory would be needed. However, without smart implementation strategies, the computational time is *N*_{E} and *N*_{t} the number of elements and time steps, respectively. Also, the method requires a linear problem, so no temperature dependence of the thermal conductivity is allowed. Further, simple medium geometries are desired, so that the Green’s function is either unbounded or can be readily found, e.g. with method of images.

SBT, in general, refers to approximating the equations or underlying physics in a specific science field when dealing with bodies that are long in one dimension and short in the other two. These approximations are typically achieved by applying matched asymptotic expansions of inner and outer solutions, although sometimes only the inner or outer solution is considered. It originated and is most commonly used in the field of fluid dynamics. Munk [9] in 1924 was the first to apply slender-body approximations in order to calculate the aerodynamic forces on airship hulls. He treated the airship as an *elongated body of revolution* in a potential flow and argued that when calculating the velocity potential and air forces, different longitudinal cross sections can be considered independent from each other. In 1927, von Kármán [10] applied a different SBT approach in solving the same problem using a distribution of monopoles and dipoles on the axis of symmetry. Other examples of potential flow problems simplified with SBT are the study of wings [11], not-so-slender wings [12], ship hydrodynamics [13,14] and fish locomotion [15,16]. Slender bodies in Stokes flow have also been widely studied using SBT methods. Examples are the general works by Cox [17,18], Keller & Rubinow [19], Batchelor [20] and Johnson [21] and the studies on flagellar hydrodynamics [22], fibre particle suspensions [23,24], and motion of a slender torus [25]. Other applications of SBT can be found in the fields of acoustic scattering [26,27], electromagnetic scattering [28], charged colloids [29] and steady-state heat transfer. A steady-state heat transfer problem often approximated using SBT is the heat transport in composites with slender inclusions [30–36], e.g. to derive the effective thermal conductivity. To the best of our knowledge, SBT has never been applied to transient heat conduction.

In §2, the slender body geometry, constraints and heat transfer problem are formally defined. A theoretical framework for solving the problem using Green’s functions and matched asymptotic expansions is derived, and the result is physically interpreted. Further, model extensions are discussed to treat various time scale, medium properties and initial and boundary conditions. In §3, the theoretical basis is translated into a numerical method and strategies are provided for efficient implementation in order to obtain high speed, accuracy and versatility. Finally, these strengths are demonstrated in §4 using two case studies in which the presented SBT numerical model is compared with an FEM. The first case study is the simulation of the heat exchange of two parallel finite cylinders at constant temperature; the second case study is the simulation of one loop of a slinky-coil heat exchanger.

## 2. Theoretical basis

### (a) Slender body geometry, constraints and heat transfer problem

Figure 1 is a sketch of the type of slender-body transient heat conduction problem studied in this work. The geometry of the slender body is set by the functions ** F**(

*s*) and

*R*(

*s*), which define, as a function of the arc-length

*s*along the centreline, the location of the centreline and the radius in parametric form in the global (

*xyz*)-coordinate system. A local (

*x*′

*y*′

*z*′)-coordinate system is also defined with

*z*′ parallel to

*s*. The total length of the slender body is

*L*

_{t}and its radius of curvature is represented by the function

*R*

_{C}(

*s*). The slender body is situated in an infinite medium with constant thermal diffusivity

*α*and constant thermal conductivity

*k*. For SBT to be applicable, the pipe radius must be small compared with its length, radius of curvature and the thermal diffusion length:

*τ*a representative time scale. Further, the pipe radius must be small compared with the distance to a nearby slender body. In §2d, a finite medium will be considered with, then the additional constraint that the pipe radius is small compared with the distance from boundaries. Also, in §2c, the model will be expanded for handling thermal diffusion lengths on the order of the pipe radius or smaller and hereby lifting the constraint posed by equation (2.3). The governing equation for transient heat conduction in the medium is

*T*

_{m}(

*x*,

*y*,

*z*,

*t*) the temperature in the medium and

*t*the time. It is assumed that the initial medium temperature is zero and the body is in perfect thermal contact with the medium.

We aim to find a relation between the temperature *T* at the surface of the slender body and the heat exchange *Q* between the slender body and the conductive medium. For the constraints given above, *T* and *Q* will show no angular dependence in the matching region to a first approximation, because higher-order multipoles will decay away, and therefore, will only be a function of *s* and *t*. The linear relationship sought between *T* and *Q* can now be written in its most general form as
*T* or *Q* is set and the other one is calculated for using equation (2.5). In more complex problems with heat transfer processes occurring *inside* the slender body, e.g. internal advective heat flow, equation (2.5) must be combined with other equations, as discussed in §2d.

### (b) Slender-body theory derivation for transient heat conduction

The idea of SBT is to develop an inner and outer solution for the problem, which are valid nearby and faraway from the body, respectively, and then impose their compatibility in the matching region. For the inner solution in our case, the temperature disturbance locally is assumed independent from nearby sections, i.e. the body is replaced with an infinite cylinder with no axial or angular heat variation and the temperature calculation simplifies to a two-dimensional axisymmetric heat transfer problem. For the outer solution, the slender body is replaced with a distribution of point sources along its centreline, and the medium temperature calculation becomes a convolution over *s* and *t* of this heat distribution with the Green’s function. This approach is equivalent to SBT for Stokes flow, where in the inner region, the local velocity has no axial dependence, whereas in the outer region, the slender body can be replaced to a first approximation with a distribution of Stokeslets [20].

#### (i) Inner solution

To estimate the temperature of the conductive medium in close proximity to the slender body, i.e. *A* at the slender body (figure 1) at a certain time *t*, the slender body is at temperature *T*(*s*_{A},*t*) and exchanges heat *Q*(*s*_{A},*t*). The medium temperature in close proximity to *A*, defined as the inner solution *T*_{in}, is estimated in the local (*r*′*z*′)-coordinate system as [1]
*k* the medium thermal conductivity, *ln*(*x*) the natural logarithm, and *R*(*s*_{A}) the local radius of the slender body. Equation (2.6) cannot be viewed as a solution of the inner region problem independent from the outer solution, because we cannot specify both *T* and *Q*. Instead, we need to match this equation to an outer solution to obtain the relationship between *T* and *Q*.

#### (ii) Outer solution

To estimate the medium temperature at time *t* at an arbitrary point *B* far away from the slender body, defined as the outer solution *T*_{out}, the slender body is replaced by a distribution of transient point sources with strength *Q*(*s*,*t*) along its centreline. Using Green’s functions, *T*_{out} is then calculated at *B* in the global (*xyz*)-coordinate system as [1]
** b**−

**(**

*F**s*) the distance vector between

*B*and every point along the centreline. Equation (2.7) represents the convolution of the heat distribution

*Q*(

*s*,

*t*) with the point source solution

*G*(

**,**

*d**τ*) (Green’s function) calculated as

*ρ*,

*c*

_{p}and

*α*the density, specific heat capacity and thermal diffusivity of the conductive medium, respectively and |

**| the length of the distance vector**

*d***.**

*d*When matching, the outer solution will be evaluated for the general point *B* approaching asymptotically the centreline of the body at an arbitrary location, e.g. point *A* at *s*=*s*_{A}. However, with the double integral in its current form, a singularity arises owing to dividing by zero at point *A* for *t**=*t*. The resolution of this singular lies in locally, around *A*, simplifying the slender body as a straight pipe (for which an analytical solution exists) and adding and subtracting the resulting straightline source with constant heat exchange *Q*(*s*_{A},*t*) to the right-hand side (RHS) of equation (2.7), an approach comparable to SBT for steady-state heat transfer [34]. The starting time of the line source, *t*_{0}, can be any time between 0 and *t* but the end time must be *t*. The straight line has a length 2*l* much smaller than the local radius of curvature (2*l*≪*R*_{C}(*s*_{A})) owing to geometry reasons, but much larger than the local radius (2*l*≫*R*(*s*_{A})), which will be beneficial for the short time-scale model extension (§2c). Equation (2.7) becomes
*H* the Heaviside step function. In equation (2.9), the singularity has been removed, because, at point *A* for *t**=*t*, *Q*(*s*,*t**)=*Q*(*s*_{A},*t*) and therefore, the argument in the double integral in term *M* is 0. Also, keeping in mind that the slender body is locally represented as a straightline section, term *K* can be calculated analytically which will produce an ln(*r*′) term that can be matched with the ln(*r*′) term of the inner solution. In preparation for matching, the location of point *B* is expressed in the local (*r*′*z*′)-coordinate system attached to point *A*. Term *K* is then calculated as
*B* asymptotically approaches point *A*. Term *K* will be evaluated for *z*′=0 and can therefore be rewritten as
*ϕ*=*l*^{2}/(*α*(*t*−*t*_{0})), *E*_{1}(*x*) the exponential integral function *erf*(*x*) the error function *x*) the complimentary error function (1−erf(*x*)) [37].

#### (iii) Matching

We now match the inner and outer solution, after adding and subtracting *Q*(*s*_{A},*t**) in term *M* of equation (2.9) to split *M* into two parts, and with further taking into account that *x* with *γ* Euler’s constant (0.5772…):
*Q*(*s*_{A},*t*) during the time period *t*_{0} to *t* and over the total length 2*l*. Careful analysis reveals this term consists of two parts. The first part, in the form *E*_{1}(*x*^{2}) for small *x* and represents the temperature change owing to an infinite line source [1]. The second part, the integral term, accounts for the finite length (2*l*) by subtracting the temperature change owing to two semi-infinite line sources, one on either end. The second and third term in equation (2.13) can be viewed as corrections to the first term by now accounting for the local heat source variation in space and time. The second term is the temperature contribution from the same FLS but now with a heat pulse strength that is still constant along the line source but varies over time, *Q*(*s*_{A},*t**), subtracted by *Q*(*s*_{A},*t*) for the time period *t*_{0} to *t*. In the third term, the line source has the actual heat source *Q*(*s*,*t**) varying over space and time, subtracted by *Q*(*s*_{A},*t**). Finally, the fourth term in equation (2.13) is the temperature contribution from all other previous and current point sources along the centreline of the slender body.

### (c) Extension of model to short time scales

Equation (2.13) is only valid for heat transfer processes occurring over long time scales, i.e. *R* and constant heat pulse *Q* from *t*_{0} to *t* is given by [1]
*J*_{1}(*x*) and *Y* _{1}(*x*) the first-order Bessel function of the first kind and second kind, respectively. Equation (2.14) readily replaces the FLS model in the first term of equation (2.13). The effect of *Q*(*s*_{A},*t*) changing over time, represented by the second term in equation (2.13), should also be calculated with the FCS model. In continuous form, this can only be written analytically in the Laplace domain. However, a solution can be found in the time domain when discretizing *Q* as a series of *n* constant heat pulses (as will be applied when numerically implementing),
*Q*_{n}(*s*_{A}) the constant heat pulse at time step *n*. No FCS model is available when *Q*(*s*,*t*) varies spatially along the cylindrical source, represented by the third term in equation (2.13). However, when considering 2*l* to be the length of a computational element on which *Q* is constant (as will be applied in §3), allowed in the case of spatially smooth variations in *Q* along the centreline of the body, this term can be neglected. Finally, the fourth term remains unaltered, because it accounts for temperature contributions caused by heat sources at a distance larger than *R*(*s*_{A}) away for the point of interest (*A*). All together, the solution for *T*(*s*_{A},*t*) applicable at short time scales becomes
*t*_{0} now corresponds to one of the time steps. Equation (2.16) is also valid for long time scales, however, from a computational point of view, one rather wants to calculate equation (2.13), as will be further discussed in §3.

### (d) Other model extensions

The equations developed so far assume an infinite, conductive medium with uniform zero initial temperature, and constant, isotropic and uniform thermal properties. This section discusses how to remove several of these constraints, as well as how to incorporate heat transfer processes inside the slender body.

The effect of a far-field or initial temperature distribution can be accounted for by modelling it as an imposed field *T*_{i}(*x*,*y*,*z*,*t*) without the slender body and then superimposing this field to the disturbance field caused by the slender body, i.e. adding *T*_{i}(*x*,*y*,*z*,*t*) to the RHS of equation (2.9) and therefore to the RHS of equation (2.13) or (2.16). In the most simple cases, the imposed field is found directly such as for a uniform non-zero initial temperature, *T*_{i}(*x*,*y*,*z*,*t*)=*T*_{0}, or for a linear far-field temperature profile, *T*_{i}(*x*,*y*,*z*,*t*)=*T*_{0}(*z*)=*az*+*a*_{0} (e.g. geothermal gradient).

For a non-infinite medium, the infinite medium model can still be an accurate representation in cases when the distance from the slender body to the medium boundary is much larger than the thermal diffusion length. An example is the simulation of the heat transmission with deep geothermal, oil or gas wells. Otherwise, a different (and more onerous) Green’s function must be used with a complete list of Green’s functions for all types of boundary conditions provided by Beck *et al.* [5]. Some can be readily found by using method of images, such as constant temperature or zero heat flux condition at a planar boundary [5,38]. Several of the solution methods provided by Beck *et al.* can be derived using the imposed field technique, introduced above for dealing with the initial condition. For example, a distribution of heat sources at a boundary must be convoluted over space and time with the Green’s function for an infinite medium, and then superimposed to the disturbance field caused by the slender body, calculated with SBT using a zero heat flux Green’s function.

Simple cases of a medium involving fluid flow, or anisotropic or non-uniform thermal properties can be simulated by using a different (and again more complicated) Green’s function. One-dimensional constant velocity fluid flow causing heat advection, instead of pure heat conduction, can be accounted for by using the moving point source Green’s function [1,39]. However, this is only directly applicable for the long time-scale solution, and should be combined with a different model, e.g. the solid cylindrical source model [40], for the short time-scale solution. Also a Green’s function exists in the case of anisotropic thermal properties, for which the thermal conductivity can now be represented by a tensor [41]. Further, for a two-layered medium with isotropic properties, an exact Green’s function is provided by Carslaw & Jaeger [1]. For any number of layers, Abdelaziz *et al.* [42] provide an approximate method using interlayer heat exchange adjustments, and laws of composite materials for estimating representative thermal properties. Again, the extensions for anisotropic and non-uniform thermal properties are only directly applicable at long time scales when the cylindrical source models can be represented by line source models. Incorporating temperature-dependent thermal properties with the SBT model is not possible since this would lead to nonlinear equations.

In simple cases, either *T* or *Q* at the surface of the slender body is specified with the other being solved for. In more complex problems with heat processes occurring inside the body, additional relations are required and calculated for simultaneously with equation (2.13) or (2.16). For example, for advective heat transfer owing to pipe fluid flow, a first additional relation is [43]
*T*_{f}(*s*,*t*) the temperature of the fluid and *R*_{t}(*s*) the thermal resistance of the pipe calculated as the sum of the thermal resistances of the conductive heat transfer through the pipe and the convective heat transfer between the fluid and the pipe wall [43]:
*R*(*s*) and *R*_{inner}(*s*) the outer and inner radius of the pipe, respectively. Further, *k*_{p}(*s*) is the thermal conductivity of the pipe material and *h*(*s*) is the convective heat transfer coefficient which can be calculated with appropriate Nusselt correlations for forced internal convection [43]. A third equation is required, which represents the energy balance of the fluid inside the pipe (with neglecting the heat conduction within the fluid)
*ρ*_{f}, *c*_{p,f} and *U*(*s*) the density, specific heat capacity and velocity of the fluid, respectively, and the fluid inlet temperature being a boundary condition. Various finite difference methods are available to handle the time and spatial derivative, with different levels of accuracy and complexity [7]. For the second case study for example (§4b), equation (2.19) is numerically implemented using an upwind differencing scheme. For other types of heat exchangers such as single-U (1U), double-U (2U) and coaxial borehole heat exchangers, similar equations can be developed.

Finally, one can derive from equation (2.13) the SBT solution for steady-state heat conduction by replacing *Q*(*s*_{A},*t*) with *Q*(*s*_{A}), and letting time approach infinity. The time convolution integral of the transient Green’s function (equation (2.8)) for any starting time *t*_{0} translates into the steady-state Green’s function [5]:
*T*(*s*_{A}) the surface temperature of the slender body at any point *s*_{A} along its centreline. The first two terms in equation (2.21) represent the temperature owing to a local straightline source with length 2*l*. The third term in equation (2.21) accounts for the line distribution along each of the other cylindrical elements making up the curved body. The second term in equation (2.13) has disappeared, because the local heat source no longer varies with time. Equation (2.21) provides a description of the interaction through steady-state heat conduction among a collection of cylindrical elements that approximate a curved slender body that is equivalent to the treatment in equation (2.18) of Mackaplow *et al.* [36] of steady-state heat conduction in a composite with many discrete straight highly conductive fibre inclusions.

## 3. Numerical implementation

In §3, we develop one possible numerical framework for the SBT model for transient heat conduction derived in the previous section. Figure 2 shows with a simplified diagram the six steps taken in developing this framework. The slender body heat transfer problem is discretized in the space (step 1), and time domain (step 2), which converts the integrals in the SBT solution (equation (2.13) or (2.16)) into summations. The most cumbersome calculation is the double summation over all previous heat pulses emitted by all neighbouring slender body elements. Therefore, simplifications are applied in the space (step 3) and time domain (step 4) when calculating these heat pulse contributions. The computational time can further be decreased by combining point sources in the space (step 5) as well as in the time domain (step 6), a method also known as the fast multipole method (FMM). The following sections explain each step in more detail and develop the SBT numerical equations, which are implemented in a computer algorithm.

### (a) Space and time discretization of slender-body theory model

A computer implementation of the SBT model requires to numerically solve equation (2.13) or (2.16) for each location *s* and time *t*. The integrals in these equations can be converted into summations by discretizing the slender body into a set of straight pipe elements *i* and discretizing the time into time steps *m* (steps 1 and 2 in figure 2). This converts the continuous functions *Q*(*s*,*t*) and *T*(*s*,*t*) into the discrete functions *Q*_{i,m} and *T*_{i,m}, respectively. *Q*_{i,m} represents the constant heat exchange (in units of W m^{−1}) for the entire pipe element *i* during time step *m*; *T*_{i,m} represents the temperature at the mid-point of pipe element *i* at the end of time step *m*.

The local geometry of the slender body sets a maximum length for the pipe elements by requiring that the set of elements still accurately represents the shape of the body. In numerical tests with the SBT algorithm, we used sample calculations with a heated toroidal ring to show that a maximum element length to radius of curvature ratio of about 0.7 results in at least 98% accuracy. On the other hand, there is no minimum length and the user can discretize the body as fine as deemed necessary to achieve the desired accuracy with element lengths on the order of the pipe radius or smaller possible. However, to simplify the numerical implementation, we choose here to set the element lengths, labelled as *L*_{i} for element *i*, equal to the length of the local constant line source (2*l*) of the outer solution (see §2b). Each element length *L*_{i} is now larger than the local pipe radius but not necessarily identical to the other element lengths, because *l* can vary along the slender body axis. For spatially smoothly varying *Q* along the centreline of the slender body, this simplification is not believed to jeopardize the accuracy. Further, it automatically sets the third term in equation (2.13) to zero, a simplification already applied when deriving equation (2.16).

There are no model constraints on the time step length. The user can choose non-uniform time steps as fine or coarse as required in obtaining the desired trade-off between accuracy and computational time. To simplify the numerical implementation, we choose to have the starting time *t*_{0} of the constant line source *Q*(*s*_{A},*t*) so that *t*−*t*_{0} corresponds to each time step length.

For the space and time discretization schemes selected above, equation (2.16) (or (2.13)) is now converted into a linear equation in *Q* with *T*_{i,m} the temperature at the selected element *i* after current time step *m* and *Q*_{j,n} the heat source emitted by element *j* during time step *n*:
*j*=*i*)). Further, term 4 in equation (2.13) or term 3 in equation (2.16) is split into two parts in equation (3.1) (third and fourth term) by distinguishing between current (*n*=*m*) and previous (*n*<*m*) heat pulses emitted by neighbouring elements (*j*≠*i*). The terms *f*_{i,j,m,n} now represent standard, easy to calculate models, which are the topic of §3b.

Other approaches exist to discretize or numerically solve equation (2.13) or (2.16). One could decouple 2*l* from the pipe element lengths to allow for very fine grids in order to capture abrupt spatial changes in *Q*, occurring over length scales on the order of the pipe radius or smaller. Another approach is to represent *Q*(*s*,*t*) as a summation of polynomials and numerically solve for the polynomial coefficients, as highlighted by Batchelor [20]. However, these methods are not further investigated here.

### (b) Equations for *f*_{i,j,m,n}

*f*_{i,j,m,n}, introduced in equation (3.1), can be interpreted as the temperature response at the mid-point of pipe element *i* at time step *m* caused by a unit heat pulse of 1 W m^{−1} emitted by pipe element *j* during time step *n*. From equations (2.13) and (2.16), it is derived that *f*_{i,j,m,n} represents simple models such as cylindrical, line or point source models, which can often be calculated analytically. To save on computational time, a method is presented in appendix A that uses decision trees and non-dimensional numbers, based on length- and time scales of the heat transfer problem, to select the most simple model that still accurately calculates each term *f*_{i,j,m,n}. For example, heat sources occurring far away in space and/or time (see steps 3 and 4 in figure 2) could be represented by point sources or even neglected rather than calculated with the full cylindrical or line source models. Appendix A provides all equations to calculate *f* for each term in equation (3.1) and explains in more detail the accuracy obtained by applying these simplifications.

### (c) Fast multipole method, other implementation strategies and computational time

The computational time of the SBT algorithm can further be decreased by combining point sources in the space and time domain (steps 5 and 6 in figure 2). This method is called the FMM and was pioneered by Greengard & Rokhlin in 1987 [44]. Although advanced decision criteria could be developed to determine which point sources to combine at which times, here a more rudimentary method is adopted. Previous pulses once represented by the most simple point source model (i.e. the PS_{3} model, introduced in appendix A) can be combined if _{L,PP,E} non-dimensional numbers defined in table 1 in appendix A, and appearing in the PS_{3} equation, is larger than a certain value, e.g. 100, and correspondingly the *f* become relatively small (but not negligible). A certain number of point sources, e.g. 5, could then be grouped together and represented by a new point source with a *Q*⋅(*t*_{n}−*t*_{n−1}) equal to the sum of the *Q*⋅(*t*_{n}−*t*_{n−1}) of the original point sources (conservation of energy) and with the new *t*_{n} and *t*_{n−1} (end and start time of a heat pulse) set in this way that the new *f* exactly matches the sum of the old *f*, e.g. at the time of combining the pulses. This simplification can be done at multiple levels, i.e. five new point sources each replacing five old point sources can later be combined themselves. Also, in the space domain, different pipe elements located at different positions can have a similar spacing, e.g. within 5%, with respect to a selected element. Once represented by a point source, the heat pulses from these different elements can simply be added together and only one *f* needs to be calculated.

Other implementation strategies to decrease the computational time are to calculate certain terms that often reoccur throughout the simulation once in advance. Examples are (*L*_{i}/4*πk*), several non-dimensional numbers, special functions such as erf(*x*) and *E*_{1}(*x*) and the integrals in the cylindrical source model. Also, because each equation (3.1) for every pipe element is coupled with each other, the solution algorithm requires solving a linear set of equations at each time step, for which either direct or iterative methods can be used. Because the matrix system in some cases is sparse and/or symmetrical, fast solution methods could be used, which are built-in into software packages such as MATLAB [45].

An SBT computer program based on the numerical framework presented requires implementing equation (3.1), the non-dimensional numbers, heat source model equations and decision trees from appendix A to calculate *f*_{i,j,m,n}, and depending on the problem possibly equations from §2d (model extensions), e.g. to include heat transfer phenomena inside the slender body or for specific boundary conditions. A simulation requires solving equation (3.1) at each time step simultaneously for each pipe element. For SBT, the total number of elements *N*_{E} is *N*_{L} the discretization of the slender body. Without any simplifications (steps 3–6 of figure 2 neglected), and with using a direct solver based on LU decomposition for solving the matrix system, the computational time is *N*_{t} the number of time steps. In contrast, for FEM, with *L* a characteristic length for the medium domain, the number of elements *N*_{E} is *N*_{E} for SBT is typically several orders of magnitude smaller than *N*_{E} for FEM which results in lower computational times in many cases. However, the SBT model could no longer be advantageous for more spatially complex but compact problems, especially in combination with a large number of time steps. With the simplifications and strategies provided, the computational time will drop but hard to determine by how much because it depends on several factors such as decision tree limits, matrix solver algorithm, and spatial and time lengths and discretizations. In terms of computer memory, the SBT model can also be superior, because an FEM simulation typically requires large amounts of data storage for meshing, matrix assembly and solving.

## 4. Case studies

Here, two case studies are presented to validate and illustrate the applicability, speed and versatility of the SBT numerical model developed in §3. The first case study is a theoretical example of two parallel cylinders in close proximity and at constant temperature; the second case study is an applied example of one loop of a slinky-coil heat exchanger. In each case, the SBT numerical model, implemented in MATLAB R2012a using all of the simplifications and strategies discussed, is compared with the FEM software package COMSOL Multiphysics v. 4.3 [46] in terms of computational time. In order to make a fair comparison, the discretization schemes are chosen such that comparable accuracy is obtained. All simulations are carried out on a standard Intel Core i3-2120 (3.3 GHz) desktop computer with 6 GB RAM and Windows 7 as operating system. Both in MATLAB and COMSOL, a direct solver is used to solve the linear system of equations.

### (a) Two parallel cylinders at constant temperature

Two parallel cylinders with length *L*_{t}=10 m, radius *R*=0.05 m and spacing *S*=1 m are placed in an infinite conductive medium with thermal conductivity *k*=1 Wm^{−1} K^{−1} and thermal diffusivity *α*=10^{−6} m^{2} s^{−1} (figure 3). The medium has an initial temperature of 0°C and both cylinders are kept at a constant temperature of 10°C. We aim to solve for the transient heat exchange *Q* (W m^{−1}) at both short and long time scales, occurring at the mid-point (*δ*_{m}=5 m) and close to the endpoint (*δ*_{e}=0.25 m) of the cylinders, averaged over the circumference. Owing to the symmetry of the problem, both cylinders will have the same heat flux.

In the SBT model, each cylinder is discretized into 20 identical pipe elements for a total number of 40 grid points, which is identical to modelling one cylinder and applying the method of images to force a zero flux boundary condition at the mid-plane between both cylinders. The time is discretized into 150 time steps from 0 to 10^{9} s. The COMSOL model is a three-dimensional model of a quarter of one cylinder (axial cross section of the half-length), with zero heat flux condition at the appropriate boundaries, in order to make full use of the symmetry embedded in this problem, and constant temperature of 0°C at the other (far-field) boundaries. This quarter of a cylinder is placed in a cuboid with height of 35 m and width and depth of 30 m, at its correct location close to one of the corners. The number of grid points along the cylindrical body is 500, the total number of mesh elements is about 130 000 and the relative tolerance is set to 0.1. With these spatial and time discretization schemes, the SBT and COMSOL model show a comparable accuracy of about 99% on average over the entire time span, with reference to a highly accurate COMSOL model with very fine mesh (1 000 000 cells) and very low relative tolerance (10^{−6}). This reference COMSOL model has been validated at short time scales with the analytical solution of a single, infinitely long cylinder at constant temperature [1].

The result for *Q* at the mid-point and endpoint as a function of time is shown in figure 4 with the heat flux and time non-dimensionalized as *Q*/(*k*Δ*T*) and *T* the temperature difference between the cylindrical surface and the far-field (Δ*T*=10°C). The non-dimensional time can be interpreted as the ratio of the diffusion distance *S*. The result for one cylinder only is included to illustrate the tube–tube interactions which start to occur at a non-dimensional time of about 0.5 corresponding to a thermal diffusion length of 0.5 m. The difference in *Q* between the mid-point and endpoint becomes evident at a non-dimensional time of about 0.2, corresponding to a thermal diffusion length of 0.2 m. The computational time for the SBT model is 7 s and for the COMSOL model about 10 min. This case study illustrates the high computational speed of the SBT model, and its versatility in easily handling short and long time scales and multiple, disconnected, interacting slender bodies of finite length. Even higher speeds in this as well as in the next case study can likely be obtained when implementing the SBT model using a compiled language such as C instead of an interpreted language as with MATLAB. However, this is not further investigated.

### (b) Slinky-coil heat exchanger

In the second case study, the outlet temperature of water flowing inside one loop of a slinky-coil heat exchanger as a function of time is simulated for (figure 3). The heat exchanger has a total length *L*_{t} of 3.325 m, a loop radius *R*_{C} of 0.37 m, a pipe radius *R* of 1 cm, a pipe wall thickness *t*_{p} of 0.1 cm and a pipe material thermal conductivity *k*_{p} of 1 Wm^{−1} K^{−1}. At the cross-over point, the pipes are within 0.5 cm of each other. The density *ρ*_{f}, specific heat capacity *c*_{p,f}, inlet temperature *T*_{in}, and mean velocity *U* of the water are 1000 kg m^{−3}, 4200 J kg^{−1} K^{−1}, 20°C, and 0.1 m s^{−1}, respectively. The convective heat transfer coefficient between the water and the pipe is set at *h*=1000 Wm^{−2} K^{−1}. Further, the infinite conductive medium has a thermal diffusivity *α* of 10^{−6} m^{2} s^{−1}, thermal conductivity *k* of 3 Wm^{−1} K^{−1}, and initial temperature of 0°C.

For the SBT equations, the heat exchanger is discretized into 30 identical pipe elements and the time is discretized into 50 linear time steps from 1 to 50 s, followed by 100 logarithmically spaced time steps from 50 s up to 10^{7} s. In the COMSOL model, the heat exchanger is placed at the centre of a cube, representing the infinite medium, with edge length of 10 m. The total number of elements is about 80 000, and the relative tolerance is set to 0.02. Furthermore, in both the SBT and COMSOL simulations, the conductive medium model is coupled to a one-dimensional advective heat flow model representing the fluid flow inside the slinky-coil heat exchanger (see §3c) with a spatial discretization of 180 and 200 elements in the SBT and COMSOL model, respectively. The spatial and time discretization schemes for both models chosen provide a comparable accuracy of about 99% on average over the entire time span, in comparison with a highly accurate COMSOL model with very low relative tolerance of 10^{−6} and a fine mesh of about 350 000 cells.

The water outlet temperature simulated with the SBT and COMSOL model is shown in figure 4. The calculation time for the SBT model is 6 s, whereas for the COMSOL model about 19 min. This case study illustrates the SBT model running at high computational speeds and easily handling short and long time scales, curved bodies, and heat transfer inside a slender body. Moreover, the slender body approaching itself at the cross-over point, which is a violation of one of the constraints provided in §2a, does not appear to have a notable impact on the outlet temperature. The reason is that the two parts only lie within the inner region of one another for a short distance *L*_{t}/*R*≫1, the heat transfer from this region is small

Slinky-coil heat exchangers are typically modelled using FEM numerical software [47–50] or g-functions [51]. This case study shows SBT can be an alternative, computationally fast approach. Other types of heat exchangers for geothermal or ground-source heat pump applications, e.g. parallel horizontal pipes, and single-U (1U) and double-U (2U) borehole heat exchangers, can also be simulated using SBT. Potential oversimplifications often applied in other models for simulating these heat exchangers, such as constant longitudinal heat exchange *Q* along a single-U borehole heat exchanger, can then be avoided.

## 5. Conclusion

The SBT was initially developed and is most often used in the field of fluid dynamics to simplify the flow and pressure field solution of a slender body problem using matched asymptotic expansions of an inner and outer solution. In this work, the SBT is developed for slender, curved bodies with circular cross section undergoing transient heat transfer in a conductive medium. The theoretical derivation is based on matching the temperature field of an infinitely long cylindrical source as inner solution with the temperature field of an FLS as outer solution. By replacing the line source model with the cylindrical source model in the final result, a solution is obtained applicable at all time scales. One possible numerical model is then developed by discretizing the SBT theoretical result in the space and time domain. High computational speeds are obtained by applying simplifications such as modelling cylindrical and line sources as point sources, combining point sources with the FMM, or even neglecting point sources when allowed. Two case studies are included to demonstrate the SBT model operating at high computational speeds, capturing tube–tube interactions and finite lengths, and handling curved bodies, short and long time scales, and heat transfer processes inside a slender body.

The SBT for transient heat conduction can be an attractive solution technique, especially in the case of complex and non-compact geometries, when no analytical solutions are available and finite-element methods require long computational times and large amount of computer memory. Examples of applications are simulating the heat transmission in deviated oil, gas or geothermal wells, designing geothermal reservoirs for underground energy storage or heat pump applications, and modelling heating or cooling elements in industrial processes.

## Data accessibility

This work does not have experimental data.

## Authors' contributions

K.F.B drafted the manuscript, developed the numerical model and implemented the two case studies. D.L.K., K.F.B. and J.W.T. derived the theoretical basis and interpreted the results. All authors gave final permission for publication.

## Competing interests

We have no competing interests.

## Funding

We acknowledge the US Department of Energy (contract nos. DE-EE0002745 and DE-EE0002852), the Verizon Foundation, Verizon Wireless, the Atkinson Center for a Sustainable Future, the Cornell Energy Institute, the Geothermal Resources Council Scholarship, the Fulbright Program, the Commission for Educational Exchange between the United States of America, Belgium and Luxembourg, Vesuvius plc and the King Baudouin Foundation for partial financial support of this work.

## Appendix A. Non-dimensional numbers, heat source model equations and decision trees for calculating *f*_{i,j,m,n}

This appendix provides the procedure and equations for calculating *f*_{i,j,m,n}, the terms that emerged in equation (3.1) after discretizing the SBT theoretical result for developing the numerical framework (§3), which was used in the case studies (§4). As mentioned in §3b, *f*_{i,j,m,n} can be interpreted as the temperature rise at the mid-point of element *i* at the end of time step *m* owing to a unit heat pulse by element *j* during time step *n*. Without simplifications, the numerical framework would require for each pipe element at each time step to evaluate the full equation (2.13) or (2.16), i.e. calculating the full finite cylindrical or line source model and the contributions from all previous and neighbouring heat pulses. However, simplifications are possible to save on computational time with no significant loss in accuracy. These simplifications are neglecting the finite radius or length of a pipe element, or neglecting the element as a whole, which translates into different models for calculating *f*_{i,j,m,n}.

Six well-known models are identified to calculate *f*_{i,j,m,n}, ordered from more to less representing the actual geometry and from more to less computationally intense to evaluate: the finite cylindrical source (FCS) model, the infinite cylindrical source (ICS) model, the finite line source (FLS) model, the infinite line source (ILS) model, the point source (PS) model and the point neglected (PN) model. The PN model simply means setting the temperature contribution to zero by neglecting the element. The equations for these heat source models, as well as the methodology to select the most appropriate model rely on dimensionless numbers, listed in table 1. In these numbers, *R*_{j} and *L*_{j} are the radius and length of pipe element *j* emitting the heat pulse, *S*_{i,j} is the spacing between the centres of pipe element *i*, where the temperature is calculated at, and pipe element *j*, the source of the heat pulse, *TDR* refers to the ratio of time durations, *SLR* refers to spacing-to-length ratio, *CP* and *PP* refer to current pulse and previous pulse and *S* and *E* refer to start and end of a previous heat pulse, respectively. Furthermore, the time steps are in general labelled as *t*_{k} with *t*_{0} the time starting point set at 0 s. *t*_{m} is the current time and *t*_{n} is the end time of heat pulse *n*. Also, (*t*_{m}−*t*_{m−1}) is the time duration of the current heat pulse, and (*t*_{m}−*t*_{n−1}) and (*t*_{m}−*t*_{n}) are the time that has past since the beginning and end of heat pulse *n*, respectively.

The equations to calculate *f*_{i,j,m,n} are provided in table 2. They are grouped into four categories, corresponding to the four terms in equation (3.1), by distinguishing between selected (*j*=*i*) and neighbouring (*j*≠*i*) pipe elements, and current (*n*=*m*) and previous (*n*<*m*) heat pulses. In the equations, *J*_{1}(*x*), *Y* _{1}(*x*), erfc(*x*) and *E*_{1}(*x*) refer to first-order Bessel function of the first kind, first-order Bessel function of the second kind, complimentary error function and exponential integral, respectively [37]. The temperature contribution from a previous heat pulse *Q*_{i,n}, emitted during the time period *t*_{n−1} to *t*_{n}, is calculated as the sum of the contributions by a heat pulse with the same strength active from *t*_{n−1} to *t*_{m} (current time) and a heat pulse with opposite strength active from *t*_{n} to *t*_{m}.

Figure 5 provides decision trees that have been developed to recommend the least computationally intense model that still accurately calculates *f*. The limits in the decision trees are calculated to obtain an accuracy of at least 99% in calculating dominant *f* and at least 95% for non-dominant *f* in equation (3.1). The dominant terms are the *f* for the current and *recent* pulses for the selected pipe and the *nearby* pipes and are characterized by having a value of up to several orders of magnitude larger than the non-dominant *f*. Different decision tree limits could be developed for low-, medium- and high-accuracy simulations with corresponding high, medium and low computational speed. The decision tree limits have been derived in MATLAB by finding the values of the non-dimensional numbers for which the solution of a certain simplified model deviates no more than e.g. 99% with respect to the solution of the most accurate model available, e.g. the FCS model in case of the first two decision trees.

- Received July 23, 2015.
- Accepted November 13, 2015.

- © 2015 The Author(s)