An analytical solution for the temperature field around a cylindrical surface subjected to a time dependent heat flux

International Journal of Heat and Mass Transfer 66 (2013) 906–910

Contents lists available at ScienceDirect

International Journal of Heat and Mass Transfer journal homepage: www.elsevier.com/locate/ijhmt

An analytical solution for the temperature field around a cylindrical surface subjected to a time dependent heat flux Enzo Zanchini ⇑, Beatrice Pulvirenti Universitá di Bologna, Dipartimento di Ingegneria Industriale (DIN), Viale Risorgimento 2, 40136 Bologna, Italy

a r t i c l e

i n f o

Article history: Received 19 June 2013 Accepted 25 July 2013

Keywords: Unsteady heat conduction Cylindrical geometry Time-dependent boundary condition Analytical solution

a b s t r a c t The analytical solution for the temperature field in an infinite solid medium which surrounds a cylindrical surface, determined by Carlslaw and Jaeger for the case of constant heat flux, is extended to the case of any time dependent heat flux. Then, with reference to a sinusoidally varying heat flux, the analytical solution is employed to determine benchmark results for the time evolution of the dimensionless temperature of the surface. These results are used to check the accuracy of the numerical solutions obtained by two different commercial codes: the finite-element code COMSOL Multiphysics and the finite-volume code FLUENT. Ó 2013 Elsevier Ltd. All rights reserved.

1. Introduction An important transient heat conduction problem is the determination of the unsteady temperature field produced by a heated or cooled cylindrical surface surrounded by an infinite solid medium. An analytical solution of this problem, for the case of uniform and constant heat flux per unit area, has been presented by Carlslaw and Jaeger [1], and has been widely used for the design of Borehole Heat Exchanger (BHE) fields for ground-coupled heat pump systems. Indeed, the design method for BHE fields recommended by ASHRAE [2] and developed by Kavanaugh and Rafferty [3] is based on this solution. In fact, it sketches each BHE as a cylindrical heat source subjected to a uniform and constant heat flux per unit area, and considers the superposition of three heat pulses, each with a constant power, with different durations: 6 h, 1 month and 10 years. This method is still widely used, but recent studies showed that, if the seasonal heat loads are unbalanced and the effects of groundwater movement are negligible, large BHE fields can reach critical conditions after some decades [4,5], so that an analysis for a period of 10 years may be insufficient. In the study of the long-term sustainability of BHE fields, time dependent heat loads must be considered. A time dependent heat load can be either obtained as a weighted superposition of constant heat loads displaced in time, or directly sketched by a suitable function of time. The second choice was performed, for instance, in Ref. [5], where the long-term performance of infinitely

⇑ Corresponding author. Tel.: +39 0512093295; fax: +39 0512093296. E-mail addresses: [email protected] (E. Zanchini), beatrice.pulvirenti@ unibo.it (B. Pulvirenti). 0017-9310/$ - see front matter Ó 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.ijheatmasstransfer.2013.07.092

long BHE lines placed in a ground with negligible groundwater movement was studied numerically by means of the finite element software package COMSOL Multiphysics. Heat loads with period of 1 year were considered, varying with time with a sinusoidal law during each season, with different degrees of compensation of winter heating with summer cooling. The numerical code was validated in the special case of constant heat load, by comparison with the analytical solution of Carlslaw and Jaeger [1]. Indeed, an analytical solution for the temperature field around a cylindrical surface subjected to a uniform but time dependent heat flux and surrounded by an infinite solid medium would be useful, to validate the results obtained by means of numerical simulation codes, but is not directly available in the literature. Two analytical solutions of the Fourier equation for the temperature field around a cylindrical surface placed in an infinite solid medium were presented, but both consider the case of a uniform and constant heat flux per unit area applied to the surface. The first is the well known solution by Carlslaw and Jaeger [1]; the second, reported in a book in German [6], is recalled in a more recent paper in English [7]. An analytical solution of the Cattaneo–Vernotte hyperbolic heat conduction equation in cylindrical geometry, with the boundary condition of a time variable heat flux at the inner surface, was determined by Barletta [8]. However, the author did not present or discuss the limit of his solution in the case of a vanishing relaxation time. Since the problem, especially in the hyperbolic heat conduction case, is very far from being elementary, the development of an independent solution of the Fourier equation, with the same geometry and boundary conditions, seems useful, to build a bridge between the simplest case, studied in Refs. [1,6],

