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

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

Al-Khoury, Rafid; Bni Lam, Noori; Arzanfudi, Mehdi M.; Saeid, Sanaz

DOI

10.1016/j.ijheatmasstransfer.2020.120517

Publication date

2020

Document Version

Final published version

Published in

International Journal of Heat and Mass Transfer

Citation (APA)

Al-Khoury, R., Bni Lam, N., Arzanfudi, M. M., & Saeid, S. (2020). A spectral model for a moving cylindrical

heat source in a conductive-convective domain. International Journal of Heat and Mass Transfer, 163,

[120517]. https://doi.org/10.1016/j.ijheatmasstransfer.2020.120517

Important note

To cite this publication, please use the final published version (if applicable).

Please check the document version above.

Copyright

Other than for strictly personal use, it is not permitted to download, forward or distribute the text or part of it, without the consent of the author(s) and/or copyright holder(s), unless the work is under an open content license such as Creative Commons. Takedown policy

Please contact us and provide details if you believe this document breaches copyrights. We will remove access to the work immediately and investigate your claim.

International Journal of Heat and Mass Transfer 163 (2020) 120517

ContentslistsavailableatScienceDirect

International

Journal

of

Heat

and

Mass

Transfer

journalhomepage:www.elsevier.com/locate/hmt

A

spectral

model

for

a

moving

cylindrical

heat

source

in

a

conductive-convective

domain

Rafid

Al-Khoury

a,∗

_,

_Noori

_BniLam

b

_,

_Mehdi

_M.

_Arzanfudi

a,c

_,

_Sanaz

_Saeid

a

a Faculty of Civil Engineering and Geosciences, Computational Mechanics Chair, Delft University of Technology, P.O. Box 5048, 2600 GA Delft, the Netherlands b University of Antwerp – imec, IDLab - Faculty of Applied Engineering, Sint-Pietersvliet 7, 20 0 0 Antwerp, Belgium

c DIANA FEA BV., Thijsseweg 11, 2629JA Delft, the Netherlands

a

r

t

i

c

l

e

i

n

f

o

Article history: Received 13 July 2020 Revised 6 September 2020 Accepted 21 September 2020 Keywords:

Moving cylindrical heat source Conduction-convection heat flow Heat flow in groundwater Ground source heat pump

a

b

s

t

r

a

c

t

Thispaperintroducesaspectralmodel foramovingcylindricalheatsource inan infinite conductive-convectivedomain.Thisphysicalprocessoccursinmanyengineeringandtechnologicalapplications in-cludingheatconduction-convectioningroundsourceheatpumpsystems,wheretheboreholeheat ex-changers likely go through layers with groundwater flow. The governing heat equationis solved for DirichletandNeumannboundaryconditionsusingthefastFouriertransformforthetimedomain,and theFourierseriesforthespatialdomain.AclosedformsolutionbasedonthemodifiedBesselfunctions isobtainedfortheDirichletboundaryconditionandanintegralformfortheNeumannboundary condi-tion.Limitingcasesofthemovingcylindricalheatsourcetorepresentamovinglineheatsourcearealso derived.Comparedto solutionsbasedontheGreen’s functionand the Laplacetransform,thespectral modelhasasimplerform,applicabletocomplicatedtime-variantinputsignals,validforawiderangeof physicalparametersandeasytoimplementincomputercodes.Themodelisverifiedagainsttheexisting infinitelineheatsourcemodelandafiniteelementmodel.

