A spectral model for a moving cylindrical heat source in a conductive-convective domain

This paper introduces a spectral model for a moving cylindrical heat source in an inﬁnite conductive- convective domain. This physical process occurs in many engineering and technological applications including heat conduction-convection in ground source heat pump systems, where the borehole heat ex- changers likely go through layers with groundwater ﬂow. The governing heat equation is solved for Dirichlet and Neumann boundary conditions using the fast Fourier transform for the time domain, and the Fourier series for the spatial domain. A closed form solution based on the modiﬁed Bessel functions is obtained for the Dirichlet boundary condition and an integral form for the Neumann boundary condi- tion. Limiting cases of the moving cylindrical heat source to represent a moving line heat source are also derived. Compared to solutions based on the Green’s function and the Laplace transform, the spectral model has a simpler form, applicable to complicated time-variant input signals, valid for a wide range of physical parameters and easy to implement in computer codes. The model is veriﬁed against the existing inﬁnite line heat source model and a ﬁnite element model.


Introduction
Heat conduction-convection in an infinite domain subjected to a heat source may be regarded either a case with a moving source of heat through a conductive medium, or a case of a convective medium passing a heat source. The moving heat source or medium arises in many technological and engineering applications, of which are the machining process [11] and the ground source heat pump (GSHP) system [2] . This paper focuses on the GSHP system, a rapidly growing renewable energy technology, primarily constituting a vertical borehole heat exchanger that collects and rejects heat from and to the ground to be used for heating and cooling of buildings. The borehole can go as deep as 200 m in the ground and likely passing soil layers with groundwater flow, giving rise to conduction-convection heat flow.
Modeling conduction-convection heat flow in a saturated porous domain subjected to cylindrical Dirichlet or Neumann boundary conditions can readily be made using numerical methods such as the finite element method, the finite volume method or a hybrid between them. Among others, Al-Khoury et al. [3] and Nam * Corresponding author.
E-mail address: r.i.n.alkhoury@tudelft.nl (R. Al-Khoury). et al. [19] utilized the finite element method to simulate heat flow in GSHP systems constituting layers with groundwater flow. Yavuzturk et al. [25] and Nabi and Al-Khoury [17,18] utilized the finite volume method for this purpose. Recently, Cimmino and Baliga [10] developed a computational model utilizing the control-volume finite element method (CVFEM), which is a hybrid between the finite element method and the finite volume method, for modeling heat flow in ground source heat pump systems under the effect of groundwater flow. Despite the versatility of the numerical methods they are computationally demanding, and the finite element method, in particular, exhibits shortcomings in simulating highly convective problems. These make the analytical methods more appealing.
However, modeling this system analytically is intricate due to the involved mathematical procedure to solve the governing partial differential equations exactly (see [24] for analytical solutions of heat transfer problems). Carslaw and Jaeger [9] were among the firsts to introduce analytical solutions to heat conductionconvection problems using the moving heat source/domain concept. They provided analytical solutions to the moving infinite line heat source using the Green's function method and to the moving infinite cylindrical heat source using the Laplace transform. The Green's function solution of Carslaw and Jaeger has been employed intensively to formulate analytical models for GSHP systems involving groundwater flow. Sutton et al. [22] adopted the Green's function solution for a moving infinite line source to simulate 2D heat flow in GSHP systems with groundwater flow. Similarly, Diao et al. [12] adopted Carslaw and Jaeger's solution to simulate the thermal interaction in a domain consisting of multiple infinite line heat sources. Capozza et al. [8] implemented the moving infinite line source model and provided three solution forms: closed analytic, asymptotic and tabulated.
Molina-Geraldo et al. [16] extended the moving infinite line source to a moving finite line source in a 3D domain using the Green's function method, and provided an elegant comparison between the two line source models. Similarly, Rivera et al. [20] introduced an analytical solution to the moving finite line source, taking into consideration the surface temperature and a spatially variable land use. More recent, Erol and Francois [14] elaborated on Molina-Geraldo et al. [16] model to simulate a moving finite line heat source in a multilayer domain including groundwater flow.
In a further development in modeling the moving heat sources, Zhang et al. [26] extended the 3D finite line source model to a moving finite cylindrical heat source using the Green's function method. Likewise, Zhou et al. [27] formulated a moving cylindrical heat source model for GSHP systems.
Apparently, the Green's function seems to be the technique that is largely being adopted for modeling convective-conductive heat flow in GSHP systems. Its formulation is relatively simple compared to that derived from the Laplace transform, making it easier to handle in computer codes. In Section 6 we discuss these two solutions and give examples of their formulations. Basically, the Green's function solution of the line heat source is in reality a solution for an instantaneous point source with a unity thermal load, where one or more points aggregated together to simulate infinite or finite line heat sources. The solution for the moving infinite line source is represented by a single integral over time, and that for the moving finite line source represented by a double integral: one over time and another over the length of the heat source. For an infinite cylindrical heat source, however, the solution requires a triple integral: one over the azimuth of the point, another over the azimuth of the cylinder and one over time. The integrals in many cases cannot be evaluated in closed forms and might impose difficulties to solve numerically, especially for relatively far distances and highly convective velocities. Yet, there are good attempts to facilitate these integrals, see for instance Zubair and Chaudhry [28] , Molina-Geraldo et al. [16] , Zhang et al. [26] and Zhou et al. [27] .
In this paper we depart from the Green's function and the Laplace transform and introduce a spectral model for heat flow in a convective infinite domain passing a cylindrical heat source. The heat flow is analyzed in a moving 2D domain in the xy −plane, shown schematically in Fig. 1 . The spectral analysis is utilized to solve the governing heat equation for prescribed Dirichlet and Neumann boundary conditions. This method relies basically on the Fourier series to solve the governing equation in the spatial domain and the fast Fourier transform (FFT) to solve it in the temporal domain. Compared to models based on the Green's function and the Laplace transform, the spectral model has a simpler formulation, applicable to complicated time-variant input signals, computationally efficient and easy to implement in computer codes. Its Fourier series summation (rather than integration) makes the solution less restrictive on the model physical and thermal parameters. Additionally, unlike the inverse Laplace transform, inverse calculation of the spectral model is rather straightforward due to the use of the inverse fast Fourier transform (IFFT) algorithm.