907

E. Zanchini, B. Pulvirenti / International Journal of Heat and Mass Transfer 66 (2013) 906–910

and the most complex case considered in Ref. [8], and to provide a cross validation of the available solutions. In this paper, an analytical solution of the Fourier equation in an infinite solid medium bounded internally by a cylindrical surface subjected to any uniform and time-dependent heat flux, is presented. It is shown that the solution reduces to that determined by Carlslaw and Jaeger [1] in the special case of constant heat flux, and coincides with the limit for vanishing relaxation time of the solution determined by Barletta [8] for hyperbolic heat conduction. Then, for the special case of a sinusoidally varying heat flux, the analytical solution is employed to determine a benchmark table of values of the dimensionless temperature and to check the accuracy of numerical solutions obtained by two different commercial codes. 2. Analytical solution Let us consider a homogeneous and infinite solid medium, bounded internally by the cylindrical surface r = r0, i.e., which occupies the whole region of space r0 6 r < + 1. Let us assume that the thermal conductivity k and the thermal diffusivity a of the medium are constants and no heat generation is present within the solid. At the initial instant of time, s = 0, the temperature field within the solid is uniform, with a value T0, and stationary. For s > 0, a uniform and time-dependent heat flux per unit area q(s) = q0F(s) is applied to the solid, at the internal surface r = r0, where F(s) is a dimensionless function of time. Under these conditions, the temperature field in the solid is axisymmetric, and the differential equation for heat conduction in a cylindrical coordinate system can be written as

@ 2 T 1 @T 1 @T þ ¼ ; @r 2 r @r a @ s

k

 @T  ¼ q0 FðsÞ; @r r¼r0

s > 0:

r ; r0

n¼

as

Tðr; sÞ  T 0 : q0 r0

~ ¼  /ðsÞ:

ð12Þ

g¼1

The general solution of Eq. (11) is

pffiffi pffiffi ~ g; sÞ ¼ c1 ðsÞI0 ðg sÞ þ c2 ðsÞK 0 ðg sÞ; hð

where c1(s) and c2(s) are arbitrary functions of s, while I0 and K0 are the modified Bessel functions of the first and second kind, with order 0. Since the dimensionless temperature must be vanishing for g ? 1, the asymptotic properties of Bessel functions imply that c1(s) = 0. By applying the boundary condition Eq. (12), and the relation K0 0(z) =  K1(z), we obtain the following solution of Eq. (11), pffiffi ~ ~ g; sÞ ¼ /ðsÞKp0 ffiffiðgpffiffisÞ ; hð K 1 ð sÞ s

ð14Þ

where K1 is the modified Bessel function of the second kind with order 1. On account of the convolution theorem, the inverse Laplace transform of Eq. (14) has the form

hðg; nÞ ¼

Z

n

vðg; uÞ/ðn  uÞdu;

ð15Þ

0

where

vðg; uÞ ¼ L1 fAðg; sÞg ¼

1 2p i

Z

cþi1

esu Aðg; sÞds;

ð3Þ

Z

cþi1

ci1

ð17Þ

I Z esu Aðg; sÞds ¼ lim lim esu Aðg; sÞds  esu Aðg; sÞds R!þ1 !0 C CR Z Z esu Aðg; sÞds  esu Aðg; sÞds  C AB  Z ð18Þ esu Aðg; sÞds :  CD

ð5Þ

@ 2 h 1 @h @h þ ¼ ; @ g2 g @ g @n

ð6Þ

hðg; 0Þ ¼ 0;

ð7Þ

with

/ðnÞ ¼ F



n > 0;

 : 2

as r0

