- Email: [email protected]
Horizontal ground heat exchangers modelling
Applied Thermal Engineering 155 (2019) 534–545
Contents lists available at ScienceDirect
Applied Thermal Engineering journal homepage: www.elsevier.com/locate/apthermeng
Research Paper
Horizontal ground heat exchangers modelling Louis Lamarche
T
École de technologie supérieure, 1100 Notre-Dame Streetr West, Montreal H3C 1K3, Canada
HIGHLIGHTS
model for horizontal ground heat exchanger. • New model which is very fast and accurate with respect to numerical simulation. • Analytical comparison with the classical trench model is performed. • AParticularity of horizontal systems are discussed possible applications of the model are discussed. • ARTICLE INFO
ABSTRACT
Keywords: Ground source heat pump Horizontal ground heat exchanger Finite line source
Geothermal heat exchangers applied to building air conditioning and heating systems are often categorized as vertical and horizontal configurations. In this article, an analytical model is presented to simulate on an hourly basis, ground heat exchangers having a horizontal configuration. The model is based on a new formalism of the finite line source approach that allows faster simulation. The model presents a novel approach to take into the effect of seasonal temperature variations especially when the inlet temperatures in the ground heat exchanger are imposed. It also takes into account the effect of the air-soil boundary condition and the thermal interference between the different portions of pipes in the buried field. The temperature as well as the heat flux distribution along the ground heat exchanger is part of the solution and not a priori assumed as in the case of previous models. The model does not consider the possible phase change due to ice formation or the moisture content variation associated with some horizontal configurations. It is based on recent studies related to vertical boreholes and adapted to horizontal systems. It is intended to be a valuable tool to simulate on an hourly basis, the time response of the ground heat exchanger coupled with a building equipped with heat pumps.
1. Introduction Ground-Coupled Heat Pumps (GCHP) using a secondary loop ground heat exchanger are often categorized as vertical and horizontal loop systems, although hybrid systems such as inclined or pile systems, are also available. Vertical borehole systems are the most popular technology in Europe, but in North America, horizontal systems continue to be largely used [1,2]. This observation notwithstanding, there is general agreement that the scientific literature on horizontal systems is much more sparse than that on vertical ones [3]. Design approaches often continue to be associated with rules-of-thumb and approximate charts [3]. Even though some common physical aspects link both technologies, there are some specific factors associated with horizontal loops that make the design of such systems unique. The following is a non-exhaustive list of such factors: 1. The ground surface has a greater influence on the thermal response
than for vertical systems. 2. Seasonal temperature variations, which are often neglected in vertical systems, are very important in horizontal ones. 3. Horizontal piping lay in trenches in the subsurface overburden where the moisture content can influence the thermal performance of the heat exchanger. 4. Phase change due to freezing and thawing may occur. 5. Several configurations are possible, including inline, slinky, etc., which cause pipe-to-pipe interference. Given the above factors, horizontal loops are more complicated to analyse, even though they may be easier to install. Several studies have looked at these systems. A well-known report on horizontal systems [4] by Claesson and Dunand is often used as a reference even though several aspects shown in the list above not addressed. Mei [5] develop a numerical model for horizontal systems, which was extended by Piechowski [6] to include mass transfer due to moisture content. More
E-mail address: [email protected]. https://doi.org/10.1016/j.applthermaleng.2019.04.006 Received 5 December 2018; Received in revised form 20 February 2019; Accepted 2 April 2019 Available online 04 April 2019 1359-4311/ © 2019 Elsevier Ltd. All rights reserved.
Applied Thermal Engineering 155 (2019) 534–545
L. Lamarche
Nomenclature A Cp D E1 h H k I L m q q’ q” r R Sp Sq t T u
X Y x−y−z
Area (m2) Specific heat capacity (J-Kg−1-K−1) Diameter (m) Exponential integral function Heat transfer coefficient (W m−2 K−1) or thermal response factor. Pipe segment length (m) Thermal conductivity (W m−1 K−1) Infinite line source function Pipe length per trench (m) Mass flow rate (kg-s−1) Heat load (W) Heat flux per unit length (W m−1) Heat flux per unit area (W m−2) Radial coordinate (m) Unit length thermal resistance (K m W−1) Group of variables defined in Eq. (A6) Group of variables defined in Eq. (A6) Time (s) Temperature (K) Fluid velocity (m s−1)
Group of variables defined in Eq. (29) Group of variables defined in Eq. (31) Cartesian coordinates system (m)
Greek letters Thermal diffusivity (m2 s−1) Reduced thermal conductance (Eq. (19)) Group of variables (defined in Eq. (11)) (K m W−1) Group of variables (defined in Eq. (19)) Time integral variable Dimensionless temperature (defined in Eq. (17)) Exit dimensionless temperature (defined in Eq. (20)) Annual frequency of temperature variations Group of variables (defined in Eq. (61))
α κ γ η τ θ θo ω ζ Subscripts p f o s
recent studies relying on numerical simulations include the work of Kupiec et al. [7], Kayaci and Demir [8] Congedo et al. [9], Gan [10] and Soni et al. [11]. Various configurations of the horizontal heat exchanger, including inline piping, slinky or spiral configurations have also recently been analysed [9,12–14]. Some experimental works on horizontal configurations are also available [13,15]. In the case of vertical boreholes, significant recent research effort has been devoted to the analytical method [16–18]. A recent study by Fontaine et al. [19] follows similar steps as in the case of horizontal loops. In their approach, they split a horizontal loop into several discrete sections, and then extend the Claesson and Dunand [4] method developed for the steady-state analysis case to the unsteady case. To that end, they use the finite line source (FLS) solution to evaluate the thermal response of the discrete pipe sections in the case of the time-varying heat exchange, but assuming that at each time step, the heat transfer rate variation along the pipe follows the steady-state Claesson and Dunand solution. In this paper, we follow similar steps, but the heat transfer rate distribution along the pipe is not known a priori; rather, it is part of the final solution. Another simplification is also used for the FLS implementation. The model addresses some, but not all, of the horizontal loop distinct characteristics taken into account. The impact of the ground boundary and the seasonal temperature variation are however covered in a rigorous manner. The soil is assumed to be homogeneous, and for now, no
At the pipe radius Fluid Far-field value Soil values
mass transfer is assumed. The possible phase change associated with freezing is not taken into account. It must be pointed out that in the Fontaine model, the freezing of the permafrost is discussed, but not solved explicitly. If no phase change occurs and the ground is assumed homogeneous, the temperature variation in the ground can easily be modelled if one assumes a sinusoidal variation of the temperature at the air-soil boundary. A well-known solution in this regard was given by Kusuda and Achenback [20]. Assuming that the surface temperature is known and given by:
Ts (t ) = To - A cos (
(t
(1)
tshift ))
The solution in the ground is:
Tground (z, t ) = To - A exp - z
365
cos
(t
tshift -
z 2
365
) (2)
where To is the mean temperature far from the boundary, ω is the annual frequency, A is the variation amplitude of the temperature at the surface and tshift is the coldest day of the year. In Kusuda’s relation, the surface temperature at z = 0 is assumed to be known and follow a periodical variation with an annual frequency. Other models exist in which the air temperature is assumed known instead of the surface
Fig. 1. Horizontal heat exchanger. 535
Applied Thermal Engineering 155 (2019) 534–545
L. Lamarche
temperature [21]. In this case, it was shown by Badache et al. [21] that the mean soil temperature To can be found from the air mean temperature by the following expression:
1 To = [hr Tamb + hrad Tsky + G he
c hconv f (1
rh )]
1T (t ) p2
1T (t ) p
Rs =
hconv :mean air convection coefficient hevap :mean evaporation coefficient
c, f
= hconv + rh hevap
1 2 ks
2. Ground models
where
Fig. 1 gives an example of a simple piping in the ground exchanging heat with a heat transfer fluid. Assuming a homogeneous ground, the solution is given by the heat conduction equation: 2T
R1 = =
(4)
,
k
r = rp
=
= qp (t )
1T (t ,
2T (t ,
x , y, z ) +
rp q I 2 ks 2 t
z2 2
+I
z1 t
I
( ) zI ( rp
2
t
z2 2
= z1 t
1T (t ) + 1T (t ) p1 p2
2
) I(
z2 + z1 2 t
)
= q Rs 0.5 I
( ) + I( ) z1 t
(5)
z)
I
+ q Rp = q (Rs + Rp ) where Rp =
log(d o/di ) 1 + 2 kp h di
1 2z1 1 2z 2 log + Rp , R2 = log + Rp , R12 2 ks rpo 2 ks rpo 1 z + z2 log 1 2 ks z2 z1
(10)
L mCp
(11)
2z1 2
t
I
R2 R12 R1 + R2
/2 2R12
q2 = 2q
R1 R12 + /2 R1 + R2 2R12
tc, v =
L2 9
,
X=
r 2
t
(12)
(13a)
whereas for the horizontal configuration, it is in the order of:
z2 + z1 2 t (6a)
1 E1 (X 2) 2
(8)
(9)
where the last two terms represent the contribution of the two images. q′ is the heat transfer rate per unit length of pipe, not the trench. Under the assumption of uniform heat flux, the response factor I(X) is known to be the classical Infinite Line Source (ILS) solution [25]:
I (X ) =
z2 t
The justification of using a steady-state analysis for horizontal ground heat exchangers is sometimes justified knowing that the time scales in horizontal configurations are much shorter than for vertical boreholes. In this latter case, Eskilson [26] showed that the time scale for the steady-state regime was:
Using the Kusuda assumption, T is simply given by Eq. (2). As stated in the Introduction, this temperature field can be evaluated differently, or even measured, the premise here is that it is known. The main problem is to find the other solution. A well-known solution is based on the assumption that the heat flux is the same for all pipes and uniform along the pipe length. Using this assumption, early ground response factors [24] were found using the theory of images, forcing the homogeneous boundary condition to be fulfilled. The solution at the periphery of the first pipe for the case shown above would then be:
=
z2 + z1 2 t
2 2R1 R2 2R12 + 2/2 = q Req2 R1 + R2 2R12
q1 = 2q
2
1T (t ) p1
I
And
Due to the linearity of the heat equation, the solution is usually found from the superposition principle. The last problem can be seen as the superposition of the following problem (see Fig. 2): is replaced by the following equivalent problem (see Fig. 3): So that the final solution is:
T (t , x , y , z ) =
2z2 2 t
The simplified assumption of a uniform heat flux along the pipes was discarded by Claesson and Dunand [4], at least in the steady-state case for some simple configurations. In their analysis, they assumed a linear variation of the temperature along the pipes and found that for the case shown in Fig. 1:
Tf = q
T (x , y , 0) = Ts (t )
I
Req1 1T p
Tf =
More details on the numerical values of the constants can be found in [21] and [22]. In the case where the meteorological data does not follow a simple sinusoidal pattern, it is possible to decompose the signal into several Fourier terms where the previous description would be the main mode [23].
T n
z1 t
Several “soil resistances” are given in graphical form in [24] for different trench configurations based on this principle. The mean fluid temperature is found from the pipe resistance:
:mean radiative coefficient = hconv + hevap + hrad
rh :mean air relative humidity :empirical constant for the evaporation flux
1 T = t
z2 2
(7)
:mean solar radiation G :mean ground solar absortion coefficient
g
+ I
The response factor of the two-pipe arrangement can then be approximated as
(3)
Tamb :mean air temperature Tsky :mean sky temperature
hr
rp q I 2 ks 2 t
(6c)
where:
hrad he
=
(6b)
And E1 is the first-order exponential integral. The solution at the second pipe would be:
Fig. 2. Horizontal heat exchanger cross-section. 536
Applied Thermal Engineering 155 (2019) 534–545
L. Lamarche
Fig. 3. Horizontal heat exchanger equivalent problem.
tc, h =
10 z o2
The solution to Eq. (14) is given based on the dimensionless temperature:
(13b)
This can be seen in Fig. 4, where the steady-state thermal resistance for a single pipe is compared with the unsteady one derived from the ILS approach. Fontaine et al. [19] extended the approach of Claesson and Dunand to the unsteady case. In their approach, they could not solve the whole problem as a global unsteady thermal response factor. Instead, they split the horizontal piping into segments. Mutual thermal responses were evaluated with the finite line source (FLS) model. In their solution, the total heat exchange through the heat exchanger varied with time, but they assumed that, at a given time, the heat transfer distribution along the pipe followed the steady distribution found with the Claesson and Dunand distribution. In this paper, an approach similar to the one followed by Fontaine et al. is used, but no assumption is made regarding the heat distribution along the pipe. Instead, this distribution is part of the solution. The model is also based on the FLS model, but uses modified mathematical expressions of the FLS that were initially solved in the context of vertical boreholes and applied here in horizontal configurations [16,17,27]. The model also takes into account seasonal temperature variations that have a significant influence on the thermal behavior of horizontal geothermal fields. The model assumes, however, that the ground is homogeneous and single-phase.
T1 Tin
dT2 dx
= q2 =
+
(To T2) R2
(T1 T2) R12
To) = q1 R1 + q2 R12
(T2
To) = q1 R12 + q2 R2
x = m
(1 )
x ))
cosh( (1 x )) cosh( ) +
2)
m sinh (
(1 sinh ( ) m
x ))
x L
=
,
1
=
1+ 2
2
R1
,
,
2
=
R2 2 m
=
,
12
+2
=
R12
(19)
12 m
It is easy to verify that in the case of 1 = 2 , the classical Zeng’s profile [29] is obtained. The exit temperature is then: o
=
T2 (0) To = Tin To
cosh( ) cosh( ) +
m sinh ( m sinh (
) )
(20)
From which an effective resistance can be found:
R eq3 =
R1 + R2 2
(1 + o) = (1 o)
cotanh( )
(21)
This is the classical result found in a U-tube borehole if both resistances are equal, but with the difference that now the heat rate is per 0.6
0.5
(14)
0.4
0.3
0.2 Unsteady Steady
0.1
(15)
Inverting the last system will give: 2 2 2 R1 R2 R12 R1 R2 R12 R1 R2 R12 R1 = , R2 = , R12 = R2 R12 R1 R12 R12
m sinh (
(18b)
This system is similar to a typical U-tube model inside a vertical heat exchanger (see, for example [28]), except that the undisturbed ground temperature replaces the borehole temperature, and that the resistance associated with each pipe is now different. The Δ formalism is used as a reference in this case. Those resistances can easily be found from the values given by Eq. (10). In fact, we have:
(T1
m sinh (
with:
(T2 T1) R12
+
cosh( (1 x )) + cosh( ) +
2)
To x( 1 = exp To 2
Soil Resistance
mCp
(To T1) R1
(17)
To x( 1 = exp To 2
T2 Tin
Before looking at the transient case, it is interesting to revisit the steady case. The simple two-pipe case (Fig. 1) will be studied here. The Claesson and Dunand solution was found, assuming a linear variation of the temperature along the pipe. It is known that this is often a valid assumption, but for long pipes with low flow rate, it may lead to errors. The exact solution can be easily found in some simple configurations solving the following differential equations system: dT
To To
(18a)
3. New model for steady-state simulations
mCp dx1 = q1 =
T Tin
=
-4
-3
-2
-1
0
1
2
log(t/tc)
Fig. 4. Comparison between the unsteady and steady-state soil resistances for a single pipe for zo/rp = 50.
(16) 537
Applied Thermal Engineering 155 (2019) 534–545
L. Lamarche
unit length of pipe, and not per unit length of borehole.
Tb (t ) =
1 do / di + log di hconv 2 kpipe
T1
2
=
q q hFLS = ( (r , t , H1, H2, d2 2 ks 2 ks
(r , t , H1, H2, × [ierf(( z + H2) s )
(22)
d1) (27a)
22 1 e r s ds z) = 2H2 1/ 4 t s2
ierf( z s ) + ierf(( z
H1) s )
ierf(( z + H2
H1) s )]
T1
2, hor
=
q q hFLS, hor ,1,2 = ( (r1, t , H1, H2, x ) 2 k 2 k (28)
(r2, t , H1, H2, x )) where
r1 = rp
The method of images described in the previous section is applicable for the unsteady case, where it is again assumed that the heat flux is uniform along the pipe. In order to take into account the variation of the heat flux, Fontaine et al. suggested splitting the pipe into segments, with the heat transfer rate varying from segment to segment. At each time step, the total heat transfer is assumed to be known, and the variation along the pipe is assumed to vary exponentially. From the known heat transfer distribution along the pipe, the fluid temperature variation is calculated with the finite line source model (FLS). In its most general form, the FLS model is given [30] by (see Fig. 8):
,segment on the same pipe section
dA2
q
t
dA1
0
(
(t
d2 4 (t )
e
))3/2
di2 4 (t )
e
r1 = (y1
y2 ) 2 + (d1
r2 = (y1
,segment on the same pipe section y2 ) 2 + (d1 + d2 )2 ,otherwise
r2 = 2di
d2 ) 2
,otherwise
4.2. Unsteady temperature history Knowing the response factor associated with a horizontal piping arrangement, it is possible to calculate the time evolution of the average pipe temperature over a pipe segment. This calculation assumes knowledge of the heat flux associated with every segment. Again, in Fontaine’s model, the total heat transfer time history is known and the variation along the pipe is assumed to follow an exponential variation. Since the heat transfer varies with time, a convolution scheme has to be performed to estimate the temperature evolution. This convolution can
d
(23) A common simplification is often done by integration over time. Doing so will result in:
T1
2 (t )
=
q 4 kH2
dA2
dA1
erfc(d/2 d
t)
erfc(di /2 di
t) (24)
which is the form used by Fontaine et al. in their model. This model involves a double integral, which has been simplified in the case of vertical boreholes by Lamarche and Beauchamp [31] and by Claesson and Javed [27]. The latter solution is of particular interest. The idea is to keep the time integral and invert the order of integration. In the case of a vertical borehole, the solution becomes (see Fig. 9): Tb (t ) =
d+H
8k H
d
dz
H d
d
q
t 0
(
(t
))3/2
e
r 2+ (z 4 (t
)2 )
e
r 2 + (z + )2 4 (t )
d
(25) with
r2
=
x2
+
(27b)
The generalization of the segment-to-segment temperature response in the general inclined case was studied by Lazarotto [33], where the author simplified the double integral involved, but could not reduce it into a single integral in the general case. Although it may not seem possible to simplify the thermal response in the general case, it can be done in the horizontal configuration. Looking at Fig. 10, it can be shown that the generalization of Cimmino’s expression is given by (Appendix A):
4.1. Finite line source model for horizontal pipes
8 k H2
ds
x2
(r , t , H1, H2, H1 + d1 + d2))
4. Unsteady analysis
=
ierf(2ds )
with ierf(x ) = x ·erf(x ) (1 e ) Details of the derivation can be found in [17]. The generalization of the previous expression for the influence of a vertical segment on another vertical segment was given by Cimmino and Bernier [32]
The results obtained by the finite element simulation are shown in Fig. 6, where they are compared with the analytic variation given by Eq. (18) and the linear variation, assuming a soil resistance given by Claesson and Dunand (Eq. (9)) and the classical method of images of Bose et al. (Eq. (8)). In the latter case, the thermal response is unsteady, and only the steady-state solution is shown. The corresponding resistances are given in Table 2, and the corresponding inlet and outlet temperatures in Table 3. A grid sensitivity was performed and is shown in Fig. 7 where it is seen that the effective resistance converge towards the value given in Table 2.
2 ( t)
2 2 2 ierf(Hs ) + 2 ierf(Hs + 2ds ) ierf(2Hs + 2ds ) e r s s2
1
A simulation using the finite element software COMSOL Multiphysics™ was done to compare the different steady-state solutions for the case just described (see Fig. 5). The geometry is the same as the one in Fig. 1. The domain was split in half to take into account the symmetry of the problem. The parameters of the simulation are given in Table 1. No thickness was simulated on the pipe. Instead the convection coefficient represent and equivalent transfer coefficient:
T1
1/ 4 t
(26)
3.1. Simulation results
1 = do h
q 4 kH
Fig. 5. Mesh generated for horizontal ground heat exchanger.
y2 538
Applied Thermal Engineering 155 (2019) 534–545
L. Lamarche
Table 1 GHE configuration. rb(m)
L(m)
u(m/s)
z1(m)
z2(m)
ρCp(MJ/ m3-K)
q(W)
ks(W/mK)
h(W/ m2K)
0.02
100
0.3
2
3
3.93
1000
8
1000
(ξi,ηi,ζi)
ground surface di
(ξ,η,ζ) A1
(x,y,z) d
A2
Fig. 8. Segment-to-segment temperature thermal influence.
product in frequency domain. In order to use the FFT approach, the heat transfer history has to be known a priori, which is not a trivial factor, especially in the case of horizontal piping, where the seasonal temperature variation can interfere with the process. The proposed model follows the path already proposed by the author in the case of mixed vertical boreholes [16], where any arrangements of vertical borehole in series and in parallel can be simulated. The idea is to split a horizontal pipe into several sections, each of which can be seen as a pipe in series with the previous section. Referring again to the typical configuration shown in Fig. 1, splitting the pipe into ns segments can be analysed by the tools given, assuming ns pipes in series. Two major differences have to be confronted in the case of horizontal systems: the response factors have to be changed and the seasonal temperature variation has to be taken into account. The first problem has been covered in the previous section, and the second one will now be discussed. The description will be done with only one parallel path, but the general case is given in [16]. It can be shown that the pipe temperature at a given time can be evaluated with the one at a previous time step in a form given by:
Fig. 6. Fluid temperature along the heat exchanger. Table 2 Equivalent resistance (m-K/W). Req1 (Eq. (8))
Req2 (Eq. (9))
Req3 (Eq. (21))
Req, comsol
0.145
0.163
0.167
0.169
Table 3 Inlet and Outlet temperatures. Req1 Tfin
Tfout
Req2 Tfin
Tfout
Req3 Tfin
Tfout
Req, comsol Tfin Tfout
−1.08
−0.40
−1.15
−0.4
−1.17
−0.50
−1.18
−0.51
ns
Tp, i = Sp, i +
Sq, ij Xj ( Tfi, j
Tp, j )
j=1
(29)
where the expressions of Sp, i, Sq, ij are given in [16] and reproduced in Appendix B. We also have:
qi = Xi (Tfi, i X=
(mCp )i Hi
where o, i
=
o, i
Tp, i )
(1
o, i )
(30)
represents the temperature variation along the segment:
Tfo, i
Tp, i
Tfi, i
Tp, i
(31)
Assuming a linear variation, we have: o, i
=
1 Yi 1 + Yi
,
Yi =
Hi 2(mCp)i Rp
(32)
Assuming an exponential variation, the result is: o, i
Fig. 7. Impact of grid refinement.
= exp( 2Yi )
(33)
Again, these assumptions are for a single segment, not the entire pipe. In the case of np parallel pipes, where all the inlet temperatures are known, Eq. (29) represents a system of np equations for the np unknowns. For the case of ns series arrangement, only the first temperature is known and the (na = ns−1) others are unknowns. The new equations that have to be added to the system are:
be very time-consuming, and Fontaine uses the Fast Fourier Transform (FFT) approach proposed by Marcotte and Pasquier [34] to accelerate the convolution evaluation. The FFT is a very powerful approach that allows the time domain convolution product to be replaced by a simple 539
Applied Thermal Engineering 155 (2019) 534–545
L. Lamarche
Fig. 9. Vertical segment-to-segment thermal response factor.
Tfi, i = Tfo, i
1
=
o, i 1 Tfi, i 1
+ Tp, i
1 (1
o, i 1)
(1 (34)
A ×T =B A =
LL =
A (ns × ns ) UR (ns × na ) LL (na × ns ) LR (na × na)
, B =
BU BL
, T =
Tp TL
1 LR =
(35)
A=
1 + X1 Sq,11 , X2 Sq,12, ... X1 Sq,21 , 1 + X2 Sq,22, .... . . ..
X2 Sq,12, X2 Sq,22,
o,2
0
0 1 o,3
0 (1
0 . .. o,2 )
(39)
0 0 . .. 0 0. .. 1 0. .. (40)
Sp,1 + X1 Sq,11 Tfin,1 BU = Sp,2 + X1 Sq,21 Tfin,1
(36)
(41)
Tfi,2 TL = T fi,3
UR =
o,1)
0
(37)
BL =
X3 Sq,13, ... X3 Sq,23, ...
Tfi,1 0
(42)
Once the inlet temperatures are known, the heat transfer distribution is evaluated with Eq. (30), the Sp coefficients are updated (Eqs.
(38)
Fig. 10. Horizontal segment-to-segment thermal response factor. 540
Applied Thermal Engineering 155 (2019) 534–545
L. Lamarche
on a pipe segment becomes:
Tf , i
qi =
Tp, i
=
Rp
( 1Tp, i +
Tf , i
2T ) p, i
=
Rp
1T p, i
Tf , i
(54)
Rp
2T where T = T and 1T are given by Kusuda equation (Eq. (2)) or any other model used for the seasonal temperature variations. System of Eqs. (29) becomes:
(55)
1T ) p,i
qi = Xi (Tfi,i
ns 1T p, i
= Sp, i +
1T ) p, j
Sq, ij Xj (Tfi, j
(56)
j=1
which means that the results given in the previous section will give the pipe temperature contribution 1T and the reduced fluid temperature. Particular attention is paid to the complementary system:
Tfi, i = Tfi, i
Fig. 11. Fluid temperature during one year (test1).
Tfi, i
o, ns
=
Tfi,1
A+
Tfo, ns
Tp, ns
Tfi, ns
Tp, ns
o, ns
×T +
A+ =
Tfi, ns
Tfi,1 (44)
Tp, ns (1
q o, ns ) = mCp
,
B+ =
BU BL
,
T+ =
UR =
0
0
0 . .. (1
0
1
0 o,1
LR =
0
0 1 o,2
(1
o, ns )
o,1)
(50)
0 . .. o, ns 0 0. .. 1 0. .. (51)
Sp,1 BU = Sp,2 (52)
BL =
q/ mCp 0
2T p, i 1)
2T p, i
2T ) p, i
o, ns
1T p, ns
Tfi, ns
(1
o, ns )
=
q + ( 2Tp, ns mCp
2T ) p,1
In the first simulation, the inlet temperature will be imposed at the inlet of the upward pipe section shown in Fig. 1. A linear profile as shown In Fig. 11 will be imposed during half of the year, and a constant temperature is imposed for the rest of the year. In the numerical simulation, a time-varying temperature profile given by Eq. (1) is imposed at the surface of the domain (z = 0), with A = 10 and tshift = 0 (January 1st is the coldest day of the year). A 5-hour time step was chosen for the numerical simulation, while a 1- hour-time step is taken during the analytic simulation. Eight segments (4 on each pipe section) were chosen for the analytic simulation. Fig. 11 shows the outlet temperature variation during the year, and Fig. 12, the total heat exchange with the ground. Fig. 13 shows the temperature variation and Fig. 14 the heat transfer per unit length along the pipe for two different simulation times: after half a year and after one year. The impact of the number of segments was studied for this particular test. The number of segments in each pipe section was varied and the difference between the total heat exchange in the ground was compared. Thus difference is defined as:
X2 Sq,12, ... X2 Sq,22, ... (49)
LL =
=(
+
2T p, i 1
6.1. Inlet temperature imposed
(48)
X1 Sq,11, X1 Sq,21,
o , i 1)
1T p, i 1
The model just described was compared with a finite element simulation on COMSOL™ Multiphysics. The geometry studied is the same as the one shown in Fig. 1, and the grid is the same as in the steady-state analysis.
Tp TL (47)
Tfi,1 TL = T fi,2
(1
o, i 1)(
o, i 1)
1T p, i 1
6. Simulation results
(46)
A (ns × ns ) UR (ns × ns ) LL (ns × ns ) LR (ns × ns )
o, i 1 Tfi, i 1
+ (1
1 (1
where the last term is zero again if the inlet and outlet are at the same depth.
(45)
=B+
o, i 1 Tfi, i 1
+ Tp, i
The last term in Eq. (57) is zero if the two segments are at the same depth, but different from zero otherwise. Eq. (45) should also be modified:
(43)
Tfo, ns )
=
o, i 1 Tfi, i 1
(57)
(A6)–(A7)), and the scheme moves to the next time step. In the case where the inlet temperature is not known, but the total heat transfer is imposed, the system is modified:
q = (mCp)(Tfi,1
2T 1p, i
(53)
Qanalytical =
qanalytical (t ) dt
Qcomsol =
qcomsol (t ) dt |Qanalytical Qcomsol |
5. Impact of seasonal temperature variation
difference (%) =
The effect of seasonal variation is quite simple when the heat flux is known, and less trivial when it is part of the solution. The heat transfer
between the analytical model and the numerical simulation was evaluated and the result is shown in Fig. 15. 541
Qcomsol
× 100
(58)
Applied Thermal Engineering 155 (2019) 534–545
L. Lamarche
Fig. 12. Total heat transfer during one year (test1).
Fig. 15. Influence of the number of segments on the difference between analytical and numerical solution.
6.2. Total heat transfer imposed As a second test, the total heat flux was imposed to test the systems of Eqs. (46)–(53). A constant heat extraction of 5000 W was imposed for the first part of the year (t < 4380 h), while no heat exchange (restitution period) was imposed for the rest of the year. All other parameters were the same as the one given in Table 4. The inlet and outlet temperatures are shown in Fig. 16. The temperature and heat flux along the pipe are now given just before the end of the heat injection (t = 4375 h), and 100 h after the power is turned off (t = 5380 h), and shown in Fig. 17. The heat transfer per unit length is shown in Fig. 18. It can be seen that when the total heat exchange is zero, the local heat transfer is not, due in part to the temperature difference between the ground temperature around the pipe. Finally, the calculation times for both approaches are compared in Table 5, where an important reduction in the calculation time is observed using the proposed model.
Fig. 13. Temperature variation along the pipe for two different times (test1).
6.3. Comparison with simpler models As stated in Section 2, the ILS model was often used to generate unsteady soil resistances for horizontal piping. In the case where the fluid temperature is imposed, it is not trivial to apply the model but in the case where the heat transfer is imposed, superposition of the solution given by the method of images and the seasonal variation is possible. For example the second test analyzed in the last subsection can be compared with the method of images. From Eq. (8): 1T f
= q Rs + Rp
With Rs given by Eq. (7) for the two-pipe arrangement. The inletoutlet temperatures are then 1T f ,o
Tf , o =
=
1T f
1T f ,o
+ +
q 2mCp
,
2T , p, o
Tf , i =
1T f ,o
=
1T f ,i
1T f
q 2mCp
+ , 1Tp, i
(59a) (59b)
Comparison of the inlet-outlet temperatures thus obtained is compared with the numerical simulation in Fig. 19. A zoom in the vicinity
Fig. 14. Heat transfer per unit length variation along the pipe for two different times (test1).
542
Applied Thermal Engineering 155 (2019) 534–545
L. Lamarche
Table 4 GHE configuration. rb(m)
L(m)
u(m/s)
z1(m)
z2(m)
(ρCp)f(MJ/m3-K)
ks(W/mK)
(ρCp)s(MJ/m3-K)
h(W/m2K)
ns
0.02
100
0.3
2
3
3.93
8
2.0
100
8
Fig. 16. Fluid temperature during one year (test2).
Fig. 18. Heat transfer per unit length variation along the pipe for two different times (test2). Table 5 CPU time. Test1 (Comsol)
Test1 (Analytic)
Test2 (Comsol)
Test2 (Analytic)
18 min
2.2 sec
21 min
2.2 sec
Fig. 17. Temperature variation along the pipe for two different times (test2).
of the end of the heat pulse is done to better see the difference in the models. Although the general trend of the temperature evolution is correctly predicted with the ILS model, larger variations are observed near the regions where the heat flux varies. Most of the differences come from the assumption that the heat transfer per unit length is uniform. The rigorous approach was described in the previous sections. However a simple variation of the method of images using the ILS can be done if we allow the heat transfer to vary from one pipe to the other. A possibility is to assume that, even in the unsteady case, the heat transfer distribution resembles the steady case given by the Claesson and Dunand [4] relation (Eq. (12)). Rewriting those equations
R R12 /2 q1 = 2q 2 = R1 + R2 2R12
1q
R R12 + /2 q2 = 2q 1 = R1 + R2 2R12
Fig. 19a. Outlet temperature given by Eq. (59)
Rs =
+
(61)
2 R2
where R1 and R2 are given by Eq. (6). Replacing this new resistance in Eq. (59) 1T f ,o
= q Rs + Rp +
2q
Tf , o =
(60) as:
1 R1
1T f ,o
+
2T p, o
q 2mCp
,
,
Tf , i =
1T f ,o
= q Rs + Rp
1T f ,i
+ , 2Tp, i
q 2mCp
(62a) (62b)
would give the results shown in Fig. 20: It is seen that even though, the unsteadiness is still not well
The unsteady soil resistance given by Eq. (7) can now be rewritten 543
Applied Thermal Engineering 155 (2019) 534–545
L. Lamarche
Fig. 19b. Inlet temperature given by Eq. (59).
Fig. 20b. Inlet temperature given by Eq. (62).
source associated to horizontal configurations was developed to simulate the heat transfer between a horizontal heat exchanger and the surrounding ground. The model was compared to a finite element simulation in the case of a simple configuration and display excellent agreement given that the model is 500 0times faster than the numerical simulation. While the configuration may be simple, it illustrates a very important aspect of horizontal systems, namely, different local ground temperatures around pipes at different heights and how this can affect the thermal behavior of the ground exchanger. The model can easily be extended to different inline configurations, which can have parallel branches as well. Extensions to slinky or spiral configurations can also be considered, but in that case, the thermal response factor between pipe sections would be more complex. An extension of the classical work of Claesson and Dunand was also presented as part of this study. It is a future goal to use it to provide potentially better guidelines for horizontal design procedures. Acknowledgements
Fig. 20a. Outlet temperature given by Eq. (62).
We would like to thank the Natural Sciences and Engineering Research Council of Canada (NSERC), which partly financed this research (RGPIN-2014-06240).
predicted, the overall comparison is better. 7. Conclusion An analytical model based on a new formulation of the finite line Appendix A Referring to Fig. we have:
q Ti j (t ) = 8k H
x + H2 x
dx
H1 0
1
t
d
0
(
))3/2
(t
r 2 + (x 4 (t
e
)2 )
e
2 rim + (x + )2 4 (t )
d (A1)
where r represents the perpendicular distance between the source line and the surface of interest. In the case where two segments are on the same pipe section, we then have:
r = rp
,
rim = 2z i
In the case where the two segments are in different pipe sections:
r2
= (yi
yj ) 2 + (z i
zj ) 2
,
2 rim = (yi
yj ) 2 + (z i + zj ) 2
Using the usual change of variable s = 1/ 4 (t
Ti
j (t )
=
q 2 4 k H2
1/ 4 t
[e
r 2s 2
e
ri2 s 2] d
s
x + H2 x
(A2)
) leads to: dx
H1 0
e
s 2 (x
)2 d
544
(A3)
Applied Thermal Engineering 155 (2019) 534–545
L. Lamarche
Ti
j (t )
=
q 4 k H2
Ti
j (t )
=
q 4 k H2
1/ 4 t
r 2s 2
[e
2 2 e ri s
22 e r s
s2
1/ 4 t
ierf(( x
e
ri2 s 2]
ds s
x + H2 x
(erf(sx )
(ierf(( x + H2) s )
H1) s ) ierf(( x + H2
erf(s (x
H1)) dx
(A4)
ierf( x s )+ (A5)
H1) s ) ds
Appendix B
t =
Sp, i =
t r p2 1 ks
Fi ( t +
Nz
e
zn2 t F ( t i
, zn) z n
, Sq, ij =
n=1
t , z n) = e
zn2 t F
i(t
, z n) +
nserie j=1
1 ks
Nz
(1
e
zn2 t
) uij (z n) zn
(A6)
n= 1
qj ( t +
t )(1
e
zn2 t
) uij (z n) (A7)
Fi ( t = 0, zn ) = 0 uij (z ) =
z
L 1 (hFLS, hor , ij )
(A8)
Details on the derivation of these expressions are given in [16]. The choice of the integral variable z is discussed in [35]. Appendix C. Supplementary material Supplementary data to this article can be found online at https://doi.org/10.1016/j.applthermaleng.2019.04.006.
ground heat exchangers, Int. J. Heat Mass Transf. 115 (2017) 354–360. [19] P.-O. Fontaine, D. Marcotte, P. Pasquier, D. Thibodeau, Modeling of horizontal geoexchange systems for building heating and permafrost stabilization, Geothermics 40 (3) (2011) 211–220. [20] T. Kusuda, P.R. Achenbach, Earth temperature and thermal diffusivity at selected stations in the United States, National Bureau of Standards Gaithersburg MD (1965). [21] M. Badache, P. Eslami-Nejad, M. Ouzzane, Z. Aidoun, L. Lamarche, A new modeling approach for improved ground temperature profile determination, Renew. Energy 85 (2016) 436–444. [22] G. Mihalakakou, On estimating soil surface temperature profiles, Energy Build. 34 (3) (2002) 251–259. [23] C. Jacovides, G. Mihalakakou, M. Santamouris, J. Lewis, On the ground temperature profile for passive cooling applications in buildings, Sol. Energy 57 (3) (1996) 167–175. [24] J.E. Bose, J.D. Parker, F.C. McQuiston, American Society of Heating R, Engineers AC. Design/data Manual for Closed-loop Ground-coupled Heat Pump Systems: American Society of Heating, Refrigerating, and Air-Conditioning Engineers, 1985. [25] H.S. Carslaw, J.C. Jaeger, Conduction of Heat in Solids, 2nd ed., 1959. [26] P. Eskilson, Thermal analysis of heat extraction boreholes, 1987. [27] J. Claesson, S. Javed, An analytical method to calculate borehole fluid temperatures for time-scales from minutes to decades, in: Conference An Analytical Method to Calculate Borehole Fluid Temperatures for Time-scales from Minutes to Decades, vol. 117, pp. 279–288. [28] L. Lamarche, S. Kajl, B. Beauchamp, A review of methods to evaluate borehole thermal resistances in geothermal heat-pump systems, Geothermics 39 (2) (2010) 187–200. [29] H. Zeng, N. Diao, Z. Fang, Heat transfer analysis of boreholes in vertical ground heat exchangers, Int. J. Heat Mass Transf. 46 (23) (2003) 4467–4481. [30] L. Lamarche, Analytical g-function for inclined boreholes in ground-source heat pump systems, Geothermics 40 (4) (2011) 241–249. [31] L. Lamarche, B. Beauchamp, A new contribution to the finite line-source model for geothermal boreholes, Energy Build. 39 (2) (2007) 188–198. [32] M. Cimmino, Fast calculation of the g-functions of geothermal borehole fields using similarities in the evaluation of the finite line source solution, J. Build. Perform. Simul. 1–14 (2018). [33] A. Lazzarotto, A methodology for the calculation of response functions for geothermal fields with arbitrarily oriented boreholes b, Renew. Energy 86 (2016) 1380–1393. [34] D. Marcotte, P. Pasquier, Fast fluid and ground temperature computation for geothermal ground-loop heat exchanger systems, Geothermics 37 (6) (2008) 651–665. [35] L. Lamarche, B. Beauchamp, A fast algorithm for the simulation of GCHP systems, ASHRAE Trans. 113 (1) (2007) 470–476.
References [1] J.W. Lund, T.L. Boyd, Direct utilization of geothermal energy 2015 worldwide review, Geothermics 60 (2016) 66–93. [2] J. Raymond, M. Malo, D. Tanguay, S. Grasby, F. Bakhteyar, Direct utilization of geothermal energy from coast to coast: a review of current applications and research in Canada, in: Conference Direct Utilization of Geothermal Energy from Coast to Coast: A Review of Current Applications and Research in Canada. [3] S. Rees, Advances in Ground-source Heat Pump Systems, Woodhead Publishing is an imprint of Elsevier, Duxford, UK, 2016. [4] J. Claesson, A. Dunand, Heat Extraction from Ground by Horizontal Pipes, Swedish Council of building Research, Stokholm, 1983. [5] V.C. Mei, Horizontal ground-coil heat exchanger theoretical and experimental analysis, United States: Report No. ORNL/CON-193, Dec, Oak Ridge National Laboratory, Oak Ridge, 1986. [6] M. Piechowski, Heat and mass transfer model of a ground heat exchanger: theoretical development, Int. J. Energy. Res. 23 (7) (1999) 571–588. [7] K. Kupiec, B. Larwa, M. Gwadera, Heat transfer in horizontal ground heat exchangers, Appl. Therm. Eng. 75 (2015) 270–276. [8] N. Kayaci, H. Demir, Numerical modelling of transient soil temperature distribution for horizontal ground heat exchanger of ground source heat pump, Geothermics 73 (2018) 33–47. [9] P. Congedo, G. Colangelo, G. Starace, CFD simulations of horizontal ground heat exchangers: A comparison among different configurations, Appl. Therm. Eng. 33 (2012) 24–32. [10] G. Gan, Dynamic thermal modelling of horizontal ground-source heat pumps, Int. J. Low-Carbon Technol. 8 (2) (2013) 95–105. [11] S.K. Soni, M. Pandey, V.N. Bartaria, Ground coupled heat exchangers: A review and applications, Renew. Sustain. Energy Rev. 47 (2015) 83–92. [12] Z. Xiong, D.E. Fisher, J.D. Spitler, Development and validation of a Slinkyb “ ground heat exchanger model, Appl. Energy 141 (2015) 57–69. [13] S. Yoon, S.-R. Lee, G.-H. Go, Evaluation of thermal efficiency in different types of horizontal ground heat exchangers, Energy Build. 105 (2015) 100–105. [14] H. Fujii, K. Nishi, Y. Komaniwa, N. Chou, Numerical modeling of slinky-coil horizontal ground heat exchangers, Geothermics 41 (2012) 55–62. [15] H. Esen, M. Inalli, M. Esen, K. Pihtili, Energy and exergy analysis of a groundcoupled heat pump system with two horizontal ground heat exchangers, Build. Environ. 42 (10) (2007) 3606–3615. [16] L. Lamarche, Mixed arrangement of multiple input-output borehole systems, Appl. Therm. Eng. 124 (2017) 466–476. [17] M. Cimmino, M. Bernier, A semi-analytical method to generate g-functions for geothermal bore fields, Int. J. Heat Mass Transf. 70 (2014) 641–650. [18] L. Lamarche, G-function generation using a piecewise-linear profile applied to
545
Contents lists available at ScienceDirect
Applied Thermal Engineering journal homepage: www.elsevier.com/locate/apthermeng
Research Paper
Horizontal ground heat exchangers modelling Louis Lamarche
T
École de technologie supérieure, 1100 Notre-Dame Streetr West, Montreal H3C 1K3, Canada
HIGHLIGHTS
model for horizontal ground heat exchanger. • New model which is very fast and accurate with respect to numerical simulation. • Analytical comparison with the classical trench model is performed. • AParticularity of horizontal systems are discussed possible applications of the model are discussed. • ARTICLE INFO
ABSTRACT
Keywords: Ground source heat pump Horizontal ground heat exchanger Finite line source
Geothermal heat exchangers applied to building air conditioning and heating systems are often categorized as vertical and horizontal configurations. In this article, an analytical model is presented to simulate on an hourly basis, ground heat exchangers having a horizontal configuration. The model is based on a new formalism of the finite line source approach that allows faster simulation. The model presents a novel approach to take into the effect of seasonal temperature variations especially when the inlet temperatures in the ground heat exchanger are imposed. It also takes into account the effect of the air-soil boundary condition and the thermal interference between the different portions of pipes in the buried field. The temperature as well as the heat flux distribution along the ground heat exchanger is part of the solution and not a priori assumed as in the case of previous models. The model does not consider the possible phase change due to ice formation or the moisture content variation associated with some horizontal configurations. It is based on recent studies related to vertical boreholes and adapted to horizontal systems. It is intended to be a valuable tool to simulate on an hourly basis, the time response of the ground heat exchanger coupled with a building equipped with heat pumps.
1. Introduction Ground-Coupled Heat Pumps (GCHP) using a secondary loop ground heat exchanger are often categorized as vertical and horizontal loop systems, although hybrid systems such as inclined or pile systems, are also available. Vertical borehole systems are the most popular technology in Europe, but in North America, horizontal systems continue to be largely used [1,2]. This observation notwithstanding, there is general agreement that the scientific literature on horizontal systems is much more sparse than that on vertical ones [3]. Design approaches often continue to be associated with rules-of-thumb and approximate charts [3]. Even though some common physical aspects link both technologies, there are some specific factors associated with horizontal loops that make the design of such systems unique. The following is a non-exhaustive list of such factors: 1. The ground surface has a greater influence on the thermal response
than for vertical systems. 2. Seasonal temperature variations, which are often neglected in vertical systems, are very important in horizontal ones. 3. Horizontal piping lay in trenches in the subsurface overburden where the moisture content can influence the thermal performance of the heat exchanger. 4. Phase change due to freezing and thawing may occur. 5. Several configurations are possible, including inline, slinky, etc., which cause pipe-to-pipe interference. Given the above factors, horizontal loops are more complicated to analyse, even though they may be easier to install. Several studies have looked at these systems. A well-known report on horizontal systems [4] by Claesson and Dunand is often used as a reference even though several aspects shown in the list above not addressed. Mei [5] develop a numerical model for horizontal systems, which was extended by Piechowski [6] to include mass transfer due to moisture content. More
E-mail address: [email protected]. https://doi.org/10.1016/j.applthermaleng.2019.04.006 Received 5 December 2018; Received in revised form 20 February 2019; Accepted 2 April 2019 Available online 04 April 2019 1359-4311/ © 2019 Elsevier Ltd. All rights reserved.
Applied Thermal Engineering 155 (2019) 534–545
L. Lamarche
Nomenclature A Cp D E1 h H k I L m q q’ q” r R Sp Sq t T u
X Y x−y−z
Area (m2) Specific heat capacity (J-Kg−1-K−1) Diameter (m) Exponential integral function Heat transfer coefficient (W m−2 K−1) or thermal response factor. Pipe segment length (m) Thermal conductivity (W m−1 K−1) Infinite line source function Pipe length per trench (m) Mass flow rate (kg-s−1) Heat load (W) Heat flux per unit length (W m−1) Heat flux per unit area (W m−2) Radial coordinate (m) Unit length thermal resistance (K m W−1) Group of variables defined in Eq. (A6) Group of variables defined in Eq. (A6) Time (s) Temperature (K) Fluid velocity (m s−1)
Group of variables defined in Eq. (29) Group of variables defined in Eq. (31) Cartesian coordinates system (m)
Greek letters Thermal diffusivity (m2 s−1) Reduced thermal conductance (Eq. (19)) Group of variables (defined in Eq. (11)) (K m W−1) Group of variables (defined in Eq. (19)) Time integral variable Dimensionless temperature (defined in Eq. (17)) Exit dimensionless temperature (defined in Eq. (20)) Annual frequency of temperature variations Group of variables (defined in Eq. (61))
α κ γ η τ θ θo ω ζ Subscripts p f o s
recent studies relying on numerical simulations include the work of Kupiec et al. [7], Kayaci and Demir [8] Congedo et al. [9], Gan [10] and Soni et al. [11]. Various configurations of the horizontal heat exchanger, including inline piping, slinky or spiral configurations have also recently been analysed [9,12–14]. Some experimental works on horizontal configurations are also available [13,15]. In the case of vertical boreholes, significant recent research effort has been devoted to the analytical method [16–18]. A recent study by Fontaine et al. [19] follows similar steps as in the case of horizontal loops. In their approach, they split a horizontal loop into several discrete sections, and then extend the Claesson and Dunand [4] method developed for the steady-state analysis case to the unsteady case. To that end, they use the finite line source (FLS) solution to evaluate the thermal response of the discrete pipe sections in the case of the time-varying heat exchange, but assuming that at each time step, the heat transfer rate variation along the pipe follows the steady-state Claesson and Dunand solution. In this paper, we follow similar steps, but the heat transfer rate distribution along the pipe is not known a priori; rather, it is part of the final solution. Another simplification is also used for the FLS implementation. The model addresses some, but not all, of the horizontal loop distinct characteristics taken into account. The impact of the ground boundary and the seasonal temperature variation are however covered in a rigorous manner. The soil is assumed to be homogeneous, and for now, no
At the pipe radius Fluid Far-field value Soil values
mass transfer is assumed. The possible phase change associated with freezing is not taken into account. It must be pointed out that in the Fontaine model, the freezing of the permafrost is discussed, but not solved explicitly. If no phase change occurs and the ground is assumed homogeneous, the temperature variation in the ground can easily be modelled if one assumes a sinusoidal variation of the temperature at the air-soil boundary. A well-known solution in this regard was given by Kusuda and Achenback [20]. Assuming that the surface temperature is known and given by:
Ts (t ) = To - A cos (
(t
(1)
tshift ))
The solution in the ground is:
Tground (z, t ) = To - A exp - z
365
cos
(t
tshift -
z 2
365
) (2)
where To is the mean temperature far from the boundary, ω is the annual frequency, A is the variation amplitude of the temperature at the surface and tshift is the coldest day of the year. In Kusuda’s relation, the surface temperature at z = 0 is assumed to be known and follow a periodical variation with an annual frequency. Other models exist in which the air temperature is assumed known instead of the surface
Fig. 1. Horizontal heat exchanger. 535
Applied Thermal Engineering 155 (2019) 534–545
L. Lamarche
temperature [21]. In this case, it was shown by Badache et al. [21] that the mean soil temperature To can be found from the air mean temperature by the following expression:
1 To = [hr Tamb + hrad Tsky + G he
c hconv f (1
rh )]
1T (t ) p2
1T (t ) p
Rs =
hconv :mean air convection coefficient hevap :mean evaporation coefficient
c, f
= hconv + rh hevap
1 2 ks
2. Ground models
where
Fig. 1 gives an example of a simple piping in the ground exchanging heat with a heat transfer fluid. Assuming a homogeneous ground, the solution is given by the heat conduction equation: 2T
R1 = =
(4)
,
k
r = rp
=
= qp (t )
1T (t ,
2T (t ,
x , y, z ) +
rp q I 2 ks 2 t
z2 2
+I
z1 t
I
( ) zI ( rp
2
t
z2 2
= z1 t
1T (t ) + 1T (t ) p1 p2
2
) I(
z2 + z1 2 t
)
= q Rs 0.5 I
( ) + I( ) z1 t
(5)
z)
I
+ q Rp = q (Rs + Rp ) where Rp =
log(d o/di ) 1 + 2 kp h di
1 2z1 1 2z 2 log + Rp , R2 = log + Rp , R12 2 ks rpo 2 ks rpo 1 z + z2 log 1 2 ks z2 z1
(10)
L mCp
(11)
2z1 2
t
I
R2 R12 R1 + R2
/2 2R12
q2 = 2q
R1 R12 + /2 R1 + R2 2R12
tc, v =
L2 9
,
X=
r 2
t
(12)
(13a)
whereas for the horizontal configuration, it is in the order of:
z2 + z1 2 t (6a)
1 E1 (X 2) 2
(8)
(9)
where the last two terms represent the contribution of the two images. q′ is the heat transfer rate per unit length of pipe, not the trench. Under the assumption of uniform heat flux, the response factor I(X) is known to be the classical Infinite Line Source (ILS) solution [25]:
I (X ) =
z2 t
The justification of using a steady-state analysis for horizontal ground heat exchangers is sometimes justified knowing that the time scales in horizontal configurations are much shorter than for vertical boreholes. In this latter case, Eskilson [26] showed that the time scale for the steady-state regime was:
Using the Kusuda assumption, T is simply given by Eq. (2). As stated in the Introduction, this temperature field can be evaluated differently, or even measured, the premise here is that it is known. The main problem is to find the other solution. A well-known solution is based on the assumption that the heat flux is the same for all pipes and uniform along the pipe length. Using this assumption, early ground response factors [24] were found using the theory of images, forcing the homogeneous boundary condition to be fulfilled. The solution at the periphery of the first pipe for the case shown above would then be:
=
z2 + z1 2 t
2 2R1 R2 2R12 + 2/2 = q Req2 R1 + R2 2R12
q1 = 2q
2
1T (t ) p1
I
And
Due to the linearity of the heat equation, the solution is usually found from the superposition principle. The last problem can be seen as the superposition of the following problem (see Fig. 2): is replaced by the following equivalent problem (see Fig. 3): So that the final solution is:
T (t , x , y , z ) =
2z2 2 t
The simplified assumption of a uniform heat flux along the pipes was discarded by Claesson and Dunand [4], at least in the steady-state case for some simple configurations. In their analysis, they assumed a linear variation of the temperature along the pipes and found that for the case shown in Fig. 1:
Tf = q
T (x , y , 0) = Ts (t )
I
Req1 1T p
Tf =
More details on the numerical values of the constants can be found in [21] and [22]. In the case where the meteorological data does not follow a simple sinusoidal pattern, it is possible to decompose the signal into several Fourier terms where the previous description would be the main mode [23].
T n
z1 t
Several “soil resistances” are given in graphical form in [24] for different trench configurations based on this principle. The mean fluid temperature is found from the pipe resistance:
:mean radiative coefficient = hconv + hevap + hrad
rh :mean air relative humidity :empirical constant for the evaporation flux
1 T = t
z2 2
(7)
:mean solar radiation G :mean ground solar absortion coefficient
g
+ I
The response factor of the two-pipe arrangement can then be approximated as
(3)
Tamb :mean air temperature Tsky :mean sky temperature
hr
rp q I 2 ks 2 t
(6c)
where:
hrad he
=
(6b)
And E1 is the first-order exponential integral. The solution at the second pipe would be:
Fig. 2. Horizontal heat exchanger cross-section. 536
Applied Thermal Engineering 155 (2019) 534–545
L. Lamarche
Fig. 3. Horizontal heat exchanger equivalent problem.
tc, h =
10 z o2
The solution to Eq. (14) is given based on the dimensionless temperature:
(13b)
This can be seen in Fig. 4, where the steady-state thermal resistance for a single pipe is compared with the unsteady one derived from the ILS approach. Fontaine et al. [19] extended the approach of Claesson and Dunand to the unsteady case. In their approach, they could not solve the whole problem as a global unsteady thermal response factor. Instead, they split the horizontal piping into segments. Mutual thermal responses were evaluated with the finite line source (FLS) model. In their solution, the total heat exchange through the heat exchanger varied with time, but they assumed that, at a given time, the heat transfer distribution along the pipe followed the steady distribution found with the Claesson and Dunand distribution. In this paper, an approach similar to the one followed by Fontaine et al. is used, but no assumption is made regarding the heat distribution along the pipe. Instead, this distribution is part of the solution. The model is also based on the FLS model, but uses modified mathematical expressions of the FLS that were initially solved in the context of vertical boreholes and applied here in horizontal configurations [16,17,27]. The model also takes into account seasonal temperature variations that have a significant influence on the thermal behavior of horizontal geothermal fields. The model assumes, however, that the ground is homogeneous and single-phase.
T1 Tin
dT2 dx
= q2 =
+
(To T2) R2
(T1 T2) R12
To) = q1 R1 + q2 R12
(T2
To) = q1 R12 + q2 R2
x = m
(1 )
x ))
cosh( (1 x )) cosh( ) +
2)
m sinh (
(1 sinh ( ) m
x ))
x L
=
,
1
=
1+ 2
2
R1
,
,
2
=
R2 2 m
=
,
12
+2
=
R12
(19)
12 m
It is easy to verify that in the case of 1 = 2 , the classical Zeng’s profile [29] is obtained. The exit temperature is then: o
=
T2 (0) To = Tin To
cosh( ) cosh( ) +
m sinh ( m sinh (
) )
(20)
From which an effective resistance can be found:
R eq3 =
R1 + R2 2
(1 + o) = (1 o)
cotanh( )
(21)
This is the classical result found in a U-tube borehole if both resistances are equal, but with the difference that now the heat rate is per 0.6
0.5
(14)
0.4
0.3
0.2 Unsteady Steady
0.1
(15)
Inverting the last system will give: 2 2 2 R1 R2 R12 R1 R2 R12 R1 R2 R12 R1 = , R2 = , R12 = R2 R12 R1 R12 R12
m sinh (
(18b)
This system is similar to a typical U-tube model inside a vertical heat exchanger (see, for example [28]), except that the undisturbed ground temperature replaces the borehole temperature, and that the resistance associated with each pipe is now different. The Δ formalism is used as a reference in this case. Those resistances can easily be found from the values given by Eq. (10). In fact, we have:
(T1
m sinh (
with:
(T2 T1) R12
+
cosh( (1 x )) + cosh( ) +
2)
To x( 1 = exp To 2
Soil Resistance
mCp
(To T1) R1
(17)
To x( 1 = exp To 2
T2 Tin
Before looking at the transient case, it is interesting to revisit the steady case. The simple two-pipe case (Fig. 1) will be studied here. The Claesson and Dunand solution was found, assuming a linear variation of the temperature along the pipe. It is known that this is often a valid assumption, but for long pipes with low flow rate, it may lead to errors. The exact solution can be easily found in some simple configurations solving the following differential equations system: dT
To To
(18a)
3. New model for steady-state simulations
mCp dx1 = q1 =
T Tin
=
-4
-3
-2
-1
0
1
2
log(t/tc)
Fig. 4. Comparison between the unsteady and steady-state soil resistances for a single pipe for zo/rp = 50.
(16) 537
Applied Thermal Engineering 155 (2019) 534–545
L. Lamarche
unit length of pipe, and not per unit length of borehole.
Tb (t ) =
1 do / di + log di hconv 2 kpipe
T1
2
=
q q hFLS = ( (r , t , H1, H2, d2 2 ks 2 ks
(r , t , H1, H2, × [ierf(( z + H2) s )
(22)
d1) (27a)
22 1 e r s ds z) = 2H2 1/ 4 t s2
ierf( z s ) + ierf(( z
H1) s )
ierf(( z + H2
H1) s )]
T1
2, hor
=
q q hFLS, hor ,1,2 = ( (r1, t , H1, H2, x ) 2 k 2 k (28)
(r2, t , H1, H2, x )) where
r1 = rp
The method of images described in the previous section is applicable for the unsteady case, where it is again assumed that the heat flux is uniform along the pipe. In order to take into account the variation of the heat flux, Fontaine et al. suggested splitting the pipe into segments, with the heat transfer rate varying from segment to segment. At each time step, the total heat transfer is assumed to be known, and the variation along the pipe is assumed to vary exponentially. From the known heat transfer distribution along the pipe, the fluid temperature variation is calculated with the finite line source model (FLS). In its most general form, the FLS model is given [30] by (see Fig. 8):
,segment on the same pipe section
dA2
q
t
dA1
0
(
(t
d2 4 (t )
e
))3/2
di2 4 (t )
e
r1 = (y1
y2 ) 2 + (d1
r2 = (y1
,segment on the same pipe section y2 ) 2 + (d1 + d2 )2 ,otherwise
r2 = 2di
d2 ) 2
,otherwise
4.2. Unsteady temperature history Knowing the response factor associated with a horizontal piping arrangement, it is possible to calculate the time evolution of the average pipe temperature over a pipe segment. This calculation assumes knowledge of the heat flux associated with every segment. Again, in Fontaine’s model, the total heat transfer time history is known and the variation along the pipe is assumed to follow an exponential variation. Since the heat transfer varies with time, a convolution scheme has to be performed to estimate the temperature evolution. This convolution can
d
(23) A common simplification is often done by integration over time. Doing so will result in:
T1
2 (t )
=
q 4 kH2
dA2
dA1
erfc(d/2 d
t)
erfc(di /2 di
t) (24)
which is the form used by Fontaine et al. in their model. This model involves a double integral, which has been simplified in the case of vertical boreholes by Lamarche and Beauchamp [31] and by Claesson and Javed [27]. The latter solution is of particular interest. The idea is to keep the time integral and invert the order of integration. In the case of a vertical borehole, the solution becomes (see Fig. 9): Tb (t ) =
d+H
8k H
d
dz
H d
d
q
t 0
(
(t
))3/2
e
r 2+ (z 4 (t
)2 )
e
r 2 + (z + )2 4 (t )
d
(25) with
r2
=
x2
+
(27b)
The generalization of the segment-to-segment temperature response in the general inclined case was studied by Lazarotto [33], where the author simplified the double integral involved, but could not reduce it into a single integral in the general case. Although it may not seem possible to simplify the thermal response in the general case, it can be done in the horizontal configuration. Looking at Fig. 10, it can be shown that the generalization of Cimmino’s expression is given by (Appendix A):
4.1. Finite line source model for horizontal pipes
8 k H2
ds
x2
(r , t , H1, H2, H1 + d1 + d2))
4. Unsteady analysis
=
ierf(2ds )
with ierf(x ) = x ·erf(x ) (1 e ) Details of the derivation can be found in [17]. The generalization of the previous expression for the influence of a vertical segment on another vertical segment was given by Cimmino and Bernier [32]
The results obtained by the finite element simulation are shown in Fig. 6, where they are compared with the analytic variation given by Eq. (18) and the linear variation, assuming a soil resistance given by Claesson and Dunand (Eq. (9)) and the classical method of images of Bose et al. (Eq. (8)). In the latter case, the thermal response is unsteady, and only the steady-state solution is shown. The corresponding resistances are given in Table 2, and the corresponding inlet and outlet temperatures in Table 3. A grid sensitivity was performed and is shown in Fig. 7 where it is seen that the effective resistance converge towards the value given in Table 2.
2 ( t)
2 2 2 ierf(Hs ) + 2 ierf(Hs + 2ds ) ierf(2Hs + 2ds ) e r s s2
1
A simulation using the finite element software COMSOL Multiphysics™ was done to compare the different steady-state solutions for the case just described (see Fig. 5). The geometry is the same as the one in Fig. 1. The domain was split in half to take into account the symmetry of the problem. The parameters of the simulation are given in Table 1. No thickness was simulated on the pipe. Instead the convection coefficient represent and equivalent transfer coefficient:
T1
1/ 4 t
(26)
3.1. Simulation results
1 = do h
q 4 kH
Fig. 5. Mesh generated for horizontal ground heat exchanger.
y2 538
Applied Thermal Engineering 155 (2019) 534–545
L. Lamarche
Table 1 GHE configuration. rb(m)
L(m)
u(m/s)
z1(m)
z2(m)
ρCp(MJ/ m3-K)
q(W)
ks(W/mK)
h(W/ m2K)
0.02
100
0.3
2
3
3.93
1000
8
1000
(ξi,ηi,ζi)
ground surface di
(ξ,η,ζ) A1
(x,y,z) d
A2
Fig. 8. Segment-to-segment temperature thermal influence.
product in frequency domain. In order to use the FFT approach, the heat transfer history has to be known a priori, which is not a trivial factor, especially in the case of horizontal piping, where the seasonal temperature variation can interfere with the process. The proposed model follows the path already proposed by the author in the case of mixed vertical boreholes [16], where any arrangements of vertical borehole in series and in parallel can be simulated. The idea is to split a horizontal pipe into several sections, each of which can be seen as a pipe in series with the previous section. Referring again to the typical configuration shown in Fig. 1, splitting the pipe into ns segments can be analysed by the tools given, assuming ns pipes in series. Two major differences have to be confronted in the case of horizontal systems: the response factors have to be changed and the seasonal temperature variation has to be taken into account. The first problem has been covered in the previous section, and the second one will now be discussed. The description will be done with only one parallel path, but the general case is given in [16]. It can be shown that the pipe temperature at a given time can be evaluated with the one at a previous time step in a form given by:
Fig. 6. Fluid temperature along the heat exchanger. Table 2 Equivalent resistance (m-K/W). Req1 (Eq. (8))
Req2 (Eq. (9))
Req3 (Eq. (21))
Req, comsol
0.145
0.163
0.167
0.169
Table 3 Inlet and Outlet temperatures. Req1 Tfin
Tfout
Req2 Tfin
Tfout
Req3 Tfin
Tfout
Req, comsol Tfin Tfout
−1.08
−0.40
−1.15
−0.4
−1.17
−0.50
−1.18
−0.51
ns
Tp, i = Sp, i +
Sq, ij Xj ( Tfi, j
Tp, j )
j=1
(29)
where the expressions of Sp, i, Sq, ij are given in [16] and reproduced in Appendix B. We also have:
qi = Xi (Tfi, i X=
(mCp )i Hi
where o, i
=
o, i
Tp, i )
(1
o, i )
(30)
represents the temperature variation along the segment:
Tfo, i
Tp, i
Tfi, i
Tp, i
(31)
Assuming a linear variation, we have: o, i
=
1 Yi 1 + Yi
,
Yi =
Hi 2(mCp)i Rp
(32)
Assuming an exponential variation, the result is: o, i
Fig. 7. Impact of grid refinement.
= exp( 2Yi )
(33)
Again, these assumptions are for a single segment, not the entire pipe. In the case of np parallel pipes, where all the inlet temperatures are known, Eq. (29) represents a system of np equations for the np unknowns. For the case of ns series arrangement, only the first temperature is known and the (na = ns−1) others are unknowns. The new equations that have to be added to the system are:
be very time-consuming, and Fontaine uses the Fast Fourier Transform (FFT) approach proposed by Marcotte and Pasquier [34] to accelerate the convolution evaluation. The FFT is a very powerful approach that allows the time domain convolution product to be replaced by a simple 539
Applied Thermal Engineering 155 (2019) 534–545
L. Lamarche
Fig. 9. Vertical segment-to-segment thermal response factor.
Tfi, i = Tfo, i
1
=
o, i 1 Tfi, i 1
+ Tp, i
1 (1
o, i 1)
(1 (34)
A ×T =B A =
LL =
A (ns × ns ) UR (ns × na ) LL (na × ns ) LR (na × na)
, B =
BU BL
, T =
Tp TL
1 LR =
(35)
A=
1 + X1 Sq,11 , X2 Sq,12, ... X1 Sq,21 , 1 + X2 Sq,22, .... . . ..
X2 Sq,12, X2 Sq,22,
o,2
0
0 1 o,3
0 (1
0 . .. o,2 )
(39)
0 0 . .. 0 0. .. 1 0. .. (40)
Sp,1 + X1 Sq,11 Tfin,1 BU = Sp,2 + X1 Sq,21 Tfin,1
(36)
(41)
Tfi,2 TL = T fi,3
UR =
o,1)
0
(37)
BL =
X3 Sq,13, ... X3 Sq,23, ...
Tfi,1 0
(42)
Once the inlet temperatures are known, the heat transfer distribution is evaluated with Eq. (30), the Sp coefficients are updated (Eqs.
(38)
Fig. 10. Horizontal segment-to-segment thermal response factor. 540
Applied Thermal Engineering 155 (2019) 534–545
L. Lamarche
on a pipe segment becomes:
Tf , i
qi =
Tp, i
=
Rp
( 1Tp, i +
Tf , i
2T ) p, i
=
Rp
1T p, i
Tf , i
(54)
Rp
2T where T = T and 1T are given by Kusuda equation (Eq. (2)) or any other model used for the seasonal temperature variations. System of Eqs. (29) becomes:
(55)
1T ) p,i
qi = Xi (Tfi,i
ns 1T p, i
= Sp, i +
1T ) p, j
Sq, ij Xj (Tfi, j
(56)
j=1
which means that the results given in the previous section will give the pipe temperature contribution 1T and the reduced fluid temperature. Particular attention is paid to the complementary system:
Tfi, i = Tfi, i
Fig. 11. Fluid temperature during one year (test1).
Tfi, i
o, ns
=
Tfi,1
A+
Tfo, ns
Tp, ns
Tfi, ns
Tp, ns
o, ns
×T +
A+ =
Tfi, ns
Tfi,1 (44)
Tp, ns (1
q o, ns ) = mCp
,
B+ =
BU BL
,
T+ =
UR =
0
0
0 . .. (1
0
1
0 o,1
LR =
0
0 1 o,2
(1
o, ns )
o,1)
(50)
0 . .. o, ns 0 0. .. 1 0. .. (51)
Sp,1 BU = Sp,2 (52)
BL =
q/ mCp 0
2T p, i 1)
2T p, i
2T ) p, i
o, ns
1T p, ns
Tfi, ns
(1
o, ns )
=
q + ( 2Tp, ns mCp
2T ) p,1
In the first simulation, the inlet temperature will be imposed at the inlet of the upward pipe section shown in Fig. 1. A linear profile as shown In Fig. 11 will be imposed during half of the year, and a constant temperature is imposed for the rest of the year. In the numerical simulation, a time-varying temperature profile given by Eq. (1) is imposed at the surface of the domain (z = 0), with A = 10 and tshift = 0 (January 1st is the coldest day of the year). A 5-hour time step was chosen for the numerical simulation, while a 1- hour-time step is taken during the analytic simulation. Eight segments (4 on each pipe section) were chosen for the analytic simulation. Fig. 11 shows the outlet temperature variation during the year, and Fig. 12, the total heat exchange with the ground. Fig. 13 shows the temperature variation and Fig. 14 the heat transfer per unit length along the pipe for two different simulation times: after half a year and after one year. The impact of the number of segments was studied for this particular test. The number of segments in each pipe section was varied and the difference between the total heat exchange in the ground was compared. Thus difference is defined as:
X2 Sq,12, ... X2 Sq,22, ... (49)
LL =
=(
+
2T p, i 1
6.1. Inlet temperature imposed
(48)
X1 Sq,11, X1 Sq,21,
o , i 1)
1T p, i 1
The model just described was compared with a finite element simulation on COMSOL™ Multiphysics. The geometry studied is the same as the one shown in Fig. 1, and the grid is the same as in the steady-state analysis.
Tp TL (47)
Tfi,1 TL = T fi,2
(1
o, i 1)(
o, i 1)
1T p, i 1
6. Simulation results
(46)
A (ns × ns ) UR (ns × ns ) LL (ns × ns ) LR (ns × ns )
o, i 1 Tfi, i 1
+ (1
1 (1
where the last term is zero again if the inlet and outlet are at the same depth.
(45)
=B+
o, i 1 Tfi, i 1
+ Tp, i
The last term in Eq. (57) is zero if the two segments are at the same depth, but different from zero otherwise. Eq. (45) should also be modified:
(43)
Tfo, ns )
=
o, i 1 Tfi, i 1
(57)
(A6)–(A7)), and the scheme moves to the next time step. In the case where the inlet temperature is not known, but the total heat transfer is imposed, the system is modified:
q = (mCp)(Tfi,1
2T 1p, i
(53)
Qanalytical =
qanalytical (t ) dt
Qcomsol =
qcomsol (t ) dt |Qanalytical Qcomsol |
5. Impact of seasonal temperature variation
difference (%) =
The effect of seasonal variation is quite simple when the heat flux is known, and less trivial when it is part of the solution. The heat transfer
between the analytical model and the numerical simulation was evaluated and the result is shown in Fig. 15. 541
Qcomsol
× 100
(58)
Applied Thermal Engineering 155 (2019) 534–545
L. Lamarche
Fig. 12. Total heat transfer during one year (test1).
Fig. 15. Influence of the number of segments on the difference between analytical and numerical solution.
6.2. Total heat transfer imposed As a second test, the total heat flux was imposed to test the systems of Eqs. (46)–(53). A constant heat extraction of 5000 W was imposed for the first part of the year (t < 4380 h), while no heat exchange (restitution period) was imposed for the rest of the year. All other parameters were the same as the one given in Table 4. The inlet and outlet temperatures are shown in Fig. 16. The temperature and heat flux along the pipe are now given just before the end of the heat injection (t = 4375 h), and 100 h after the power is turned off (t = 5380 h), and shown in Fig. 17. The heat transfer per unit length is shown in Fig. 18. It can be seen that when the total heat exchange is zero, the local heat transfer is not, due in part to the temperature difference between the ground temperature around the pipe. Finally, the calculation times for both approaches are compared in Table 5, where an important reduction in the calculation time is observed using the proposed model.
Fig. 13. Temperature variation along the pipe for two different times (test1).
6.3. Comparison with simpler models As stated in Section 2, the ILS model was often used to generate unsteady soil resistances for horizontal piping. In the case where the fluid temperature is imposed, it is not trivial to apply the model but in the case where the heat transfer is imposed, superposition of the solution given by the method of images and the seasonal variation is possible. For example the second test analyzed in the last subsection can be compared with the method of images. From Eq. (8): 1T f
= q Rs + Rp
With Rs given by Eq. (7) for the two-pipe arrangement. The inletoutlet temperatures are then 1T f ,o
Tf , o =
=
1T f
1T f ,o
+ +
q 2mCp
,
2T , p, o
Tf , i =
1T f ,o
=
1T f ,i
1T f
q 2mCp
+ , 1Tp, i
(59a) (59b)
Comparison of the inlet-outlet temperatures thus obtained is compared with the numerical simulation in Fig. 19. A zoom in the vicinity
Fig. 14. Heat transfer per unit length variation along the pipe for two different times (test1).
542
Applied Thermal Engineering 155 (2019) 534–545
L. Lamarche
Table 4 GHE configuration. rb(m)
L(m)
u(m/s)
z1(m)
z2(m)
(ρCp)f(MJ/m3-K)
ks(W/mK)
(ρCp)s(MJ/m3-K)
h(W/m2K)
ns
0.02
100
0.3
2
3
3.93
8
2.0
100
8
Fig. 16. Fluid temperature during one year (test2).
Fig. 18. Heat transfer per unit length variation along the pipe for two different times (test2). Table 5 CPU time. Test1 (Comsol)
Test1 (Analytic)
Test2 (Comsol)
Test2 (Analytic)
18 min
2.2 sec
21 min
2.2 sec
Fig. 17. Temperature variation along the pipe for two different times (test2).
of the end of the heat pulse is done to better see the difference in the models. Although the general trend of the temperature evolution is correctly predicted with the ILS model, larger variations are observed near the regions where the heat flux varies. Most of the differences come from the assumption that the heat transfer per unit length is uniform. The rigorous approach was described in the previous sections. However a simple variation of the method of images using the ILS can be done if we allow the heat transfer to vary from one pipe to the other. A possibility is to assume that, even in the unsteady case, the heat transfer distribution resembles the steady case given by the Claesson and Dunand [4] relation (Eq. (12)). Rewriting those equations
R R12 /2 q1 = 2q 2 = R1 + R2 2R12
1q
R R12 + /2 q2 = 2q 1 = R1 + R2 2R12
Fig. 19a. Outlet temperature given by Eq. (59)
Rs =
+
(61)
2 R2
where R1 and R2 are given by Eq. (6). Replacing this new resistance in Eq. (59) 1T f ,o
= q Rs + Rp +
2q
Tf , o =
(60) as:
1 R1
1T f ,o
+
2T p, o
q 2mCp
,
,
Tf , i =
1T f ,o
= q Rs + Rp
1T f ,i
+ , 2Tp, i
q 2mCp
(62a) (62b)
would give the results shown in Fig. 20: It is seen that even though, the unsteadiness is still not well
The unsteady soil resistance given by Eq. (7) can now be rewritten 543
Applied Thermal Engineering 155 (2019) 534–545
L. Lamarche
Fig. 19b. Inlet temperature given by Eq. (59).
Fig. 20b. Inlet temperature given by Eq. (62).
source associated to horizontal configurations was developed to simulate the heat transfer between a horizontal heat exchanger and the surrounding ground. The model was compared to a finite element simulation in the case of a simple configuration and display excellent agreement given that the model is 500 0times faster than the numerical simulation. While the configuration may be simple, it illustrates a very important aspect of horizontal systems, namely, different local ground temperatures around pipes at different heights and how this can affect the thermal behavior of the ground exchanger. The model can easily be extended to different inline configurations, which can have parallel branches as well. Extensions to slinky or spiral configurations can also be considered, but in that case, the thermal response factor between pipe sections would be more complex. An extension of the classical work of Claesson and Dunand was also presented as part of this study. It is a future goal to use it to provide potentially better guidelines for horizontal design procedures. Acknowledgements
Fig. 20a. Outlet temperature given by Eq. (62).
We would like to thank the Natural Sciences and Engineering Research Council of Canada (NSERC), which partly financed this research (RGPIN-2014-06240).
predicted, the overall comparison is better. 7. Conclusion An analytical model based on a new formulation of the finite line Appendix A Referring to Fig. we have:
q Ti j (t ) = 8k H
x + H2 x
dx
H1 0
1
t
d
0
(
))3/2
(t
r 2 + (x 4 (t
e
)2 )
e
2 rim + (x + )2 4 (t )
d (A1)
where r represents the perpendicular distance between the source line and the surface of interest. In the case where two segments are on the same pipe section, we then have:
r = rp
,
rim = 2z i
In the case where the two segments are in different pipe sections:
r2
= (yi
yj ) 2 + (z i
zj ) 2
,
2 rim = (yi
yj ) 2 + (z i + zj ) 2
Using the usual change of variable s = 1/ 4 (t
Ti
j (t )
=
q 2 4 k H2
1/ 4 t
[e
r 2s 2
e
ri2 s 2] d
s
x + H2 x
(A2)
) leads to: dx
H1 0
e
s 2 (x
)2 d
544
(A3)
Applied Thermal Engineering 155 (2019) 534–545
L. Lamarche
Ti
j (t )
=
q 4 k H2
Ti
j (t )
=
q 4 k H2
1/ 4 t
r 2s 2
[e
2 2 e ri s
22 e r s
s2
1/ 4 t
ierf(( x
e
ri2 s 2]
ds s
x + H2 x
(erf(sx )
(ierf(( x + H2) s )
H1) s ) ierf(( x + H2
erf(s (x
H1)) dx
(A4)
ierf( x s )+ (A5)
H1) s ) ds
Appendix B
t =
Sp, i =
t r p2 1 ks
Fi ( t +
Nz
e
zn2 t F ( t i
, zn) z n
, Sq, ij =
n=1
t , z n) = e
zn2 t F
i(t
, z n) +
nserie j=1
1 ks
Nz
(1
e
zn2 t
) uij (z n) zn
(A6)
n= 1
qj ( t +
t )(1
e
zn2 t
) uij (z n) (A7)
Fi ( t = 0, zn ) = 0 uij (z ) =
z
L 1 (hFLS, hor , ij )
(A8)
Details on the derivation of these expressions are given in [16]. The choice of the integral variable z is discussed in [35]. Appendix C. Supplementary material Supplementary data to this article can be found online at https://doi.org/10.1016/j.applthermaleng.2019.04.006.
ground heat exchangers, Int. J. Heat Mass Transf. 115 (2017) 354–360. [19] P.-O. Fontaine, D. Marcotte, P. Pasquier, D. Thibodeau, Modeling of horizontal geoexchange systems for building heating and permafrost stabilization, Geothermics 40 (3) (2011) 211–220. [20] T. Kusuda, P.R. Achenbach, Earth temperature and thermal diffusivity at selected stations in the United States, National Bureau of Standards Gaithersburg MD (1965). [21] M. Badache, P. Eslami-Nejad, M. Ouzzane, Z. Aidoun, L. Lamarche, A new modeling approach for improved ground temperature profile determination, Renew. Energy 85 (2016) 436–444. [22] G. Mihalakakou, On estimating soil surface temperature profiles, Energy Build. 34 (3) (2002) 251–259. [23] C. Jacovides, G. Mihalakakou, M. Santamouris, J. Lewis, On the ground temperature profile for passive cooling applications in buildings, Sol. Energy 57 (3) (1996) 167–175. [24] J.E. Bose, J.D. Parker, F.C. McQuiston, American Society of Heating R, Engineers AC. Design/data Manual for Closed-loop Ground-coupled Heat Pump Systems: American Society of Heating, Refrigerating, and Air-Conditioning Engineers, 1985. [25] H.S. Carslaw, J.C. Jaeger, Conduction of Heat in Solids, 2nd ed., 1959. [26] P. Eskilson, Thermal analysis of heat extraction boreholes, 1987. [27] J. Claesson, S. Javed, An analytical method to calculate borehole fluid temperatures for time-scales from minutes to decades, in: Conference An Analytical Method to Calculate Borehole Fluid Temperatures for Time-scales from Minutes to Decades, vol. 117, pp. 279–288. [28] L. Lamarche, S. Kajl, B. Beauchamp, A review of methods to evaluate borehole thermal resistances in geothermal heat-pump systems, Geothermics 39 (2) (2010) 187–200. [29] H. Zeng, N. Diao, Z. Fang, Heat transfer analysis of boreholes in vertical ground heat exchangers, Int. J. Heat Mass Transf. 46 (23) (2003) 4467–4481. [30] L. Lamarche, Analytical g-function for inclined boreholes in ground-source heat pump systems, Geothermics 40 (4) (2011) 241–249. [31] L. Lamarche, B. Beauchamp, A new contribution to the finite line-source model for geothermal boreholes, Energy Build. 39 (2) (2007) 188–198. [32] M. Cimmino, Fast calculation of the g-functions of geothermal borehole fields using similarities in the evaluation of the finite line source solution, J. Build. Perform. Simul. 1–14 (2018). [33] A. Lazzarotto, A methodology for the calculation of response functions for geothermal fields with arbitrarily oriented boreholes b, Renew. Energy 86 (2016) 1380–1393. [34] D. Marcotte, P. Pasquier, Fast fluid and ground temperature computation for geothermal ground-loop heat exchanger systems, Geothermics 37 (6) (2008) 651–665. [35] L. Lamarche, B. Beauchamp, A fast algorithm for the simulation of GCHP systems, ASHRAE Trans. 113 (1) (2007) 470–476.
References [1] J.W. Lund, T.L. Boyd, Direct utilization of geothermal energy 2015 worldwide review, Geothermics 60 (2016) 66–93. [2] J. Raymond, M. Malo, D. Tanguay, S. Grasby, F. Bakhteyar, Direct utilization of geothermal energy from coast to coast: a review of current applications and research in Canada, in: Conference Direct Utilization of Geothermal Energy from Coast to Coast: A Review of Current Applications and Research in Canada. [3] S. Rees, Advances in Ground-source Heat Pump Systems, Woodhead Publishing is an imprint of Elsevier, Duxford, UK, 2016. [4] J. Claesson, A. Dunand, Heat Extraction from Ground by Horizontal Pipes, Swedish Council of building Research, Stokholm, 1983. [5] V.C. Mei, Horizontal ground-coil heat exchanger theoretical and experimental analysis, United States: Report No. ORNL/CON-193, Dec, Oak Ridge National Laboratory, Oak Ridge, 1986. [6] M. Piechowski, Heat and mass transfer model of a ground heat exchanger: theoretical development, Int. J. Energy. Res. 23 (7) (1999) 571–588. [7] K. Kupiec, B. Larwa, M. Gwadera, Heat transfer in horizontal ground heat exchangers, Appl. Therm. Eng. 75 (2015) 270–276. [8] N. Kayaci, H. Demir, Numerical modelling of transient soil temperature distribution for horizontal ground heat exchanger of ground source heat pump, Geothermics 73 (2018) 33–47. [9] P. Congedo, G. Colangelo, G. Starace, CFD simulations of horizontal ground heat exchangers: A comparison among different configurations, Appl. Therm. Eng. 33 (2012) 24–32. [10] G. Gan, Dynamic thermal modelling of horizontal ground-source heat pumps, Int. J. Low-Carbon Technol. 8 (2) (2013) 95–105. [11] S.K. Soni, M. Pandey, V.N. Bartaria, Ground coupled heat exchangers: A review and applications, Renew. Sustain. Energy Rev. 47 (2015) 83–92. [12] Z. Xiong, D.E. Fisher, J.D. Spitler, Development and validation of a Slinkyb “ ground heat exchanger model, Appl. Energy 141 (2015) 57–69. [13] S. Yoon, S.-R. Lee, G.-H. Go, Evaluation of thermal efficiency in different types of horizontal ground heat exchangers, Energy Build. 105 (2015) 100–105. [14] H. Fujii, K. Nishi, Y. Komaniwa, N. Chou, Numerical modeling of slinky-coil horizontal ground heat exchangers, Geothermics 41 (2012) 55–62. [15] H. Esen, M. Inalli, M. Esen, K. Pihtili, Energy and exergy analysis of a groundcoupled heat pump system with two horizontal ground heat exchangers, Build. Environ. 42 (10) (2007) 3606–3615. [16] L. Lamarche, Mixed arrangement of multiple input-output borehole systems, Appl. Therm. Eng. 124 (2017) 466–476. [17] M. Cimmino, M. Bernier, A semi-analytical method to generate g-functions for geothermal bore fields, Int. J. Heat Mass Transf. 70 (2014) 641–650. [18] L. Lamarche, G-function generation using a piecewise-linear profile applied to
545