Governing equations
The equation governing heat conduction-convection in an infinite, homogeneous, isotropic domain, moving with velocity −U along the x -axis, can be expressed as is the temperature (K) in an xy −plane, and α = λ/ (ρc) is the material thermal diffusivity (m 2 /s) with λ the thermal conductivity (W/(m.K)), ρ the mass density (kg/m 3 ) and c the specific heat capacity (J/(kg.K)). The domain is initially at zero temperature, and for times t > 0 it is subjected to Dirichlet or Neumann boundary conditions at x 2 + y 2 = a 2 , with a the radius of the heat source (see Fig. 1 ). The Dirichlet boundary condition entails T (a, θ , t ) = T in (t) (2) where T in is any time-dependent input temperature signal. The Neumann boundary condition implies where Q in is any time-dependent input heat flux signal, and r, θ are the polar coordinates of the cylindrical heat source. Since the cylindrical cross section is periodic (continuous) in the azimuth ( θ − direction), the following physical constraints are maintained:

Solving the governing equations
The solution to Eq. (1) can be expressed in a more convenient way [9] , as (5) where u ( x, y, t ) is some function in space and time, and κ is a constant, needs to be determined. Substituting Eq. (5) into Eq. (1) gives The spatial derivative in the third term of this equation can be and Transforming this equation to the polar coordinate system in xy − plane, yields for which x = r cos θ , y = r sin θ and r = x 2 + y 2 have been employed. Considering Eqs. (2) and (7) , the relevant Dirichlet boundary Similarly, considering Eqs. (3) and (7) , the relevant Neumann boundary condition can be expressed as In the same way, the physical constraints in Eq. (4) can be denoted as