Since there is no pole of A(g, s) within the region bounded by C, the integral on C is zero, on account of Cauchy integral theorem. Byffi pffiffiffiffiffiffiffiffi considering the asymptotic expression of K n ; K n ðsÞ  es = 2ps,

ð8Þ

ð9Þ

Eq. (6) can be solved by using the Laplace transform of h(g,n) with respect to n,

~ g; sÞ ¼ hð

Z

þ1

esn hðg; nÞdn:

ð10Þ

0

In the Laplace transformed domain, Eq. (6) becomes

ð16Þ

ci1

Thus, Eqs. (1)–(3) can be rewritten as

 @h  ¼ /ðnÞ; @ gg¼1

ð13Þ

The integral in Eq. (16) can be evaluated by considering the closed path C shown in Fig. 1 as

and the dimensionless temperature

hðg; nÞ ¼ k

 @ h~  @ g

ð2Þ

ð4Þ

r 20

with the boundary condition

pffiffi K 0 ð g sÞ p ffiffi pffiffi : Aðg; sÞ ¼ K 1 ð sÞ s

Let us introduce the dimensionless coordinates

g¼

ð11Þ

ð1Þ

with initial and boundary conditions given by

Tðr; 0Þ ¼ T 0 ;

@ 2 h~ 1 @ h~ þ  s h~ ¼ 0; @ g2 g @ g

Fig. 1. Integration path.

908

E. Zanchini, B. Pulvirenti / International Journal of Heat and Mass Transfer 66 (2013) 906–910

one can prove that the integral on CR vanishes for R ? 1. Moreover, A(g, s) and the integral on Ce vanish for e ? 0. Then, Eq. (16) becomes

pffiffi  Z  1 K 0 ðg sÞ pffiffi pffiffi ds : vðg; uÞ ¼  lim lim 2 esu 2pi R!þ1 !0 K 1 ð sÞ s AB

p 2

ðiÞ

nþ1

½J n ðzÞ  iY n ðzÞ;

ð20Þ

where Jn is the Bessel function of the first kind with order n, one obtains the following expression of the real part of v(g,u):

vðg; uÞ ¼

1

Z

p

þ1

Bðg; yÞeyu dy;

ð21Þ

0

where B(g, y) is the function

Bðg; yÞ ¼

pffiffiffi pffiffiffi pffiffiffi pffiffiffi Y 0 ðg yÞJ 1 ð yÞ  J 0 ðg yÞY 1 ð yÞ : 2 pffiffiffi 2 pffiffiffi pffiffiffi ½J 1 ð yÞ þ Y 1 ð yÞ y

ð22Þ

Thus, the solution of Eqs. (6)–(8) is

hðg; nÞ ¼

Z

1

p

þ1

Bðg; yÞ

Z

0

hðg; nÞ ¼

Z

2

eyu /ðn  uÞdu dy:

p

ð23Þ

0

þ1

Cðg; zÞ 0

Z

n

pffiffiffi y, one can rewrite the solution as

2

ez u /ðn  uÞdu dz;

ð24Þ

0

where C(g, z) is the function

Cðg; zÞ ¼

Y 0 ðgzÞJ 1 ðzÞ  J 0 ðgzÞY 1 ðzÞ J 21 ðzÞ þ Y 21 ðzÞ

ð25Þ

:

hðg; nÞ ¼

Z

p

þ1

Cðg; zÞ 0

Z

n

z2 u

e

du dz:

ð26Þ

0

By solving the second integral one obtains the well known solution [1]

hðg; nÞ ¼

Z

2

p

þ1 0

Y 0 ðgzÞJ 1 ðzÞ  J 0 ðgzÞY 1 ðzÞ 1  ez z2 J 21 ðzÞ þ Y 21 ðzÞ

dz:

s0

as0 r 20

ð28Þ

:

ð29Þ

;

we can rewrite F(s) in the dimensionless form

/ðnÞ ¼ sin

  2p n ; Fo

ð30Þ

where Fo takes the role of a dimensionless time period. With this boundary condition, the general solution Eq. (24) becomes

hðg; nÞ ¼

2

p

Z