© 2020TheAuthor(s).PublishedbyElsevierLtd. ThisisanopenaccessarticleundertheCCBYlicense(http://creativecommons.org/licenses/by/4.0/ )

1. Introduction

Heatconduction-convectioninan infinitedomainsubjectedto aheatsourcemayberegardedeitheracasewithamovingsource of heat through a conductive medium, or a case of a convec-tive medium passing a heat source. The moving heat source or mediumarisesinmanytechnologicalandengineeringapplications, of which are the machining process [11] andthe groundsource heat pump (GSHP) system [2]. This paper focuses on the GSHP system, arapidlygrowing renewableenergytechnology,primarily constituting a vertical borehole heat exchanger that collects and rejects heat from andto the groundto be used forheating and coolingofbuildings.Theboreholecangoasdeepas200minthe groundandlikelypassingsoillayerswithgroundwaterflow,giving risetoconduction-convectionheatflow.

Modeling conduction-convection heat flow in a saturated porous domain subjected to cylindrical Dirichlet or Neumann boundaryconditionscanreadilybemadeusingnumericalmethods such asthefiniteelementmethod,thefinitevolumemethodora hybridbetweenthem.Amongothers,Al-Khouryetal.[3]andNam

∗_{Corresponding author.}

E-mail address: r.i.n.alkhoury@tudelft.nl (R. Al-Khoury).

etal.[19]utilizedthefiniteelementmethodtosimulateheatflow inGSHPsystemsconstitutinglayerswithgroundwaterflow. Yavuz-turket al.[25] andNabi andAl-Khoury[17,18] utilizedthe finite volume method for this purpose. Recently, Cimmino and Baliga [10]developedacomputationalmodelutilizingthecontrol-volume finiteelementmethod(CVFEM),whichisahybridbetweenthe fi-niteelement methodandthefinitevolumemethod,formodeling heatflowingroundsourceheatpumpsystemsundertheeffectof groundwater flow. Despitethe versatility of thenumerical meth-ods they are computationally demanding, and the finite element method,in particular, exhibits shortcomings in simulatinghighly convectiveproblems.Thesemaketheanalyticalmethodsmore ap-pealing.

However, modeling thissystem analytically isintricate due to the involvedmathematical procedure tosolve the governing par-tialdifferentialequationsexactly(see [24]foranalytical solutions of heat transfer problems). Carslaw and Jaeger [9] were among the firsts to introduce analytical solutions to heat conduction-convection problems using the moving heat source/domain con-cept. They provided analytical solutions to the moving infi-nite line heat source using the Green’s function method and to the moving infinite cylindrical heat source using the Laplace transform.

https://doi.org/10.1016/j.ijheatmasstransfer.2020.120517

The Green’s function solution of Carslaw andJaeger hasbeen employedintensivelytoformulateanalyticalmodelsforGSHP sys-tems involving groundwater flow. Sutton et al.[22] adopted the Green’sfunction solutionforamovinginfinitelinesourceto sim-ulate2DheatflowinGSHPsystemswithgroundwaterflow. Simi-larly,Diaoetal.[12]adoptedCarslawandJaeger’ssolutionto sim-ulate the thermal interaction in a domain consisting ofmultiple infinitelineheatsources.Capozzaetal.[8]implementedthe mov-ing infinitelinesource modelandprovidedthree solutionforms: closedanalytic,asymptoticandtabulated.

Molina-Geraldo et al. [16] extended the moving infinite line source to a moving finite line source in a 3D domain using the Green’sfunctionmethod,andprovidedan elegantcomparison be-tween thetwo linesource models.Similarly,Rivera etal.[20] in-troduced an analytical solution to the moving finite line source, taking intoconsideration the surfacetemperature anda spatially variable land use. More recent,Erol andFrancois [14] elaborated on Molina-Geraldoet al. [16] model to simulate a moving finite line heat source in a multilayer domain including groundwater flow.

Inafurtherdevelopmentinmodelingthemovingheatsources, Zhang etal. [26] extended the 3D finite line source model to a moving finite cylindrical heat source using the Green’s function method. Likewise, Zhou etal.[27] formulated a moving cylindri-calheatsourcemodelforGSHPsystems.

Apparently,theGreen’sfunctionseemstobethetechniquethat is largely beingadopted formodeling convective-conductive heat flow in GSHP systems. Its formulation is relatively simple com-paredtothatderivedfromtheLaplacetransform,makingiteasier to handle in computer codes. In Section 6 we discuss thesetwo solutions and give examples of their formulations. Basically, the Green’sfunctionsolutionofthelineheatsourceisinreality a so-lutionforaninstantaneouspointsourcewithaunitythermalload, where one or more points aggregated together to simulate infi-niteorfinitelineheatsources.Thesolutionforthemovinginfinite linesourceisrepresentedby asingleintegralovertime,andthat forthemovingfinitelinesourcerepresentedbyadoubleintegral: oneovertimeandanotheroverthelengthoftheheatsource.For aninfinitecylindricalheatsource,however,thesolutionrequiresa tripleintegral:oneovertheazimuthofthepoint,anotheroverthe azimuthofthecylinderandone overtime.Theintegralsinmany cases cannotbe evaluated inclosedformsandmightimpose dif-ficultiestosolvenumerically,especiallyforrelativelyfardistances and highly convective velocities. Yet, there are good attempts to facilitate these integrals, see for instance Zubair and Chaudhry [28],Molina-Geraldoetal. [16], Zhangetal.[26] andZhou etal. [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 infinitedomain passing a cylindrical heatsource. Theheatflowisanalyzedinamoving2Ddomaininthe xy −plane, shown schematically in Fig. 1. The spectral analysis is utilized to solve thegoverning heatequation forprescribed Dirichletand Neumannboundaryconditions.Thismethodreliesbasicallyonthe Fourier series to solve the governing equation in the spatial do-main andthe fastFouriertransform(FFT)to solveit inthe tem-poral domain. Compared to models based on the Green’s func-tion and the Laplace transform, the spectral model has a sim-plerformulation,applicabletocomplicatedtime-variantinput sig-nals,computationallyefficientandeasytoimplementincomputer codes.ItsFourierseriessummation(ratherthanintegration)makes the solution less restrictive on the model physical and thermal parameters. Additionally,unlike theinverseLaplacetransform, in-verse calculation of the spectral model is rather straightforward due to the use of the inverse fast Fouriertransform (IFFT) algo-rithm.

Fig. 1. A schematic representation of the physical domain

2. Governingequations

The equation governing heat conduction-convection in an in-finite, homogeneous, isotropic domain, moving with velocity −U alongthe x -axis,canbeexpressedas

∂

2_T

∂

x2 +

∂

2_T

∂

y2 + U

α

∂

T

∂

x − 1

α

∂

T

∂

t = 0 (1)

inwhich T ≡ T(x, y, t )isthetemperature(K)inan xy −plane,and

α

=

λ

/

(

ρ

c

)

is the material thermal diffusivity (m2_{/s) with}

_λ

_the

thermalconductivity(W/(m.K)),

ρ

themassdensity(kg/m3_{)and c}

thespecificheatcapacity(J/(kg.K)).

Thedomainisinitiallyatzerotemperature,andfortimes t _>0 it is subjected to Dirichlet or Neumann boundary conditions at x 2_+y 2_{= a} 2_{,with a theradiusoftheheatsource(seeFig.1).The}

Dirichletboundaryconditionentails

T

(

a,

θ

,t

)

= Tin

(

t

)

(2)

where T _in is any time-dependent input temperature signal. The Neumannboundaryconditionimplies

−2

π

a

λ∂

T

(

r,

θ

,t

)

∂

r





r=a = Qin

(

t

)

(3)

where Q in is anytime-dependentinput heat flux signal, and r,

θ

arethepolarcoordinatesofthecylindricalheatsource.

Since the cylindrical cross section is periodic (continuous) in theazimuth(

θ

− direction), thefollowingphysicalconstraints are maintained: T

(

a,

θ

,t

)

= T

(

a,

θ

+ 2

π

,t

)

−2

π

a

λ

∂T(_∂r,_rθ,t)



_r =a= −2

π

a

λ

∂T(r,θ+2π,t) ∂r



r=a (4)

3. Solvingthegoverningequations

ThesolutiontoEq.(1)canbe expressedinamoreconvenient way[9],as

T

(

x,y,t

)

= u

(

x,y,t

)

eκx ₍₅₎

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

∂

2_u

∂

x2 +

∂

2_u

∂

y2 +

(

2

κ

+ U

α

)

∂

u

∂

x+

(

κ

2₊ U

α κ

)

u− 1

α

∂

u

∂

t = 0 (6)

Thespatialderivativeinthethirdtermofthisequation canbe eliminatedif

κ

=−U/2

α

,leadingto

R. Al-Khoury, N. BniLam, M.M. Arzanfudi et al. International Journal of Heat and Mass Transfer 163 (2020) 120517 and

∂

2_u

∂

x2 +

∂

2_u

∂

y2 − U2 4

α

2u− 1

α

∂

u

∂

t = 0 (8)

Transforming this equation to the polar coordinate systemin xy − plane, yields

∂

2_u

∂

r2 + 1 r

∂

u

∂

r + 1 r2

∂

2_u

∂

θ

2− U2 4

α

2u− 1

α

∂

∂

ut = 0 (9)

for which x =r cos

θ

, y =r sin

θ

and r =



x 2_+y 2_{have been}

em-ployed.

Considering Eqs. (2) and (7), the relevant Dirichlet boundary conditiontoEq.(9)at x =a cos

θ

is

u

(

a,

θ

,t

)

e−Uacosθ/2α= Tin

(

t

)

(10)

Similarly, considering Eqs. (3) and (7), the relevant Neumann boundaryconditioncanbeexpressedas

−2

π

a

λ



∂

u

(

r,

θ

,t

)

∂

r − U cos

θ

2

α

u

(

r,

θ

,t

)





r=a e−Uacosθ/2α_{= Q} in

(

t

)

(11) Inthesameway,thephysicalconstraintsinEq.(4)canbe de-notedas u

(

a,

θ

,t

)

≡ u

(

a,

θ

+ 2

π

,t

)

∂u(r,θ,t) ∂r



r=a≡ ∂u(r,θ+2π,t) ∂r



r=a (12) 3.1. Spectral analysis

ApplyingtheFouriertransformtoEq.(9)yields

∂

2_uˆ

∂

r2 + 1 r

∂

uˆ

∂

r + 1 r2

∂

2_uˆ

∂

θ

2− q 2_{uˆ = 0 ; q}2₌ U2 4

α

2+ i

ω

α

(13)

where u ˆ≡ ˆu

(

r,

θ

,

ω

)

, with

ω

theangularfrequency, i =√−1,and thetransformisrepresentedby

u ⇒ ˆ u ∂u ∂t ⇒ i

ω

uˆ ∂p_u ∂xp ⇒ ∂ p_uˆ ∂xp (14)

Eq.(13)canbesolved usingtheseparationofvariables,which canbeexpressedas

ˆ

u

(

r,

θ

,

ω

)

= R

(

r,

ω

)

(

θ

)

(15) SubstitutingthisequationintoEq.(13),dividingby R (r,

ω

)

(

θ

) andequatingbothtermstoaconstant, n 2_,gives

r2 R d2_R dr2 + r R dR dr − r 2_q2_{= −}1

d2

d

θ

2 = n 2 ₍₁₆₎ 3.1.1.



− Dependency The

terminEq.(16)is

d2

d

θ

2 + n

2

_{= 0 (17)}

The solution to this ordinary differential equation can be ex-pressedas

n

(

θ

)

= acos n

θ

+ bsin n

θ

(18)

where a and b aretheintegrationconstants,whichneedtobe de-terminedformtheboundaryconditions.

The periodicityinEq.(12) givesrisetoan eigenvalueproblem with n beinginteger,leadingto

(

θ

)

= 

nancos n

θ

; n = 0 , 1 , 2 ,· · · (19)

ItisnoteworthyindicatingthatthesolutiontoEq.(17)canalso beexpressedas

(

θ

)

=  _naneinθ ; n= 0 , ∓1 , ∓2 ,· · · (20)

Comparing the two solutions it can readily be seen that the summationinEq.(20)requiresdoublethatinEq.(19),andhence, inwhatfollows,thecosineformulationwillbepursued.

3.1.2. R − Dependency

The R terminEq.(16)gives

r2d2R dr2 + r

dR dr −

(

r

2_q2_{+ n}2

₎

_{R = 0 (21)}

The solution to this ordinary differential equation can be ex-pressedin termsofthemodifiedBesselfunction offirstand sec-ondkinds[1],butasthefirstkindisunboundedatfardistances, thesolutionreducesto

R

(

r,

ω

)

= Kn

(

qr

)

(22)

where K nisthemodifiedBesselfunctionofthesecondkindof

or-der n .

3.1.3. Solution of u ˆ

(

r,

θ

,

ω

)

SubstitutingEqs.(19)and(22)intoEq.(15)leadsto ˆ

u

(

r,

θ

,

ω

)

= 

nancos n

θ

Kn

(

qr

)

; n = 0 , 1 , 2 ,· · ·

(23)

3.2. Solving for Dirichlet boundary condition

Applying the Fourier transform to the boundary condition in Eq.(10)gives ˆ u

(

a,

θ

,

ω

)

= _Tˆ _in

₍

_ω

₎

_eUacosθ/2α ₍₂₄₎ At r ₌ a ,Eq.(23)becomes ˆ u

(

a,

θ

,

ω

)

=  nancos n

θ

Kn

(

qa

)

(25)

EquatingEq.(24)toEq.(25)yields ˆ

Tin

(

ω

)

eUacosθ/2α=



nancos n

θ

Kn

(

qa

)

(26)

ThisfunctionisaFouriercosineserieswithits coefficients ex-pressedas a0 = ˆ Tin(ω) K0(qa) 1 2π 2π 0 eUacosθ/2αd

θ

; n= 0 an = ˆ Tin(ω) Kn(qa) 1 π  2π 0 eUacosθ/2αcos n

θ

d

θ

; n ≥ 1 (27)

Theintegralsinthesecoefficientshaveaclosedformsolutionin termsofthemodifiedBessel functionofthe firstkindof0order (I 0)and n order(I n)([1],p.376),leadingto

a0= ˆ Tin(ω) K0(qa)I0

(

Ua 2α

)

; n = 0 an = 2_KTˆ_nin_(qa(ω)₎In

(

Ua_2α

)

; n ≥ 1 (28)

Thesolutioncanthusbeexpressedas ˆ u

(

r,

θ

,

ω

)

= n

ε

nT_Kˆin(ω) n(qa)In

(

Ua 2α

)

cos n

θ

Kn

(

qr

)

; n=0 , 1 , 2 ,· · ·

ε

0 = 1 for n= 0 ;

ε

n = 2 for n≥ 1 (29) ApplyingtheinverseFouriertransform,gives

u

(

r,

θ

,t

)

= 

m



nuˆ n

(

r,

θ

,

ω

m

)

e

iωmt ₍₃₀₎

wherethesummingindex m representsthefastFouriertransform (FFT)samplenumber.

Havingthesolutionfor u (r,

θ

, t ),usingEq.(7),thetemperature distributioninthedomaincanreadilybedeterminedas

T

(

r,

θ

,t

)

= u

(

r,

θ

,t

)

e−Urcosθ/2α (31) Shouldwehavepursuedtheexponentialform,Eq.(20),the so-lutionfortheDirichletboundaryconditionwouldleadto

ˆ u

(

r,

θ

,

ω

)

=  _n Tˆ in

(

ω

)

Kn

(

qa

)

in_J n

(

Ua 2 i

α

)

e inθ_K n

(

qr

)

; n = 0 ,∓1 ,∓2 ,· · · (32) inwhich J nistheBesselfunctionofthefirstkindoforder n .

3.3. Solving for Neumann boundary condition

Taking the derivative ofEq.(23)withrespect to r and substi-tutingintothetransformedformofEq.(11)gives



nancos_αnθe −Uacosθ/2α

(

K n

(

qa

)

Ua cos

θ

+2K n+1

(

qa

)

α

qa − 2K n

(

qa

)

α

n

)

=Qˆin (ω)

πλ

(33) Thisequationcanbeexpressedas



nancos n

θ

= f

(

θ

)

;

f

(

θ

)

= αQˆin(ω)eUacosθ/2α

πλ(Kn(qa)Uacosθ+2Kn+1(qa)αqa−2Kn(qa)αn)

(34)

ThecoefficientsofthisFouriercosineseriesare

a0 = ₂1_π  2π 0 f

(

θ

)

d

θ

; n = 0 an = _π1  2π 0 f

(

θ

)

cos n

θ

d

θ

; n ≥ 1 (35)

As f (

θ

) in Eq. (34) cannot be arranged in a form similar to that ofEq.(27),itseems that thereisnoclosedformsolutionto Eq.(35)andtheintegralshavetobesolvednumerically.However, theintegralscanreadilybesolvedusinganyappropriatenumerical solver. Here,we usethethedefaultintegralalgorithmofMATLAB [15],whichisbasedonavectorisedadaptivequadraturegivenby Shampine[21].

Thesolutioncanthenbeexpressedas ˆ

u

(

r,

θ

,

ω

)

=  _nancos n

θ

Kn

(

qr

)

; n = 0 , 1 , 2 ,· · · (36)

SimilartoEq.(30),applyingtheinverseFouriertransform,gives

u

(

r,

θ

,t

)

= 

m



nuˆ n

(

r,

θ

,

ω

m

)

e

iωmt ₍₃₇₎

Using Eq.(7), thetemperature distributionin the domaincan thenbedeterminedas

T

(

r,

θ

,t

)

= u

(

r,

θ

,t

)

e−Urcosθ/2α ₍₃₈₎

Shouldwe havepursuedtheexponentialformulation,Eq.(20), thesolutionbecomes

ˆ u

(

r,

θ

,

ω

)

=  _naneinθKn

(

qr

)

; n = 0 , ∓1 , ∓2 ,· · · (39) with an = ₂1_π  2π 0 f

(

θ

)

e−inθd

θ

f

(

θ

)

= 2aαQˆin(ω)eUacosθ/2α λ(Kn(qa)Uacosθ+2Kn+1(qa)αqa−2Kn(qa)αn) (40) 4. LimitingCases

As the radius ofthe cylinder approacheszero,the solution of themovinginfinitecylindricalheatsourceshouldleadtothe solu-tionofthemovinginfinitelinesource(movingpointsource).

The Neumannboundary condition for the lineheat source in thefrequencydomainis

−2

π

r

λ∂

Tˆ

(

r,

θ

,

ω

)

∂

r





r→0 = Qˆ in

(

ω

)

(41)

Similar toderiving Eq.(34),Eq.(41) leadsto a Fouriercosine seriesoftheform



nancos n

θ

= f

(

θ

)

|

r→0; f

(

θ

)

= αQˆin(ω)eUrcosθ/2α

πλ(Kn(qr)Urcosθ+2Kn+1(qr)αqr−2Kn(qr)αn)

(42)

To determine the a 0 coefficient, we solve Eq. (42) for n =0,

yielding

f

(

θ

)

=

α

Qˆ in

(

ω

)

eUrcosθ/2α

πλ

(

K0

(

qr

)

Urcos

θ

+ 2 K1

(

qr

)

α

qr

)

(43) As r →0,thisequationleadsto

lim

r→0f

(

θ

)

=

ˆ

Qin

(

ω

)

2

λπ

(44)

wheretheselimitingformsareutilized[23]: lim

r→0rK0

(

qr

)

Ucos

θ

∼ 0

lim

r→02 K1

(

qr

)

α

qr ∼ 2

α

(45)

Followingthis,the a 0 coefficientreads a0 = Qˆ₂in_λπ(ω)₂1_π  2π 0 d

θ

= Qˆin(ω) 2λπ (46)

TakingthelimitofEq.(42)as r →0,leadsto

an = 0 (47)

wherethislimitingcasehasresultedfrom: lim

r→02 Kn

(

qr

)

α

n∼ ∞ (48)

HavingtheFouriercoefficientsEqs.(46)and(47),thesolution tothemovinginfinitelinesourcecanthenbeexpressedas

u

(

r,

θ

,t

)

=  m

ˆ

Qin(ωm)

2λπ K0

(

qr

)

eiωmt

T

(

r,

θ

,t

)

= u

(

r,

θ

,t

)

e−Urcosθ/2α (49)

Inthe same way,the limitingcasefor theDirichlet boundary conditionyields u

(

r,

θ

,t

)

=  m− ˆ Tin(ωm) ln(q a)K0

(

qr

)

eiωmt T

(

r,

θ

,t

)

= u

(

r,

θ

,t

)

e−Urcosθ/2α

(

qa

)

is a small argument (50) wheretheselimitingformsareutilized[1]:

lim z→0K0

(

z

)

∼ − ln

(

z

)

lim z→0I0

(

z

)

∼ 1

z is a small argument (51)

5. Modelverificationandapplication

Themodelisverifiedagainstanalyticalsolutionsandnumerical solutions,viz.:

1 Comparingthe spectralmoving infinitecylindricalheat source modeltoitslimitingcasesandacommonlyusedGreen’s func-tioninfinitelineheatsourcemodelbymakingtheradiusofthe cylinderverysmall.

2 Comparingthe spectralmoving infinitecylindricalheat source modeltoafiniteelementmodel.

A two-dimensional, fully saturated porous domain, constitut-ing a solid phase and a water phase, subjected to a cylindrical heat source isconsidered forthispurpose. The geometry resem-bles a soil layer with groundwater flow passing a borehole heat exchanger.Thedomainisinlocalthermalequilibriumandtheheat equationisaveragedovertheconstituents,suchthat

∂

2_T

∂

x2 +

∂

2_T

∂

y2 + Ueff

α

eff

∂

T

∂

x − 1

α

eff

∂

T

∂

t = 0 (52)

R. Al-Khoury, N. BniLam, M.M. Arzanfudi et al. International Journal of Heat and Mass Transfer 163 (2020) 120517

Table 1

Physical and thermal parameters.

Parameter Magnitude Unit

Solid mass density ( ρs )

Solid specific heat ( c s )

Solid thermal conductivity ( λs )

Water mass density ( ρw )

Water specific heat ( c w )

Water thermal conductivity ( λw )

Porosity ( ϕ) Water velocity Heat source radius ( a )

2650 900 2.5 1000 4180 0.56 0.2

variable: 1e-5; 1e-4; 5e-4 variable: 0.0001; 0.075 kg/m 3 J/kg.K W/m.K kg/m 3 J/kg.K W/m.K - m/s m where

α

eff = _ρλceffeff

λ

eff =

ϕ

λ

w +

(

1 −

ϕ

)

λ

s

ρ

ceff =

ϕ

ρ

wcw +

(

1 −

ϕ

)

ρ

scs (53) and Ueff= n

ρ

wcw

ρ

ceff U (54)

where the subscripts w ,s ,eff denote the water phase, the solid phaseandtheeffectiveparameter,respectively,and

ϕ

isthe poros-ity. The thermo-physical parameters of the porous domain are giveninTable1.

Theinitialandboundaryconditionsare:

Tinitial= 0 ◦C ; t = 0

Tin= 10 ◦C ; t>0 for Dirichlet case Qin = 100 W /m ; t>0 for Neumann case

(55)

5.1. Verification against analytical solutions

Computationalresultsobtainedfromthe spectral moving infinite cylindrical source (SMICS) modelandits limitingcase:the spectral moving infinite line source (SMILS) model, are compared to each other andtothose obtainedfromthe Green’s function moving in- finite line source (GMILS)model,givenbySuttonetal.[22].

To mimic the line source, the radius of the cylinder in the SMICS model is made too small, namely 0.0001m. Both Dirich-letandNeumannboundaryconditionsaretested.Inliterature,the GMILSmodelisgivenforaNeumannboundaryconditiononly.To compare it with the spectral models for the Dirichlet boundary conditioncase,wefirstrun theGMILSmodelandoutput its tem-perature at0.0001mfromthelinesource,thenwe prescribethis temperatureonthecylindricalheatsourceintheSMICSmodel.

Fig.2showsthetemperaturecontoursfortheDirichlet bound-aryconditioncase,forthreewaterflow velocities:1e-5,1e-4and 5e-4m/s,andtwophysicaltimes:10and100hr.Thefigureshows nearly perfectmatchingbetweenthe threemodelsexceptforthe case with relatively fast fluid velocity and longer time (Fig. 2c, 100 h), where theGMILS modelexhibited some spurious oscilla-tions.

Fig.3showsthetemperaturecontoursfortheNeumann bound-ary conditioncasefortheabove mentioned waterflowvelocities andphysical times.Thethree modelsare ingood agreement, ex-cept fortherelatively fastfluidvelocity andlongertime,(Fig.3c, 100 h), where theGMILS modelexhibited some spurious oscilla-tions,similartothatfortheDirichletcase.

The three models were implemented in MATLAB, and while conducting the calculationsthe spectral models haveexhibited a robust capabilityinhandlinga widerrangeoffluidvelocitiesand physical times.Therunswere conductedfora stepfunction, and thus there was no significant CPU time difference between the

three models. However, if a random input signal has been em-ployed,therewouldhavebeenasignificantdifferencebetweenthe spectral models and the Green’s function model. Using the FFT, the spectral models are capable ofdealing withanyinput signal inthe sameefficiency asthat fora stepfunction; whereas using theGreen’sfunction wouldrequireatemporalsuperposition pro-cedurethat entailsdiscretizingtheinput signalintomultiplestep functions,foreach theinvolvedintegralsneed tobesolved sepa-rately.

5.2. Verification against numerical solutions

To examine the accuracy of the spectral model in simulat-ing a typicalcylindrical heatsource asexisting inGSHP systems, computationalresultsobtainedfromtheSMICSmodelwere com-paredto those obtained from thefinite element package: Multi-physics, Multidomain, Multiphase finite element analysis (MMM-FEM); a comprehensive 3D thermo-hydrodynamic-mechanical fi-nite element software for engineering applications in geother-malenergy,CO2 geo-sequestration,soil freezingandthawingand

poromechanicsin general;developed at DelftUniversity of Tech-nology[4–6].

Thefiniteelementdomainisahalfcylinder,20minradius, en-compassingacylindricalheatsourcewithradius0.075m,atypical radius of a boreholeheat exchanger. The mesh consistsof 9900, 8-node hexahedron elementswith a minimum size of 0.0026 m along the heat source circumference. The top view of the mesh withazoomaroundtheheat sourceisshowninFig.4.Fluid ve-locities: U =1e − 5m/s,and U =1e − 4m/swereexamined.

Fig.5showsthetemperaturecontoursafter48hforthe Dirich-let case. The figureclearlyshows avery good matchingbetween theSMICSmodelandthefiniteelementmodel.

ThegoodmatchingintheDirichletboundaryconditioncaseis duetotheexactdescriptionoftheboundaryconditionatthe ele-ments nodesdefiningthe heatsource.Thiscannot be guaranteed fortheNeumannboundaryconditioncaseastheheat fluxinthe finiteelementmethodisbydefinitioninterelement-discontinuous. Theheatfluxisdiscretizedattheelementsintegrationpointsand then mappedto the nodes;leadingeventually to an errorin de-scribingthe boundary condition. Tocircumvent thisproblem, we first conductedspectral analysis fora Neumann boundary condi-tion and output the temperature profile at the source boundary (a =0.075m) versustime.Thistemperatureprofilewasthen pre-scribedatthesourceboundaryinthefiniteelementmodel.Fig.6 showsthetemperaturecontoursafter48hr.Itshowsthatthereis agood matchingbetweenthe SMICSmodelforNeumann bound-arycondition andthe equivalentDirichlet casecomputed by the finiteelementmodel.

6. ComparingspectralmodeltoLaplacetransformandGreen’s functionsolutions

Existing solution techniques for a moving cylindrical heat source in a convective-conductive domain, applicable mainly to groundsourceheatpumpsystems,canbeputintotwocategories: LaplacetransformandGreen’sfunction.Inthissectionwepresent theLaplacetransformsolutiongivenbyCarslawandJaeger[9]and theGreen’sfunctionsolutiongivenbyZhangetal.[26]tohighlight thedifference in complexityofthe mathematicalformulations of thethreemodels.

6.1. Laplace transform solution

Carslaw andJaeger provided, in their seminal monograph on Heat Conduction in Solids [9], an analytical solution to an infi-nite solid initially at constant temperature T 0, moving along the

R. Al-Khoury, N. BniLam, M.M. Arzanfudi et al. International Journal of Heat and Mass Transfer 163 (2020) 120517

Fig. 3. Temperature contour plots of SMICS and SMILS models vs. GMILS model for Neumann boundary condition

Fig. 4. Top view of the finite element mesh with a zoom around the heat source

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,theysolvedEq.(1)as

T

(

r,

θ

,t

)

= u

(

r,

θ

,t

)

e−Urcosθ/2α_{+ T} 0 (56) with u

(

r,

θ

,t

)

= −T 0 ∞  n=0 εnIn(Ua/2α)Kn(Ur/2α) Kn(Ua/2α) cos n

θ

−2T0 π ∞  n=0 e−U2t/4α

_ε

ncos n

θ

In



_Ua 4α



×  _∞ 0 e−αu2t_[J n(ur)Yn(ua)−Yn(ur)Jn(ua)]udu [ u2_+(U2_/4α2_{)] [ J}2 n(ua)+Yn2(ua)]

ε

0 =1 if n=0 ,

ε

n =2 if n≥ 1 (57)

in which Y n isthe Besselfunction ofthesecond kindof order n .

Allotherparametersaredefinedabove.

Obviously, the presence of the semi-infinite integral of this transcendentalfunctionwithoscillatoryBesselfunctionsmakesthe problemcomplicatedandthecalculation,ifpossible,mustbe con-ducted using numerical algorithms. Carslaw and Jaeger [9] pro-vidednonumericalresultstothiscase,neithersimplifiedthe func-tion for short or long times,as they usually did formany other cases.Theyalsodidn’tprovideasolutionfortheNeumann bound-aryconditioncase.

6.2. Green’s function solution

Zhang et al. [26] provided an analytical solution to the infi-nite and finite cylindrical heat sources in convective-conductive domains.TheirsolutionisbasedontheGreen’sfunction formula-tionofaninstantaneous pointsourcegivenbyCarslawandJaeger [9]. Due to the unsymmetrical nature of heat convection arising fromgroundwaterflow,heatemittedfromapointisdescribedby integrating over the azimuth of the point. To consider the com-binedeffectsofallpointsconstitutingthecylinder,thesolutionof apoint isintegratedoverthe polarangles ofthecylinder.As the heatsourceismovingwithtimewitha certainvelocity,the solu-tioniscompletedbyintegrationovertime.Fig.7shows schemati-callyZhangetal.conceptualmodel.

Fortheinfinitecylindricalsourcecase, thismodelisdescribed byatripleintegral,suchthat

T(t) = Qin 8π3_λ 2π  0 dϕ 2_π 0

exp U(acosφ−acosϕ)

2α

 ·

4ατ/ ( acosφ−acosϕ₎ 2

+( asinφ−asinϕ₎ 2

 0 1 η× exp  −1 η−

U2 ( _{acosφ−acosϕ)} 2

+( asinφ−asinϕ₎ 2

16α2



dφdη

(58) inwhich

φ

isthepolarangleofthecylindricalheatsource,and

ϕ

 isthe polarangleofthe temperaturedistributionaround a point source.

Calculatingthisequationrequires:

1 Apropernumericalintegrationscheme.Though, suchintegrals might be restrictive, especially for far distances and longer times.

2 Q in is a step function, but ifit is a function of time, as it is

usually the case, its time signal must be divided into several stepfunctions,foreach,thisequationmustberecalculated. 3 Fordetailedanalysisto producecontour plotsoftemperature,

this equation can take rather long CPU time. Zhang et al. [26]providednonumericalresultsforaconvectivecase show-ing isothermal contours for an x – y plane; rather providing contourlinesinthe z _−directionfortwoanglesonly:

φ

₌0and

φ

=

π

.

Comparing these two solution techniques to those given in Eq.(29),fortheDirichletcase,andEqs.(34)-(36)fortheNeumann,

R. Al-Khoury, N. BniLam, M.M. Arzanfudi et al. International Journal of Heat and Mass Transfer 163 (2020) 120517

Fig. 6. Contour plots of Neumann SMICS model vs. equivalent Dirichlet MMM-FEM

Fig. 7. A schematic representation of Zhang et al. [26] conceptual model

it can readily be deduced that the spectral solution has a rather simplerformulationwhichmakes itcomputationally efficientand robust.

7. Conclusions

Heat conduction-convection in porous solids occurs in many engineering applications and technologies, including the ground sourceheatpumpsystems;atechnologydesignedtoharvest ther-mal energy from shallow depths of the earth to be utilized for heatingand coolingof buildings.In thepresence ofgroundwater flow,theamountofharvestedenergydependssignificantlyonthe amountofheatthatistransportedbythegroundwater,andthusit isimportanttobeconsideredinthedesignandanalysisofthe sys-tem. This paperaddressesthis issueandintroduces an analytical solution toconductive-convective heat flowinan infinitedomain subjectedtoa cylindricalheatsource.Thespectral analysisis uti-lizedtosolvethegoverningheatequationforprescribed Dirichlet andNeumannboundaryconditions.

Thesolution fortheDirichletboundary conditionhasledto a closedformfunction,andthat fortheNeumannboundary condi-tion has led to an integral which is relatively easy to solve nu-merically.ComparedtosolutionsbasedontheGreen’sfunctionor theLaplacetransform,thespectralsolutionhasasimpler formula-tion,asevident inSection6.Thismakesthespectral model com-putationallyefficientandrobust.Additionally,thespectralanalysis isinherently applicableto anyrandominput signal andits series summation (rather than integration) makes the solution less re-strictiveonthemodelphysicalandthermalparameters.

Thespectralmodelpresentedinthispapercanbeincorporated within the spectral element method [13] and the superposition principletostudydetailedheat flowingroundsourceheatpump systemsconstitutingmultipleboreholeheatexchangersembedded inmultilayersystemswithlayersencompassinggroundwaterflow. Inparticular,thegivenspectralmodelcanreadilybeincorporated in the spectral element model of BniLam et al. [7]. This will be undertakeninafuturework.

Authorshipcontributions

Pleaseindicate thespecificcontributionsmadeby eachauthor (listtheauthors’initialsfollowedbytheirsurnames,e.g.,Y.L. Che-ung).Thenameofeachauthormustappearatleastonceineach ofthethreecategoriesbelow.

DeclarationofCompetingInterest None.

Acknowledgment

Theauthors acknowledgetheinsightofProfessorSorinPopof theUniversityofHasseltonEq.(5).

References

[1] M. Abramowitz , I.A. Stegun , Handbook of Mathematical Functions, Dover Pub- lications, 1972 .

[2] R. Al-Khoury , Computational Modeling of Shallow Geothermal Systems, CRC Press, 2012 .

[3] R. Al-Khoury , P.G. Bonnier , R. Brinkgreve , Efficient finite element formulation for geothermal heating systems. Part I: steady state, Int. J. Numer. Methods Eng. 63 (2005) 988–1013 .

[4] M.M. Arzanfudi , R. Al-Khoury , Thermo-hydrodynamic-mechanical multiphase flow model for crack propagation in deformable porous media, Int. J. Numer. Methods Fluids 84 (2017) 635–674 .

[5] M.M. Arzanfudi , R. Al-Khoury , Freezing-thawing of porous media: An extended finite element approach for soil freezing and thawing, Adv. Water Res. 119 (2018) 210–226 .

[6] M.M Arzanfudi , R. Al-Khoury , L.J. Sluys , G.M.A. Schreppers , A thermo-hy- dro-mechanical model for energy piles under cyclic thermal loading, Comput. Geotech. 125 (2020) 103560 .

[7] N. BniLam , R. Al-Khoury , A. Shiri , L.J. Sluys , A semi-analytical model for de- tailed 3D heat flow in shallow geothermal systems, Int. J. Heat Mass Transfer 123 (2018) 911–927 .

[8] A. Capozza , M. De Carli , A. Zarrella , Investigations on the influence of aquifers on the ground temperature in ground-source heat pump operations, Appl. En- ergy 107 (2013) 350–363 .

[9] H.S. Carslaw , J.C. Jaeger , Conduction of Heat in Solids, second ed., Oxford Uni- versity Press, 1959 .

[10] M. Cimmino , B.R Baliga , A hybrid numerical-semi-analytical method for com- puter simulations of groundwater flow and heat transfer in geothermal bore- hole fields, Int. J. Therm. Sci. 142 (2019) 366–378 .

[11] N.R. DesRuisseaux , R.D. Zerkle , Temperature in semi–infinite and cylindrical bodies subjected to moving heat sources and surface cooling, J. Heat Transfer 92 (1970) 456–464 .

[12] N. Diao , Q. Li , Z. Fang , Heat transfer in ground heat exchangers with ground- water advection, Int. J. Therm. Sci. 43 (2004) 1203–1211 .

[13] J.F. Doyle , Wave Propagation in Structures: Spectral Analysis Using Fast Dis- crete Fourier Transforms, Springer-Verlag, New York, 1997 .

[14] S. Erol , B. Francois , Multilayer analytical model for vertical ground heat ex- changer with groundwater flow, Geothermics 71 (2018) 294–305 .

[15] MATLAB®. https://www.mathworks.com/help/matlab/ref/rand.html . [16] N. Molina-Geraldo , P. Blum , K. Zhu , P. Bayer , Z. Fang , A moving finite line

source model to simulate borehole heat exchangers with groundwater advec- tion, Int. J. Therm. Sci. 50 (2011) 2506–2513 .

[17] M. Nabi , R. Al-Khoury , An efficient finite volume model for shallow geother- mal systems. Part I: model formulation, Comput. Geosci. 49 (2012) 290–296 .

[18] M. Nabi , R. Al-Khoury , An efficient finite volume model for shallow geothermal systems—Part II: verification, validation and grid convergence, Comput. Geosci. 49 (2012) 297–307 .

[19] Y. Nam , R. Ooka , S. Hwang , Development of a numerical model to predict heat exchange rates for a ground-source heat pump system, Energy Build. 40 (2008) 2133–2140 .

[20] J.A. Rivera , P. Blum , P. Bayer , Analytical simulation of groundwater flow and land surface effects on thermal plumes of borehole heat exchangers, Appl. En- ergy 146 (2015) 421–433 .

[21] L.F. Shampine, Vectorized adaptive quadrature in MATLAB®, J. Comput. Appl. Math. 211 (2008) 131–140 See also https://nl.mathworks.com/help/matlab/ref/ integral.html .

[22] M.G. Sutton , D.W. Nutter , R.J. Couvillion , A ground resistance for vertical bore heat exchangers with groundwater flow, J. Energy Res. Technol. 125 (2003) 183–189 .

[23] Wolframalpha: https://www.wolframalpha.com/

[24] B. Weigand , Analytical Methods for Heat Transfer and Fluid Flow Problems (2015) .

[25] C. Yavuzturk , J.D. Spitler , S.J. Rees , A transient two-dimensional finite volume model for the simulation of vertical U-tube ground heat exchangers, ASHRAE Trans. 105 (2) (1999) 527–540 .

[26] W. Zhang , H. Yang , L. Lu , Z. Fang , The analysis on solid cylindrical heat source model of foundation pile ground heat exchangers with groundwater flow, En- ergy 55 (2013) 417–425 .

[27] Y. Zhou , C. Xu , D. Sego , D. Zhang , Analytical solution for solid cylindrical heat source model with convective boundary condition, J. Heat Transfer 141 (2019) 121701-1–121701-14 .

[28] S.M. Zubair , M.A. Chaudhry , Temperature solutions due to time-dependent moving-line-heat sources, Heat Mass Transf. 31 (1996) 185–189 .