Spectral analysis
Applying the Fourier transform to Eq. (9) yields where ˆ u ≡ˆ u (r, θ , ω) , with ω the angular frequency, i = √ −1 , and the transform is represented by Eq. (13) can be solved using the separation of variables, which can be expressed as ˆ u (r, θ , ω) = R (r, ω) (θ ) (15) Substituting this equation into Eq. (13) , dividing by R ( r, ω) ( θ ) and equating both terms to a constant, n 2 , gives The solution to this ordinary differential equation can be expressed as n (θ ) = a cos nθ + b sin nθ (18) where a and b are the integration constants, which need to be determined form the boundary conditions.
The periodicity in Eq. (12) gives rise to an eigenvalue problem with n being integer, leading to (θ ) = n a n cos nθ ; It is noteworthy indicating that the solution to Eq. (17) can also be expressed as (θ ) = n a n e i n θ ; Comparing the two solutions it can readily be seen that the summation in Eq. (20) requires double that in Eq. (19) , and hence, in what follows, the cosine formulation will be pursued.

R − Dependency
The R term in Eq. (16) gives The solution to this ordinary differential equation can be expressed in terms of the modified Bessel function of first and second kinds [1] , but as the first kind is unbounded at far distances, the solution reduces to R (r, ω) = K n (qr) (22) where K n is the modified Bessel function of the second kind of order n .

Solving for Dirichlet boundary condition
Applying the Fourier transform to the boundary condition in Eq.
Having the solution for u ( r, θ , t ), using Eq. (7) in which J n is the Bessel function of the first kind of order n .

Solving for Neumann boundary condition
Taking the derivative of Eq. (23) with respect to r and substituting into the transformed form of Eq. (11) gives n a n cos nθ This equation can be expressed as n a n cos nθ = f (θ ) ; The coefficients of this Fourier cosine series are As f ( θ ) in Eq. (34) cannot be arranged in a form similar to that of Eq. (27) , it seems that there is no closed form solution to Eq. (35) and the integrals have to be solved numerically. However, the integrals can readily be solved using any appropriate numerical solver. Here, we use the the default integral algorithm of MATLAB [15] , which is based on a vectorised adaptive quadrature given by Shampine [21] . The solution can then be expressed as ˆ u (r, θ , ω) = n a n cos nθ K n (qr) ; n = 0 , 1 , 2 , · · · Using Eq. (7) , the temperature distribution in the domain can then be determined as Should we have pursued the exponential formulation, Eq. (20) , the solution becomes ˆ u (r, θ , ω) = n a n e inθ K n (qr) ; n = 0 , ∓1 , ∓2 , · · · (39) with a n = 1

Limiting Cases
As the radius of the cylinder approaches zero, the solution of the moving infinite cylindrical heat source should lead to the solution of the moving infinite line source (moving point source).
The Neumann boundary condition for the line heat source in the frequency domain is Similar to deriving Eq. (34) , Eq. (41) leads to a Fourier cosine series of the form n a n cos nθ = f (θ ) | r→ 0 ; f (θ ) = αˆ Q in (ω) e Ur cos θ/ 2 α πλ( K n (qr) Ur cos θ +2 K n +1 (qr ) αqr −2 K n (qr ) α n ) (42) To determine the a 0 coefficient, we solve Eq. (42) for n = 0 , yielding f (θ ) = αˆ Q in (ω) e Ur cos θ/ 2 α πλ( K 0 (qr) Ur cos θ + 2 K 1 (qr) αqr ) As r → 0, this equation leads to where these limiting forms are utilized [23] : Following this, the a 0 coefficient reads Taking the limit of Eq. (42) as r → 0, leads to a n = 0 where this limiting case has resulted from: Having the Fourier coefficients Eqs. (46) and (47) , the solution to the moving infinite line source can then be expressed as In the same way, the limiting case for the Dirichlet boundary condition yields where these limiting forms are utilized [1] : z is a small argument (51)

Model verification and application
The model is verified against analytical solutions and numerical solutions, viz.: 1 Comparing the spectral moving infinite cylindrical heat source model to its limiting cases and a commonly used Green's function infinite line heat source model by making the radius of the cylinder very small. 2 Comparing the spectral moving infinite cylindrical heat source model to a finite element model.
A two-dimensional, fully saturated porous domain, constituting a solid phase and a water phase, subjected to a cylindrical heat source is considered for this purpose. The geometry resembles a soil layer with groundwater flow passing a borehole heat exchanger. The domain is in local thermal equilibrium and the heat equation is averaged over the constituents, such that and where the subscripts w , s , eff denote the water phase, the solid phase and the effective parameter, respectively, and ϕ is the porosity. The thermo-physical parameters of the porous domain are given in Table 1 .
The initial and boundary conditions are:

Verification against analytical solutions
Computational results obtained from the spectral moving infinite cylindrical source (SMICS) model and its limiting case: the spectral moving infinite line source (SMILS) model, are compared to each other and to those obtained from the Green's function moving infinite line source (GMILS) model, given by Sutton et al. [22] .
To mimic the line source, the radius of the cylinder in the SMICS model is made too small, namely 0.0 0 01 m. Both Dirichlet and Neumann boundary conditions are tested. In literature, the GMILS model is given for a Neumann boundary condition only. To compare it with the spectral models for the Dirichlet boundary condition case, we first run the GMILS model and output its temperature at 0.0 0 01 m from the line source, then we prescribe this temperature on the cylindrical heat source in the SMICS model. Fig. 2 shows the temperature contours for the Dirichlet boundary condition case, for three water flow velocities: 1e-5, 1e-4 and 5e-4 m/s, and two physical times: 10 and 100 hr. The figure shows nearly perfect matching between the three models except for the case with relatively fast fluid velocity and longer time ( Fig. 2 c,  100 h), where the GMILS model exhibited some spurious oscillations. Fig. 3 shows the temperature contours for the Neumann boundary condition case for the above mentioned water flow velocities and physical times. The three models are in good agreement, except for the relatively fast fluid velocity and longer time, ( Fig. 3 c,  100 h), where the GMILS model exhibited some spurious oscillations, similar to that for the Dirichlet case.
The three models were implemented in MATLAB, and while conducting the calculations the spectral models have exhibited a robust capability in handling a wider range of fluid velocities and physical times. The runs were conducted for a step function, and thus there was no significant CPU time difference between the three models. However, if a random input signal has been employed, there would have been a significant difference between the spectral models and the Green's function model. Using the FFT, the spectral models are capable of dealing with any input signal in the same efficiency as that for a step function; whereas using the Green's function would require a temporal superposition procedure that entails discretizing the input signal into multiple step functions, for each the involved integrals need to be solved separately.

Verification against numerical solutions
To examine the accuracy of the spectral model in simulating a typical cylindrical heat source as existing in GSHP systems, computational results obtained from the SMICS model were compared to those obtained from the finite element package: Multiphysics, Multidomain, Multiphase finite element analysis (MMM-FEM); a comprehensive 3D thermo-hydrodynamic-mechanical finite element software for engineering applications in geothermal energy, CO 2 geo-sequestration, soil freezing and thawing and poromechanics in general; developed at Delft University of Technology [4][5][6] .
The finite element domain is a half cylinder, 20 m in radius, encompassing a cylindrical heat source with radius 0.075 m, a typical radius of a borehole heat exchanger. The mesh consists of 9900, 8-node hexahedron elements with a minimum size of 0.0026 m along the heat source circumference. The top view of the mesh with a zoom around the heat source is shown in Fig. 4 . Fluid velocities: U = 1 e − 5 m / s , and U = 1 e − 4 m / s were examined. Fig. 5 shows the temperature contours after 48 h for the Dirichlet case. The figure clearly shows a very good matching between the SMICS model and the finite element model.
The good matching in the Dirichlet boundary condition case is due to the exact description of the boundary condition at the elements nodes defining the heat source. This cannot be guaranteed for the Neumann boundary condition case as the heat flux in the finite element method is by definition interelement-discontinuous. The heat flux is discretized at the elements integration points and then mapped to the nodes; leading eventually to an error in describing the boundary condition. To circumvent this problem, we first conducted spectral analysis for a Neumann boundary condition and output the temperature profile at the source boundary ( a = 0 . 075 m ) versus time. This temperature profile was then prescribed at the source boundary in the finite element model. Fig. 6 shows the temperature contours after 48 hr. It shows that there is a good matching between the SMICS model for Neumann boundary condition and the equivalent Dirichlet case computed by the finite element model.

Comparing spectral model to Laplace transform and Green's function solutions
Existing solution techniques for a moving cylindrical heat source in a convective-conductive domain, applicable mainly to ground source heat pump systems, can be put into two categories: Laplace transform and Green's function. In this section we present the Laplace transform solution given by Carslaw and Jaeger [9] and the Green's function solution given by Zhang et al. [26] to highlight the difference in complexity of the mathematical formulations of the three models.