þ1

Cðg; zÞ 0

Z 0

n

2p e

pn pn  2pcos 2Fo þ Fo z2 cos 2Fo

4p2 þ Fo2 z4

dz:

ð32Þ

Values of the dimensionless temperature h yielded by Eq. (32), for

g= 1, Fo = 5000 and n/Fo (number of periods) variable in the range from 0 to 10, were determined by means of the software package Mathematica. To obtain results as accurate as possible, the integral which appears in Eq. (32) was evaluated numerically only in the range 0 6 z 6 1000, while an analytical expression of the integral was determined in the range 1000 6 z < + 1. In fact, for g = 1 and z > 1000, the value of C(1, z) is very close to 1 (0.999999 for z = 1000) and the integral can be solved analytically by assuming C(1, z) = 1. In the range 1000 6 z < +1 the first term which appears at the numerator of the fraction in Eq. (32) is negligible, and the relevant solution is

Z

þ1

pn pn 2pcos 2Fo þ Fo z2 cos 2Fo

dz 4p2 þ Fo2 z4 " ( rffiffiffiffi#   1=4  3=4  1 1 1 2 2pn 1 2 cos ¼  pffiffiffi  2 coth 500  2 Fo Fo p 2 2 Fo Fo v 2 3 q ffiffi ffi u u 500 p2 7 2pn u 1 6 þ tFocoth 4

1=4 5sin Fo  Fo12 9 2 3 qffiffiffi "   rffiffiffiffiffiffiffiffiffiffiffiffi  #> 2 500 2pn 1 2pn = p 7 6 þ cot 1 4

: þ  2 Fosin 1=4 5 cos > Fo Fo Fo ;  Fo12

2

ez u sin



By the method described above, benchmark values of h were obtained from Eq. (32), for g= 1, Fo = 5000 and n/Fo variable in the range from 0 to 10 with steps of 0.05. These values are reported in Table 1 for 0.05 6 n/Fo 6 5 and in Table 2 for 5.05 6 n/Fo 6 10. We checked that, for the number of digits reported in the table, the results did not change when the numerical integration param-

Table 1 Benchmark values of h in g = 1 for sinusoidal heat flux, for Fo = 5000 and n/Fo 6 5.

By introducing the Fourier number

Fo ¼

z2 n



Cðg; zÞ

0

ð27Þ

Let us now consider the special case in which the heat flux per unit area is a sinusoidal function of time, with period s0, i.e. the dimensionless function F(s) is given by

2ps

þ1

2n

3. Analytical solution and benchmark results for sinusoidal heat flux

FðsÞ ¼ sin

p

Z



ð33Þ

If the heat flux is constant, i.e. /(n) =1, Eq. (24) becomes

2

2Fo

1000

n

By introducing the variable z ¼

hðg; nÞ ¼

ð19Þ

On AB, one has s = eip y, with y real and non-negative. By employing the variable y = eip s and by considering the following property of Bessel functions

K n ðeip=2 zÞ ¼

By solving the second integral, one can rewrite Eq. (31) as

2p ðn  uÞ dudz: Fo

ð31Þ

n/Fo

h

n/Fo

h

n/Fo

h

n/Fo

h

0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0.55 0.6 0.65 0.7 0.75 0.8 0.85 0.9 0.95 1.0 1.05 1.1 1.15 1.2 1.25

0.8338 1.7957 2.6628 3.3218 3.6950 3.7376 3.4405 2.8288 1.9598 0.9161 0.2015 1.2854 2.2304 2.9451 3.3604 3.4364 3.1662 2.5770 1.7269 0.6995 0.4041 1.4756 2.4097 3.1147 3.5213

1.3 1.35 1.4 1.45 1.5 1.55 1.6 1.65 1.7 1.75 1.8 1.85 1.9 1.95 2.0 2.05 2.1 2.15 2.2 2.25 2.3 2.35 2.4 2.45 2.5

3.5895 3.3122 2.7166 1.8605 0.8277 0.2809 1.3570 2.2954 3.0044 3.4147 3.4863 3.2124 2.6197 1.7666 0.7365 0.3696 1.4433 2.3794 3.0862 3.4945 3.5642 3.2883 2.6939 1.8391 0.8074

2.55 2.6 2.65 2.7 2.75 2.8 2.85 2.9 2.95 3.0 3.05 3.1 3.15 3.2 3.25 3.3 3.35 3.4 3.45 3.5 3.55 3.6 3.65 3.7 3.75

0.3003 1.3754 2.3130 3.0212 3.4307 3.5017 3.2270 2.6338 1.7800 0.7494 0.3572 1.4313 2.3679 3.0752 3.4838 3.5539 3.2784 2.6843 1.8298 0.7984 0.3089 1.3838 2.3211 3.0290 3.4383

3.8 3.85 3.9 3.95 4.0 4.05 4.1 4.15 4.2 4.25 4.3 4.35 4.4 4.45 4.5 4.55 4.6 4.65 4.7 4.75 4.8 4.85 4.9 4.95 5.0

3.5090 3.2342 2.6407 1.7868 0.7560 0.3508 1.4252 2.3619 3.0693 3.4781 3.5483 3.2730 2.6791 1.8247 0.7934 0.3138 1.3886 2.3257 3.0335 3.4428 3.5133 3.2384 2.6448 1.7908 0.7599

E. Zanchini, B. Pulvirenti / International Journal of Heat and Mass Transfer 66 (2013) 906–910 Table 2 Benchmark values of h in g = 1 for sinusoidal heat flux, for Fo = 5000 and 5 < n/Fo 6 10. n/Fo

h

n/Fo

h

n/Fo

h

n/Fo

h

5.05 5.1 5.15 5.2 5.25 5.3 5.35 5.4 5.45 5.5 5.55 5.6 5.65 5.7 5.75 5.8 5.85 5.9 5.95 6.0 6.05 6.1 6.15 6.2 6.25

0.3469 1.4214 2.3582 3.0657 3.4746 3.5449 3.2696 2.6758 1.8214 0.7902 0.3169 1.3916 2.3287 3.0365 3.4457 3.5162 3.2412 2.6476 1.7935 0.7625 0.3444 1.4188 2.3557 3.0632 3.4722

6.3 6.35 6.4 6.45 6.5 6.55 6.6 6.65 6.7 6.75 6.8 6.85 6.9 6.95 7.0 7.05 7.1 7.15 7.2 7.25 7.3 7.35 7.4 7.45 7.5

3.5425 3.2673 2.6735 1.8192 0.7880 0.3191 1.3938 2.3308 3.0386 3.4477 3.5182 3.2432 2.6495 1.7954 0.7644 0.3425 1.4170 2.3539 3.0615 3.4704 3.5408 3.2656 2.6718 1.8175 0.7864

7.55 7.6 7.65 7.7 7.75 7.8 7.85 7.9 7.95 8.0 8.05 8.1 8.15 8.2 8.25 8.3 8.35 8.4 8.45 8.5 8.55 8.6 8.65 8.7 8.75

0.3207 1.3953 2.3324 3.0401 3.4492 3.5197 3.2446 2.6510 1.7969 0.7658 0.3411 1.4156 2.3525 3.0601 3.4691 3.5395 3.2643 2.6705 1.8163 0.7851 0.3219 1.3966 2.3336 3.0413 3.4504

8.8 3.85 8.9 8.95 9.0 9.05 9.1 9.15 9.2 9.25 9.3 9.35 9.4 9.45 9.5 9.55 9.6 9.65 9.7 9.75 9.8 9.85 9.9 9.95 10.0

3.5209 3.2458 2.6521 1.7980 0.7669 0.3400 1.4145 2.3515 3.0591 3.4681 3.5384 3.2633 2.6695 1.8153 0.7842 0.3229 1.3975 2.3345 3.0422 3.4513 3.5218 3.2467 2.6530 1.7989 0.7678

Fig. 2. Plot of h versus n/Fo for sinusoidal heat flux, with Fo = 5000 and 0 6 n/Fo 6 10.