Laplace transform solution
Carslaw and Jaeger provided, in their seminal monograph on Heat Conduction in Solids [9] , an analytical solution to an infinite solid initially at constant temperature T 0 , moving along the   x − axis with velocity −U and crossing a cylindrical heat source x 2 + y 2 = a 2 maintained at zero temperature. Using the Laplace transform, they solved Eq. (1) as e −U 2 t / 4 α ε n cos nθ I n Ua in which Y n is the Bessel function of the second kind of order n . All other parameters are defined above. Obviously, the presence of the semi-infinite integral of this transcendental function with oscillatory Bessel functions makes the problem complicated and the calculation, if possible, must be conducted using numerical algorithms. Carslaw and Jaeger [9] provided no numerical results to this case, neither simplified the function for short or long times, as they usually did for many other cases. They also didn't provide a solution for the Neumann boundary condition case.

Green's function solution
Zhang et al. [26] provided an analytical solution to the infinite and finite cylindrical heat sources in convective-conductive domains. Their solution is based on the Green's function formulation of an instantaneous point source given by Carslaw and Jaeger [9] . Due to the unsymmetrical nature of heat convection arising from groundwater flow, heat emitted from a point is described by integrating over the azimuth of the point. To consider the combined effects of all points constituting the cylinder, the solution of a point is integrated over the polar angles of the cylinder. As the heat source is moving with time with a certain velocity, the solution is completed by integration over time. Fig. 7 shows schematically Zhang et al. conceptual model.
For the infinite cylindrical source case, this model is described by a triple integral, such that in which φ is the polar angle of the cylindrical heat source, and ϕ is the polar angle of the temperature distribution around a point source.
Calculating this equation requires: 1 A proper numerical integration scheme. Though, such integrals might be restrictive, especially for far distances and longer times. 2 Q in is a step function, but if it is a function of time, as it is usually the case, its time signal must be divided into several step functions, for each, this equation must be recalculated. 3 For detailed analysis to produce contour plots of temperature, this equation can take rather long CPU time. Zhang et al. [26] provided no numerical results for a convective case showing isothermal contours for an x -y plane; rather providing contour lines in the z −direction for two angles only: φ = 0 and φ = π .
Comparing these two solution techniques to those given in Eq. (29) , for the Dirichlet case, and Eqs. (34) -(36) for the Neumann,  it can readily be deduced that the spectral solution has a rather simpler formulation which makes it computationally efficient and robust.

Conclusions
Heat conduction-convection in porous solids occurs in many engineering applications and technologies, including the ground source heat pump systems; a technology designed to harvest thermal energy from shallow depths of the earth to be utilized for heating and cooling of buildings. In the presence of groundwater flow, the amount of harvested energy depends significantly on the amount of heat that is transported by the groundwater, and thus it is important to be considered in the design and analysis of the system. This paper addresses this issue and introduces an analytical solution to conductive-convective heat flow in an infinite domain subjected to a cylindrical heat source. The spectral analysis is utilized to solve the governing heat equation for prescribed Dirichlet and Neumann boundary conditions. The solution for the Dirichlet boundary condition has led to a closed form function, and that for the Neumann boundary condition has led to an integral which is relatively easy to solve numerically. Compared to solutions based on the Green's function or the Laplace transform, the spectral solution has a simpler formulation, as evident in Section 6 . This makes the spectral model computationally efficient and robust. Additionally, the spectral analysis is inherently applicable to any random input signal and its series summation (rather than integration) makes the solution less restrictive on the model physical and thermal parameters.
The spectral model presented in this paper can be incorporated within the spectral element method [13] and the superposition principle to study detailed heat flow in ground source heat pump systems constituting multiple borehole heat exchangers embedded in multilayer systems with layers encompassing groundwater flow. In particular, the given spectral model can readily be incorporated in the spectral element model of BniLam et al. [7] . This will be undertaken in a future work.

Authorship contributions
Please indicate the specific contributions made by each author (list the authors' initials followed by their surnames, e.g., Y.L. Cheung). The name of each author must appear at least once in each of the three categories below.

Declaration of Competing Interest
None.