909

In each case, a 2D computational domain was considered, obtained by drawing a circle with dimensionless radius equal to 1000, centered in the axis of the BHE, to which the circle occupied by the BHE (0 6 g 6 1) was subtracted. The condition of vanishing heat flux was imposed at the external boundary, where the temperature changes were checked to be quite negligible, and the dimensionless-time interval 0 6 n/Fo 6 10 was selected. For the numerical simulations performed by means of COMSOL Multiphysics, unstructured triangular meshes were chosen. Three kinds of meshes were considered. The meshes of the first kind, which will be called kind COMSOL 1, were built through the options free mesh parameters, global, and predefined mesh size ‘‘Normal’’, and by applying refinements until the number of triangular elements reached 448512. Five meshes of kind COMSOL 1 were obtained, with 1752, 7008, 28032, 112128 and 448512 triangular elements respectively. The meshes of the second kind, which will be called kind COMSOL 2, were built through the options free mesh parameters, global, and predefined mesh size ‘‘Fine’’, and by applying refinements until the number of triangular elements reached 528896. Five meshes of kind COMSOL 2 were obtained, with 2066, 8264, 33056, 132224 and 528896 triangular elements respectively. The meshes of the third kind, which will be called kind COMSOL 3, were built through the options free mesh parameters, global, and predefined mesh size ‘‘Extra fine’’, and by applying refinements until the number of triangular elements reached 608512. Four meshes of kind COMSOL 3 were obtained, with 9508, 38032, 152028 and 608512 triangular elements respectively. The time dependent solver was chosen, with the following accuracy parameters: relative tolerance 105, absolute tolerance 106. For each mesh, the root mean square deviation r from the values of the dimensionless temperature obtained by the analytical solution (Eq. (32)) was evaluated. The results are illustrated in Fig. 3, where, for each kind of mesh, r is plotted versus the natural logarithm of the number n on triangular elements. The figure shows an excellent convergence to the analytical values for all kinds of meshes considered. In particular, the meshes of kind COMSOL 2 and COMSOL 3 reached a mean square deviation from the analytical results lower than 5  106 with 33056 and 38032 triangular elements, respectively. The lowest value of r was obtained by the mesh of kind COMSOL 3 with 152128 elements: it is 2.74  106 and corresponds to a relative root mean square deviation of 1.21  106 with respect to the arithmetic mean of the absolute values of the analytical results, which is 2.2600.

eters ‘‘Minrecursion’’, ‘‘Maxrecursion’’ and ‘‘Workingprecision’’ in Mathematica were changed from default to higher values. A plot of the dimensionless temperature h in g = 1 for sinusoidal heat flux, Fo = 5000 and 0 6 n/Fo 6 10 is reported in Fig. 2. The plot shows that, with this boundary condition, the temperature field becomes time-periodic after 3–4 periods of the heat flux oscillation. Fo = 5000 is a typical value of the Fourier number for U-tube BHEs subjected to a time periodic heat load with period of one year. It corresponds to: r0 = 0.08 m, a = 106 m2/s,s0 = 3.1536  107 s (365 days). 4. Accuracy of the solutions obtained by a finite-element and a finite-volume code The benchmark results obtained through Eq. (32), reported in Tables 1 and 2, were employed to check the accuracy of numerical solutions of Eqs. (6)–(8) and (30), obtained by two different numerical codes: the finite-element code COMSOL Multiphysics and the finite-volume code FLUENT.

Fig. 3. Plots of the root mean square deviation r from the analytical values versus ln n, for meshes of kind COMSOL 1, COMSOL 2 and COMSOL 3, with Fo = 5000 and 0 6 n/Fo 6 10.

910

E. Zanchini, B. Pulvirenti / International Journal of Heat and Mass Transfer 66 (2013) 906–910

sults are illustrated in Fig. 4, where, for each kind of mesh, r is plotted versus the natural logarithm of the number n on elements. The figure shows an excellent convergence to the analytical values for all kinds of meshes considered. In particular, the meshes of kind FLUENT 2 and FLUENT 3 reached a mean square deviation from the analytical results lower than 104 with 128000 and 152000 elements, respectively. The lowest value of r was obtained by the mesh of kind FLUENT 2 with 512000 elements: it is 1.00  105 and corresponds to a relative root mean square deviation of 4.42  106 with respect to the arithmetic mean of the absolute values of the analytical results. 5. Conclusions

Fig. 4. Plots of the root mean square deviation r from the analytical values versus ln n, for meshes of kind FLUENT 1, FLUENT 2 and FLUENT 3, with Fo = 5000 and 0 6 n/Fo 6 10.

For the numerical simulations performed by means of FLUENT, structured quadrilateral meshes were chosen. Three kinds of meshes were considered. The meshes of the first kind, which will be called kind FLUENT 1, were built by dividing the cylindrical surface in equal parts, and by applying an exponential increase of the element size along the radial direction with coefficient 1.19. Five meshes of kind FLUENT 1 were obtained, with 1700, 6800, 27200, 108800 and 435200 elements respectively. The meshes of the second kind, which will be called kind FLUENT 2, were built by dividing the cylindrical surface in equal parts, and by applying a method named first length on the radial direction, with the first length equal to 0.08, 0.04, 0.02, 0.01 respectively in the five cases. Five meshes of kind FLUENT 2 were obtained, with 2000, 8000, 32000, 128000 and 512000 structured elements respectively. The meshes of the third kind, which will be called kind FLUENT 3, were built by dividing the cylindrical surface in equal parts, and by applying a method named successive ratio on the radial direction, with the successive ratio equal to 1.08, 1.04, 1.02, 1.01 respectively in the five cases. Four meshes of kind FLUENT 3 were obtained, with 9500, 38000, 152000 and 608000 structured elements respectively. The time dependent solver was chosen, with a dimensionlesstime step equal to 0.005. For each mesh, the root mean square deviation r from the values of the dimensionless temperature obtained by the analytical solution (Eq. (32)) was evaluated. The re-

The analytical solution for the temperature field in an infinite solid medium which surrounds a cylindrical surface, determined by Carlslaw and Jaeger, has been extended to the case of any time dependent heat flux. The analytical solution has been employed to determine benchmark values of the dimensionless temperature of the surface, in the case of sinusoidally varying wall heat flux. The benchmark results have been compared with those obtained through the finite-element code COMSOL Multiphysics and the finite-volume code FLUENT. The results show that both numerical codes, if employed with good choices of the meshes and of either the accuracy parameters (COMSOL) or the time step (FLUENT), yield results in excellent agreement with the analytical ones. The benchmark analytical results can be useful to optimize the choice of the mesh and of either the accuracy parameters or the time step of the numerical codes. References [1] H.S. Carlslaw, J.C. Jaeger, Conduction of Heat in Solids, Oxford University Press, Oxford, 1959. [2] 2007 ASHRAE Handbook, HVAC Applications, ASHRAE, Atlanta, 2007. Chapter 32. [3] S.P. Kavanaugh, K. Rafferty, Ground-Source Heat Pumps. Design of Geothermal Systems for Commercial and Institutional Buildings, ASHRAE, Atlanta, 1997. [4] S. Signorelli, T. Kohl, L. Rybach, Sustainability of production from borehole heat exchanger fields, in: World geothermal congress 2005, Antalya, Turkey, 24–29 April 2005. [5] S. Lazzari, A. Priarone, E. Zanchini, Long-term performance of BHE (borehole heat exchanger) fields with negligible groundwater movement, Energy 35 (2010) 4966–4974. [6] H. Tautz, Warmeleitung und Temperaturausgleich, Verlag Chemie GmbH Weinheim/Bergstr, 1971. [7] A. Graßmann, F. Peters, Experimental investigation of heat conduction in wet sand, Heat Mass Transfer 35 (1999) 289–294. [8] A. Barletta, Hyperbolic propagation of an axisymmetric thermal signal in an infinite solid medium, Int. J. Heat Mass Transfer 39 (1996) 3261–3271.