Heat transfer in a borehole heat exchanger: Frequency domain modeling

International Journal of Heat and Mass Transfer 69 (2014) 129–139

Contents lists available at ScienceDirect

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

Heat transfer in a borehole heat exchanger: Frequency domain modeling Griet Monteyne a,⇑, Saqib Javed b, Gerd Vandersteen a a b

Vrije Universiteit Brussel, Department of ELEC, Pleinlaan 2, 1050 Brussel, Belgium Chalmers University of Technology, SE-412 96 Göteborg, Sweden

a r t i c l e

i n f o

Article history: Received 4 April 2013 Received in revised form 8 October 2013 Accepted 9 October 2013

Keywords: Heat transfer Modeling Borehole heat exchanger Frequency domain Modeling Ground Source Heat Pump

a b s t r a c t This paper proposes a new frequency domain method to model the heat transfer between the injected/ extracted heat and the temperature of the fluid exiting a borehole heat exchanger. The method is based on in situ measurements and focuses particularly on the short-term borehole heat transfer. It uses a rational function of the Warburg variable in the Laplace domain to model the borehole heat transfer. The rational model is transformed to a time domain model using inverse Laplace transformation. This time domain model makes it possible to calculate the temperature response on a random heat variation signal. The paper also demonstrates a new way to perform the classical thermal response test. Instead of injecting a constant amount of heat, the experiments have been performed using multiple short-duration heat injections. In this way, the obtained rational heat transfer model contains information about both the short- and the long-term heat transfer. The results obtained using the proposed modeling method are compared with those obtained from a state-of-the-art analytical method. The time domain model can be used to design a controller to optimize the performance of a Ground Source Heat Pump system, the efficiency of which depends strongly on the temperature of the fluid exiting the borehole. Ó 2013 Elsevier Ltd. All rights reserved.

1. Introduction A Ground Source Heat Pump (GSHP) system provides space heating and cooling by heat extraction from or injection into the ground. The heat exchange with the ground is often realized through a multitude of vertical borehole heat exchangers, which form a so-called borehole field. The energy efficiency of a GSHP system depends on the temperature of the fluid entering the GSHP at the source side. In other words, the temperature of the fluid exiting the borehole heat exchanger (BHE) influences the energy efficiency of the GSHP system. Thus, knowing this temperature is important to develop a control strategy that minimizes the energy consumption. Consequently, this paper aims to find an accurate model for the relation between the ground thermal loads (heating or cooling) and the temperature exiting the borehole. Over the years, numerous controllers and control schemes have been developed for adaptive control and dynamic optimization of GSHP systems. The objectives of these developments range from improving thermal comfort, optimizing borehole heat injection/ extraction, and maximizing heat pump performance to reducing operational costs and minimizing energy consumption. Gao et al. [1] showed that an effectively controlled intermittent heat injection/extraction process can optimize the efficiency of a GSHP ⇑ Corresponding author. Tel.: +32 2 6292949; fax: +32 2 6292850. E-mail address: [email protected] (G. Monteyne). 0017-9310/$ - see front matter Ó 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.ijheatmasstransfer.2013.10.015

system. Salque et al. [2] used a neural predictive controller to minimize the energy consumption while maintaining an acceptable level of thermal comfort. They used an artificial neural network (ANN) model to predict fluid temperature exiting the BHE. Gang and Wang [3] also used ANNs to develop static and dynamic models of BHEs. Verhelst [4] presented a model predictive control (MPC)-based controller to optimize the operation of a hybrid-GSHP system. She modeled the borehole heat transfer by white- and gray-box representations of a 1-dimensional finite-difference model using data reduction and parameter estimation approaches, respectively. Beck et al. [5] and de Paly et al. [6] used linear programming based control strategies to maximize the heat pump performance. They optimized the operation of a multiple BHE using optimized loads for each borehole at each time-step. They used infinite line-source solution [7] coupled with spatial and temporal superposition techniques to model the BHE. De Ridder et al. [8] used a dynamic programming controller for optimizing energy delivered by the borehole system while simultaneously minimizing the operational costs. They simulated a borehole thermal storage system using a simple first order response model. Most controllers and algorithms developed for GSHP system applications use fairly simple approaches to model borehole heat transfer. The heat transfer in a BHE is too complicated a process to be captured accurately by a simple mathematical model. The quality of the heat transfer model is vital to the quality of the controller using the model. A controller based on an oversimplified

130

G. Monteyne et al. / International Journal of Heat and Mass Transfer 69 (2014) 129–139

borehole heat transfer model will be unable to achieve the desired performance. This paper presents a new method to model the borehole heat transfer, which can be used to design a controller that optimizes the performance of a geothermal system. The proposed frequency domain model (Section 4) focuses particularly on the short-term borehole heat transfer. It approximates the relation between the heat and the temperature using a rational Frequency Response Function (FRF) model in the Laplace variable s or in the Warburg pffiffi variable s. The choice of the Warburg rational FRF model is based on a comparison with a state-of-the-art analytical solution [9] for the heat transport near a BHE. The FRF model is used to calculate the temperature variation of the fluid exiting the borehole for a prescribed sequence of heat injection pulses. The rational model is transformed to an impulse response in the time domain using inverse Laplace transformation. This paper also explores a new way of performing thermal response tests (TRTs). These tests are often used to acquire knowledge of the thermal characteristics of the borehole and its surrounding ground. The thermal characteristics, which relate the thermal load (heat) to the temperature, are required by traditional borehole heat transfer models as inputs. The TRTs in this study have been performed using a series of short-duration heat injections. In this way, the obtained measurement data contain more information about the short-term heat transfer than obtained from a classical TRT. This paper first presents an overview of several noteworthy borehole heat transfer models in Section 2. Next, the theory for the current state-of-the-art analytical model and the new model proposed in this paper are discussed in Sections 3 and 4, respectively. This theory is then applied to estimate the best rational model for the heat transfer near a borehole. The experimental setup and three experiments performed to estimate and validate the models are discussed in Section 5. Section 6 discusses and compares the results.

2. Borehole heat transfer models Several models have been developed to model the relation between ground loads and temperature of the fluid exiting a BHE. These models vary from complex numerical three-dimensional models to simple analytical solutions. The numerical methods require extensive computation and are used mostly in research and for complex design problems. The analytical models require short computation times and are mainly used for simple design purposes. Noteworthy numerical methods include the work of [10–14], among others. Eskilson [15] developed non-dimensional thermal response functions (a.k.a. g-functions) for multiple borehole fields that relate the temperature response to a unit step heat input based on a numerical approach. These g-functions are only valid for long time scales (t > 5r2b =a, where rb is the radius of the borehole and a is the thermal diffusivity of the surrounding soil). Yavuzturk et al. [16] extended the g-functions to shorter time scales. Classical analytical solutions include infinite line-source [7], and cylindrical-source [17] solutions. However, in recent years there has been an increased interest in analytical modeling of borehole heat transfer and many new analytical and semi-analytical solutions have been developed and presented. These include works of [18–21], among others. These analytical methods model the borehole as a line- or a cylindrical-source to obtain simple analytical solutions. This makes these techniques well suited to model the long-term response of the ground surrounding the borehole. However, the short-term response from these models is less

accurate. Recently, a new analytical model has been presented that also accounts for the short-term response [9]. This state-of-the-art analytical model will be briefly discussed in Section 3. Besides purely analytical and numerical models, models based on electrical analogy, with lumped capacities and thermal resistances, have also been developed [22–25]. These models require less computation time than the numerical models and in some cases, also account for the short-term response of the BHE. 3. Current state-of-the-art analytical model for heat transfer When designing controllers for GSHP system applications, analytical methods, because of their superior computational efficiencies, are often used to model the borehole heat transfer. This section describes the current state-of-the-art analytical model for heat transfer near a BHE [9]. The frequency domain model proposed in this paper will be compared to this analytical model later. The state-of-the-art model is based on the analytical solution of the radial heat equation in the two-dimensional infinite plane perpendicular to the borehole axis, viz:

  1 @ @Tðr; tÞ @Tðr; tÞ ¼ qðrÞC p ðrÞ kðrÞr r @r @r @t

ð1Þ

where r represents the radial distance to the center of the borehole, t the time, T(r, t) the temperature, k the thermal conductivity, q the density, and Cp the heat capacity of the surrounding medium. 3.1. Basic idea of current state-of-the-art model The analytical solution of the heat conduction Eq. (1) is determined for the following initial and boundary conditions (Fig. 1): 1. The U-tube is modeled as an equivalent-diameter pipe of radius rp. 2. The circulating fluid in the equivalent-diameter pipe has temperature Tf(t) and thermal capacity Cp in the pipe. 3. The combined thermal resistance of the circulating fluid and the pipe is Rp. 4. A grout region surrounds the equivalent-diameter pipe and has a thermal conductivity kb and thermal diffusivity ab. 5. The borehole is surrounded by homogeneous ground with thermal conductivity k and thermal diffusivity a. 6. All heat fluxes and temperatures are at rest at time t = 0. 7. A constant heat flux qinj is injected to the circulating fluid from time t = 0. 8. Heat flux qp(t) is transferred from the fluid to the grout region through the pipe wall, and heat flux qb(t) is transferred across the borehole boundary to the surrounding ground. 9. The resulting outer boundary temperatures of the equivalent-diameter pipe and the borehole are Tp(t) and Tb(t), respectively.

Ground λ, a

T(r,t)

Grout λb , ab qp(t)

Fluid Tb(t)

Tf (t) Cp qinj

rp

qb(t)

rb

r

Rp Tp(t) Fig. 1. Heat transfer in the equivalent diameter pipe and its surrounding regions.

G. Monteyne et al. / International Journal of Heat and Mass Transfer 69 (2014) 129–139

The relation between the increase of the mean fluid temperature ðDT f ;mean Þ and the injected heat (qinj) is determined in the Laplace domain. The solution first computes the Laplace transforms of the heat transfer for the pipe, grout and the ground regions and obtains equations between the Laplace transforms of the boundary temperatures and the boundary fluxes for each region. The equations are then represented in the form of a thermal network for a borehole in the ground (Fig. 2). The network involves two absorptive resistances Rp ðsÞ and Rb ðsÞ, one transmittive resistance Rt ðsÞ and a ground thermal resistance Rg ðsÞ in the Laplace domain. The relation between the Laplace transform of the fluid temperature increase ðT f ;mean ðsÞÞ and the injected heat (qinj/s) is then obtained by solving the following thermal network.

T f ;mean ðsÞ ¼

qinj s

1

ð2Þ

1

Cps þ

1

Rp ðsÞ þ 1 Rp ðsÞ

1 Rg ðsÞRb ðsÞ

þ Rt ðsÞ þ

Rg ðsÞ þ Rb ðsÞ

2

p

Z 0

1

1  eu u

2 t=t

0

LðuÞdu

ð3Þ

with,

  LðuÞ ¼ Imag sT f ;mean ðsÞ C

C : t 0 s ¼ u2 Here, T0 is the mean fluid temperature at the start of the experiment, t0 is an arbitrary time constant and the function L(u) is obtained by taking the Laplace transform of the fluid temperature increase (2) for s on the negative real axis ðCÞ. More details can be found in [9]. The solution accounts for time varying loads by considering piece-wise constant heat pulses for each time step. The mean fluid temperature at any given time is obtained by superposition of the solution from each of the preceding pulses. The borehole exit fluid temperature is calculated from the mean fluid temperature as:

T f ;out ðtÞ ¼ T f ;mean ðtÞ þ

Q in ðtÞ 2C f V f ðtÞ

ð4Þ

Here, Qin(t) is the prescribed injection rate in Watts, Cf(J/m3 K) is the volumetric heat capacity of the circulating fluid and Vf(t) is the volumetric flow rate in m3/s. Limitations of the analytical solution include the assumption of radial heat transfer in a borehole. Also, similar to other 1- and 2dimensional models, the solution does not account for residence time of the fluid in the borehole. Hence, time steps shorter than the residence time of the fluid in the borehole introduce an error in the output of the analytical model. Residence time of the fluid in a BHE depends on its velocity. It is of the order of a few minutes qinj s

Tf (s)

Rp

Tp (s)

qb (s)

qp (s) 1 Cp·s 0

Rg (s)

Rt (s) Tb (s)

Rp (s)

0

and typically ranges between 5 and 20 min for 100–300 m deep boreholes. 4. Proposed modeling technique for heat transfer This section describes the new heat transfer modeling method in the frequency domain. The basic idea and theory behind this proposed frequency domain modeling method are explained in Sections 4.1 and 4.2. This frequency domain model is then inverse-transformed into an impulse response in the time domain (Section 4.3). The impulse response is then used to calculate the fluid temperature variations caused by the heat injection variations. This impulse response can be used to optimize a controller, for example, using Model Predictive Control [26]. Note that this impulse response is equivalent to the derivative of the g-functions from Eskilson’s numerical model (Section 2). However, the proposed method obtains this impulse response directly via measurements, while for the numerical methods, the g-functions are derived from numerical simulations. 4.1. Basic idea of proposed method

Finally, the mean fluid temperature increase in time domain is obtained by using the following standard inversion formula.

DT f ;mean ðtÞ ¼ T f ;mean ðtÞ  T 0 ¼

131

0

Rb (s)

0

Fig. 2. Thermal network representation of the heat transfer from the fluid to the surrounding ground (in Laplace domain).

The basic idea of the proposed method is based on the fact that the relation between the injected/extracted heat and the temperature of the fluid (2) can be approximated by a rational Frequency Response Function (FRF) model (G1 or G2) in the Laplace variable pffiffi s (5) or in the Warburg variable s (6), viz:

 Qn z    s  zig G1 s; hzg ; hp ¼ c1 Qi¼1 np i¼1 ðs  pi Þ

ð5Þ

 Qnz pffiffi pffiffi  s  zig i¼1  G2 s; hzg ; hp ¼ c2 Qnp pffiffi s  pi i¼1

ð6Þ

where c1 and c2 represent constants, zig the zeros of the rational model, pi the poles of the rational model, and np and nz the number of poles and zeros of the rational models G1(s, h) or G2(s, h). The poles are grouped in the parameter vector hp, while the constant and the pffiffi zeros are grouped in hzg. In Eq. (6), a rational model in s is considpffiffi ered since Rt ; Rp ; Rg and Rb are functions of s. The advantage of these rational frequency domain models is that they are not based on strict approximations as opposed to the current state-of-the-art modeling methods. In fact, these models try to find the best relation between the injected/extracted heat and the temperature increase of the fluid by evaluating the measurement data. These rational models are applicable to different borehole types (e.g. single U-tube, double U-tube, coaxial pipe) and any borehole field configuration. The model also implicitly includes convective heat transfer in a BHE and/or advective heat transfer around the borehole and even natural convection in annulus region of groundwater-filled boreholes. However, the rational model is only valid for operational setups that correspond to the experimental setup. In other words, the obtained rational model only reflects what was present during the experiment used for modeling. The transport time (tdelay) of the fluid in the pipe from the position of the heat input to the position of the temperature measurement has to be included into the rational models, viz:

Gdelay ðs; hzg ; hp Þ ¼ G1 ðs; hÞetdelay s 1 pffiffi  pffiffi  s; hzg ; hp ¼ G2 s; h etdelay s Gdelay 2

ð7Þ

In order to simplify the analysis, the mass flow rate is not integrated into this model. However, the model can be extended to take variations of mass flow rate into account. In that case, the flow

132

G. Monteyne et al. / International Journal of Heat and Mass Transfer 69 (2014) 129–139

rate would act as a varying parameter in the obtained heat transfer model. But, the experiment must be redesigned in order to include flow rate variations. 4.2. Model estimation via proposed method in practice

ð8Þ

where XðkÞ ¼ Q ðkÞ or DT f ;out ðkÞ represents the DFT spectrum of q(t) or (DTf,out(t) = Tf,out(t)  Tf(0)) at frequency line k (s = j2pk), and N represents the number of measurement samples. It is assumed that there is no initial temperature gradient (i.e. Tf(0) is the same everywhere at the start of the experiment). The relation between the injected/extracted heat and the temperature of the fluid exiting the borehole can be written in the frequency domain as follows [27]:

DT f ;out ðkÞ ¼ Gdelay ðkÞQðkÞ þ Lx ðkÞ þ NT f ;out ðkÞ x Gdelay x

Gdelay 1

ð9Þ

Gdelay 2

where ¼ or represent the rational heat transfer model, Lx the transient, and N T f ;out the noise on the measured temperature signal. The heat transfer model Gdelay corresponds to the x complementary function of the partial differential heat equation in the time domain. The transient Lx is related to the difference between the initial and final conditions of the experiment [28]. For example, the fluid temperature at the start of the experiment differs from the fluid temperature at the end of the experiment. The best rational (parametric) model for each type (Gdelay and 1 delay G2 ) is found by estimating the transient and the heat transfer system simultaneously. A parametric rational model for both the transient and the heat transfer system can be represented as

  Gdelay s; hzg ; hp 1

t delay s

¼ ce

Qnzg 



ð10Þ

or,

 Qnzg pffiffi pffiffi  s  zig i¼1 t delay s  s; hzg ; hp ¼ ce Qnp pffiffi s  pi i¼1  Qnzt pffiffi pffiffi  s  zit  L2 s; hzt ; hp ¼ cl etdelay s Qi¼1 np pffiffi s  pi i¼1

ð12Þ



e

^hzg ; ^hzt ; ^hp





DT f ;out ðkÞ  Gdelay k; ^hzg ; ^hp QðkÞ  Lx k; ^hzt ; ^hp x ¼ sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  ffi 

2

½k r2 ðkÞ þ Gdelay k; h^zg ; h^p  r2 ðkÞ  2real Gdelay k; h^zg ; h^p r2



DT f ;out

x

Q

x

Q DT f ;out

ð13Þ where Q ðkÞ and DT f ;out ðkÞ represent the DFT spectra (8) of the input q(t) and output Tf(t) signals at frequency line k. r2Q ðkÞ and r2DT ðkÞ f ;out

represent the variance and

r2Q DT

the covariance of the input and f ;out

output signal. It is assumed that the noise on the temperature signal is Gaussian distributed. The minimum of the cost function is found using the method of Levenberg–Marquardt [29]. The different orders are estimated for both Gdelay and Gdelay . The 1 2 best model is selected out of these models by taking the following criteria into account: 1. 2. 3. 4.

The model is stable. The mean model error is small. The standard deviation of the rational model is small. Akaike Information Criterion (AIC) and Minimum Description Length (MDL) are optimal. The AIC minimizes the one-step ahead prediction error and is optimal for control purposes. However, AIC tends to select model orders that are too high. The MDL is better suited to obtain a physical model [31]. Hence, MDL is expected to be superior for long-term predictions.

Out of the models that satisfy these criteria, the model with the smallest order is selected. This proposed modeling technique has the advantage that it can handle different heat input signals. As long as the input signal cannot be presented easily by a rational form with an order similar to the order of the transient, a rational heat transfer model can be estimated. Typically, the on/off cycle of a GSHP under normal working conditions has enough information to obtain a rational heat transfer model. 4.3. Time domain equivalent of the frequency domain models

s  zig Qi¼1 np ð i¼1 s  pi Þ

Qnzt   ðs  zit Þ L1 s; hzt ; hp ¼ cl etdelay s Qi¼1 np i¼1 ðs  pi Þ







V SML h^zg ; ^hzt ; ^hp ¼ eH ^hzg ; ^hzt ; ^hp e ^hzg ; ^hzt ; ^hp with,

The following paragraphs explain the procedure to estimate the rational models (7) in the frequency domain. These rational model should describe the relation between the injected/extracted heat and the temperature of the fluid exiting the borehole. First, the initial temperature value Tf(0) is subtracted from the measured temperature. Second, the time domain data are transformed into the frequency domain using the Discrete Fourier Transform (DFT). The DFT of the heat and temperature signals is calculated as N1 1 X XðkÞ ¼ pffiffiffiffi xðtÞej2pkt=N N t¼0

the frequency band of interest w.r.t. the parameters hzg, hzt, and hp [29]

 Gdelay 2

ð11Þ

where cl represents a constant, and zit the zeros of the transient model. The poles are grouped in the parameter vector hp, while the constant and the zeros are grouped in hzg and hzt. Note that the poles of the transient model correspond to the poles of the heat transfer model [29]. The rational FRF model is estimated in two steps. First, the time delay is estimated as the time difference between the start of the first heat injection and the first observed fluid temperature increase. Second, the rational model that describes the relation between the noisy input and output spectra (Gdelay or Gdelay ) is 1 2 estimated, using the Fdident MATLAB Toolbox [30]. The following Sample Maximum Likelihood (SML) cost function is minimized in

The best rational models for both types are also evaluated in the time domain. These time domain models can be compared with models obtained using standard time domain methods, for example, an analytical model based on data coming from a thermal response test [32–34]. This section first explains the transformation of the two rational pffiffi FRF models (in s and s) to an impulse response in the time domain. Next, the relation between the impulse response, the applied heat input signal and the temperature response is described. Finally, the evaluation and validation of the models in the time domain is explained. 4.3.1. Inverse-transformation of rational model in Laplace variable s to an impulse response in the time domain The rational model found in the Laplace variable s (5) is inversetransformed to an impulse response in the time domain. First, the rational model is split into partial fractions. Next, every partial fraction is inverse-transformed to the time domain and summedup in order to find the total impulse response:

G1 ðs; hÞ ¼

np np X X Ai ! g 1 ðt; hÞ ¼ Ai epi t s  pi i¼1 i¼1

ð14Þ

G. Monteyne et al. / International Journal of Heat and Mass Transfer 69 (2014) 129–139

Note: The order of the numerator of the rational model has to be smaller than the order of the denominator to apply the inverse Laplace transform. This requirement is expected, consistent with the physics of the heat transfer. 4.3.2. Inverse-transformation of rational model in Warburg variable pffiffi s to an impulse response in the time domain The impulse response in the time domain, corresponding to the pffiffi rational model in the Warburg variable s, is found using similar steps as for the rational model in the Laplace variable. First, the rational model (6) is split into partial fractions. Then, every partial fraction is inverse-transformed to the time domain by use of the following relationship [35]:  np np pffiffi  X pffiffi  X A 1 2 pffiffi i ! g 2 ðt;hÞ ¼ G2 s;h ¼ Ai pffiffiffiffiffiffi  ai eai t erfc ai t ð15Þ s þ ai pt i¼1 i¼1 Note: This inverse Laplace transform implies that the order of the numerator of the rational model has to be smaller than the order of the denominator. This inverse Laplace transformation is calculated in MATLAB by use of the code developed by Weideman [36]. This code calculates 2 the Faddeeva function wðzÞ ¼ ez erfcðjzÞ [37] in a simple and efficient way. Thus, the inverse Laplace transform is found as

g 2 ðt; hÞ ¼

  np X 1 Ai pffiffiffiffiffiffi  ai wðzÞ pt i¼1

ð16Þ

g d ðk  lÞ ¼

Z

133

ðklÞT s

gðsÞds

ð22Þ

½ðklÞ1T s

It is assumed that the system is causal and that all heat fluxes and temperatures are at rest at time t = 0. The first heat injection is then at time t = 0. After these assumptions Eq. (21) can be rewritten as

T f ;out ðkT s Þ  T 0 ¼

k X g d ðlÞqMLBS ðk  lÞ d

ð23Þ

l¼1

Thus, the measured discrete temperature Tf,out(kTs) can be calculated by use of the known discrete heat input qMLBS and the estid mated discrete impulse response gd. A time delay is coupled to the temperature signal in order to take the transport time of the fluid (7) into account. k X   T delay g d ðlÞqMLBS ðk  lÞ d f ;out kT s þ t delay  T 0 ¼

ð24Þ

l¼1

Eq. (24) is then used to calculate the temperature response of the system on an arbitrary MLBS applied heat input signal. In order to avoid leakage errors when calculating the temperature response, the heat input signal is oversampled (Toversample = 3 s) w.r.t. the switching time (Tswitch = 6 min). The calculated temperature response is then compared to the measured temperature response in order to evaluate and validate the model.

with

pffiffi z ¼ jai t

ð17Þ

For jimag (z)j > 26.6, the code for calculating the Faddeeva function returns an infinite value. In this case, the Faddeeva function is approximated by its asymptotic value [36], namely,

pffiffiffiffi  wðzÞ  j= pz

ð18Þ

4.3.3. Temperature response on the heat input signal in the time domain Since the considered system is assumed to be Linear Time Invariant (LTI), the temperature signal Tf,out(t) is calculated by convolving the heat signal q(t) with the impulse response g(t) of the system, viz:

Z

T f ;out ðtÞ  T 0 ¼ gðtÞ  qðtÞ ¼

1

gðsÞqðt  sÞds

ð19Þ

1

where ⁄ represents the convolution operator, and T0 the initial homogeneous temperature. The calculated impulse responses g1(t, h) and g2(t, h) are discrete-time approximations of the continuous-time impulse response g(t). That is why Eq. (19) has to be rearranged in a discrete-time form. We assume that the heat input signal is a Maximum Length Binary Sequence (MLBS), thus, it can be written as a sum of Zero-Order Hold (ZOH) functions, as follows:

qMLBS ðtÞ ¼

1 X

qMLBS ðlÞZOHðt  lT s Þ d

ð20Þ

l¼1

where Ts represents the sampling time. Combining Eqs. (19) and (20) gives Z 1 1 X T f ;out ðkT s Þ  T 0 ¼ gðsÞ qMLBS ðlÞZOHððk  lÞT s  sÞds d 1

¼

1 X

l¼1

qMLBS ðlÞ d

l¼1

¼

1 X l¼1

qMLBS ðlÞ d

Z

1

gðsÞZOHððk  lÞT s  sÞds

5. Experiment This section describes three experiments performed at a GSHP test facility at Chalmers University of Technology, Sweden. First, Section 5.1 discusses the experimental setup. Second, Section 5.2 explains the three different experiments. 5.1. Experimental setup The test facility has nine 80 m deep boreholes but only one borehole is used for three tests of this study. The borehole used for this study has a diameter of 110 mm and a single polyethylene U-tube (DN40 PN6) is inserted in it. The borehole is not grouted but is, instead, filled with groundwater, which is within 1 m of the ground surface. The test setup includes a variable capacity electric resistance heater, a variable speed circulation pump, and temperature and flow meters (Fig. 3). The temperature measurements of the circulating fluid are taken using Pt100 sensors at two different instances: first, when entering (Tf,in) or exiting (Tf,out) the borehole; and second, at the input to and output from the electric heater. The flow rate of the fluid is also measured twice: first using an installed vortex flow meter at the input to the electric heater; and second over the balancing valve before entering the borehole. The input power (q) to the electric heater is measured by means of a high-accuracy meter. A state-of-the-art data acquisition and storage system is used for taking measurements. All the data, including input power and temperature measurement, are acquired every minute. For the range of conditions used in this study, the power and temperature measurements have accuracies within 1% and ±0.4 °C, respectively. The three experiments are performed several weeks apart to ensure that the borehole is fully recovered from the effects of the earlier test before starting the new test.

1

Z

ðklÞT s

½ðklÞ1T s

gðsÞds ¼

1 X

qMLBS ðlÞg d ðk  lÞ ð21Þ d

l¼1

where the discrete impulse response gd(l) is related to the continuous impulse response g(t) as

5.2. Description of the three experiments During the first experiment, heat is injected into the ground at a constant injection rate. This experiment will be used to demonstrate the state-of-the-art modeling. The second and third

134

G. Monteyne et al. / International Journal of Heat and Mass Transfer 69 (2014) 129–139

mated equivalent-radius rp. The grout (kb, ab), pipe (Rp) and fluid (Cp) properties are known by default since the materials are known. The obtained analytical model is used in Section 6.1 to calculate the temperature of the fluid exiting the borehole for different heat input signals of experiments 2 and 3. 5.2.2. Second and third experiment with heat input according to MLBS The two other experiments were performed similarly. The system (heat transfer) was excited by switching the electric heating on and off, consistent with a Maximum Length Binary Sequence (MLBS) (Fig. 5). This MLBS was repeated in order to measure two periods of the heat input and temperature signal. A MLBS looks like the on/off cycle of an operational heat pump. However, an MLBS corresponds to a random binary signal in the time domain, and has the advantage that its power can be concentrated in the frequency band of interest. The frequency band is chosen by taking the following criteria into account:

Fig. 3. Schematic diagram of experimental setup. 1: Borehole heat exchanger; 2: temperature sensor; 3: electric heater and power meter; 4: vortex flow meter; 5: circulating pump; 6: balancing valve (used for differential pressure and flow measurements); Tf,in: temperature of the fluid entering the borehole; Tf,out: temperature of the fluid exiting the borehole.

experiments apply a MLBS heat input signal. The second experiment is used to demonstrate the proposed modeling method. The third experiment is used for validation. The second and third experiments are also used for comparing the performance of the state-of-the-art analytical and the proposed frequency domain methods.

5.2.1. First experiment with constant heat input A thermal response test was conducted in accordance with [38] recommendations. The test was carried out for 65 h with an input power of 5.55 ± 0.09 kW. Flow regime was kept turbulent throughout the test. The measured mean fluid temperature is shown in Fig. 4. The test was evaluated using line-source approximation approach [33]. Estimations of ground thermal conductivity (W (m K)1) and borehole thermal resistance (m K W1) were obtained. The estimated values are then used as inputs to the analytical model described in Section 3. Ground thermal conductivity estimation k is a direct input to the analytical method whereas the borehole resistance estimation is implicitly added to the esti-

20

°

Temperature ( C)

25

15

10

5

0

20

40

60

Time (h) Fig. 4. Measured mean fluid temperature during the first experiment.

1. The highest required frequency is determined by the speed of the heat transfer process. Since the borehole heat transfer is a slow process (in the order of hours), high frequency (frequency > 0.28 mHz) heat variations are suppressed. Taking a safety margin on the highest frequency into account, the power of the heat input signal is chosen to be concentrated at frequencies smaller than 1 mHz. 2. The smallest frequency is determined by the limited measurement time. A typical thermal response test takes 50 h [32]. It has been shown that longer test times (up to 100 h) lead to more accurate results [34]. However, the total test time should remain as small as possible in practice. In order not to exceed the test time of 100 h, the length of one period of the MLBS was chosen to approximate 50 h. The power (2.75 kW) of the heat input signal was chosen to be concentrated in the frequency band between 5.4  103 mHz and 1.4 mHz (Fig. 6). This was realized by choosing a switching frequency of the MLBS (Tswitch) of 6 min and a period length of 50 h. Consequently, the chosen MLBS comprised 511 binary (on/off) signals. This makes that the chosen MLBS consists out of 511 binary (on/off) signals. The dips in the power spectral density (Fig. 6) correspond to the switching frequency of the MLBS (fswitch = 1/ Tswitch = 2.8 mHz). The measured temperature of the fluid exiting the borehole is shown for the two experiments in Fig. 7. The temperature increases slowly to a steady state value. This temperature increase is proportional to the mean value of the injected heat. The fast temperature variations are caused by the on–off switching of the heater. 6. Results and discussion This section discusses and compares the measured temperature with the temperature calculated using the different estimated models (Sections 3 and 4). First, Section 6.1 validates the model obtained from the state-of-the-art analytical method. Second, Section 6.2 discusses and validates the models obtained from the proposed method. Finally, Sections 6.3 and 6.4 compare the fluid temperatures obtained from the state-of-the-art method and the proposed method for the short- and long-term, respectively. Since the temperature response is measured every minute (Ts = 1 min), the calculated temperature response is compared with the measured temperature response at the measurement instances (every minute). The temperature error is calculated as meas T error f ;out ðkT s Þ ¼ T f ;out ðkT s Þ  T f ;out ðkT s Þ

T meas f ;out ðkT s Þ

ð25Þ

where represents the measured temperature of the fluid exiting the borehole at the kth minute.

135

Electric Heater Input

G. Monteyne et al. / International Journal of Heat and Mass Transfer 69 (2014) 129–139

on

off 0

5

10

15

20

25

Time (h) Fig. 5. Part of the input signal for the electric heater (MLBS).

Amplitude Input (dB)

60 band 40

of interest

20 0 −20 0

2

4

6

Frequency (Hz)

8 −3

x 10

Fig. 6. Power spectral density of the input signal for the electric heater (MLBS).

in the BHE. As the analytical model does not account for axial heat transfer in the borehole, a time step shorter than the residence time of the fluid in the BHE will lead to an error in the model output. The absolute difference between the measured and calculated temperature is smaller than 2 °C (Fig. 8: second experiment; Fig. 9: third experiment). The standard deviation of the temperature difference corresponds to 0.42 °C and 0.41 °C for the second and third experiments, respectively. The errors are rather high at the beginning of the experiments. This is because, at the beginning of the experiments, the temperature increase is fast. For the rest of the experiment, the temperature error contains some fast variations. This indicates that some of the faster dynamics are not perfectly simulated due to the inability of the model to handle time steps shorter than the residence time of the fluid in the BHE. 6.2. Temperature errors: proposed method

In this section the analytical model (Section 3) is used to calculate the fluid temperature variation when the heating is switched on and off consistent with a MLBS (second and third experiment). The analytical model uses thermal conductivity and borehole resistance values obtained from the thermal response test (first experiment in Section 5). For both experiments 2 and 3, the chosen switching time (6 min) of the MLBS signal is in the range of the residence time of the fluid in the borehole. The sampling time (1 min) is even shorter than the residence time of the fluid in the borehole. One of the reasons behind choosing a fast sampling and switching time is to highlight limitations of the state-of-the-art model when simulating time scales that are smaller than the residence time of the fluid in the BHE. In reality, the influence of any change in the heat injection rate is delayed and dampened by the fluid residence time

This section discusses the temperature response of the BHE obtained via the proposed method (Section 4). The data obtained from the second experiment consists of two periods. Since the proposed method only requires one period to obtain a rational heat transfer model, the first measured period of the second experiment is selected. This means that only 51.1 h of measurement data are used instead of the original 102.2 h. The second period is then used to test the extrapolation capabilities of the obtained heat transfer model. The rational model is estimated in the low frequency band from 0 Hz to 0.28 mHz in order to avoid the influence of the resonance at 2.8 mHz due to the round-trip time of the water. Avoiding this resonance has the advantage that the rational FRF model becomes simpler (lower model order). However, the rational FRF model will not contain this resonance information any longer. This is not a

Fluid temperature ( °C) experiment 3

Fluid temperature ( °C) experiment 2

6.1. Temperature errors: State-of-the-art method

20 16 12

10

20

30

40

10

20

30

40

50

60

70

80

90

100

50 60 Time (h)

70

80

90

100

20 16 12

Fig. 7. Measured temperature of the fluid exiting the borehole during the second and third experiment.

G. Monteyne et al. / International Journal of Heat and Mass Transfer 69 (2014) 129–139

°

Temperature error ( C)

136

0 −1

20

40

60

80

100

Time (h)

°

Temperature error ( C)

Fig. 8. Difference between the measured and calculated (using state-of-the-art method) temperature of the fluid exiting the borehole (25) for the second experiment.

0 −1

20

40

60

80

100

Time (h) Fig. 9. Difference between the measured and calculated (using state-of-the-art method) temperature of the fluid exiting the borehole (25) for the third experiment.

problem since the main goal is to model the low frequency behavior. Thus, the advantage of a simpler model outweighs the exclusion of the resonance information. The model order selection procedure (Section 4) leads to a rational s model of order 2/3 and pffiffi a s model of order 1/4, viz:

Gdelay s; ^h ¼ e300s 1 pffiffi

s; ^h ¼ e300s Gdelay 2

8:86  103 s2 þ 24:5s þ 0:0025 5:65  109 s3

þ 6:50  107 s2

 4

C

þ 3:44  10 s þ 1 W

pffiffi models discussed in this paper, the s model estimated from the proposed method, using 51.1 h of data, is best at predicting the pffiffi temperature variation over 102.2 h. As a consequence, the s model obtained from the proposed method is selected as the best model for the heat transfer between the injected heat and the temperature of the fluid.

ð26Þ

pffiffi  6:03 s þ 0:0059 C pffiffi 6 2 6 pffiffi 3:245  10 s þ 4:69  10 s s þ 3:61  105 s þ 2:97  103 s þ 1 W

ð27Þ

The impulse responses that correspond to the best rational models (Section 4) are used to calculate the temperature variations caused by the injected heat in the time domain (Section 4.3). Eq. (23) is used to calculate the temperature variation of the fluid caused by the applied heat input signal. The calculated temperature response is then compared with the measured temperature response in order to evaluate the model. The temperature errors (25) between the measured and calculated fluid temperatures obtained from the proposed method are pffiffi shown in Fig. 10 (s model) and Fig. 11 ( s model). Recall that only one measurement period (of 51.1 h) was used for the estimation of the models by the proposed method. However, the model has been used to calculate the fluid temperature for two periods in order to check its extrapolation capability. The temperature error obtained from the best rational model in the s variable (Fig. 10) never exceeds 1.6 °C and has a standard deviation of 0.36 °C, but shows a strong variation. The temperature error during the first measured period indicates that the transient was not well modeled. During the second period this transient becomes less dominant and thus the calculated temperature is more accurate. pffiffi The s model (Fig. 11) delivers a smaller temperature error than the s model. The temperature error never exceeds 0.8 °C and has a standard deviation of 0.20 °C. The temperature is slightly underestimated at the start due to the fact that the transient is not completely captured by the proposed model. After one period the calculated temperature corresponds to the measured temperature. During the second period the temperature difference increases. However, this increase is expected to converge. Out of the three

6.2.1. Validation of the best heat transfer model pffiffi In this section the rational s model obtained via the proposed method (Section 4) is validated using the data from the third experiment (Section 5). The heat input signal is combined with the impulse response of this rational model in order to calculate the temperature response. The absolute difference between the calculated temperature response and the measured one never exceeds 0.65 °C and has a pffiffi standard deviation of 0.17 °C. This result validates the rational s model obtained from the proposed method. 6.3. Comparison of state-of-the-art and proposed method for shortterm heat transfer pffiffi The rational s model obtained with the proposed method results in a smaller temperature error (Fig. 12) than the state-ofthe-art analytical model (Fig. 9) when the heat is switched on and off consistent with a MLBS. The temperature error obtained with the state-of-the-art method contains fast variations due to its inability to handle time-steps shorter than residence time of the fluid in the BHE. These faster variations are better modeled with the proposed method which results in a smaller temperature error.

6.4. Comparison of state-of-the-art and proposed method for longterm heat transfer The most suitable application of rational FRF model is to provide a short-term heat transfer model for a borehole, which can be used to develop, for example, a controller for a GSHP system. However, the model can also predict the long-term response of a BHE with reasonable accuracy.

137

Temperature error ( °C)

G. Monteyne et al. / International Journal of Heat and Mass Transfer 69 (2014) 129–139

0.8 0 −0.8 20

40

60

80

100

Time (h)

Temperature error ( °C)

Fig. 10. Difference between the measured and calculated temperature of the fluid exiting the borehole (25) (second experiment). The fluid temperature was calculated by use of the s model with order 2/3.

0.8 0 −0.8 20

40

60

80

100

Time (h)

Temperature error ( °C)

Fig. 11. Difference between the measured and calculated temperature of the fluid exiting the borehole (25) (second experiment). The fluid temperature was calculated by use pffiffi of the s model with order 1/4.

0.8 0 −0.8 20

40

60

80

100

Time (h) pffiffi s model obtained via the proposed method by use of the temperature difference between calculated and measured fluid temperature of

Heating load (kW)

Fig. 12. Validation of the rational the third experiment.

4 2

50

100

150 200 Time (days)

250

300

350

300

350

Fluid temperature ( °C)

Fig. 13. Heat extraction profile during one year.

10 7.5 5 2.5 0

1

50

100

150

2

3

200

4

250

Time (days) pffiffi Fig. 14. Fluid temperature variations simulated by the state-of-the-art analytical method (gray dashed line) and the FRF method ( s model with order 1/4, black solid line) during year 10 of the periodic heat extraction. The fast temperature variations during the first five days are enlarged.

G. Monteyne et al. / International Journal of Heat and Mass Transfer 69 (2014) 129–139

Fluid temperature difference ( °C)

138

0.5 0

1

2

3

4

−0.5 50

100

150 200 Time (days)

250

300

350

pffiffi Fig. 15. Difference between the temperature variations simulated by the state-of-the-art analytical method and the FRF method ( s model with order 1/4) during year 10. The difference between the fast temperature variations during the first five days is enlarged.

In this section the long-term performance of the FRF model is investigated using a simulation study. A typical heat extraction profile for a Swedish single family house [39] is considered with annual hourly loads (Fig. 13). This yearly profile is then repeated over 10 years. The case with only heat extraction is chosen intentionally to avoid recovery of ground temperature in the heat injection mode. Fluid temperatures are simulated using the FRF model and the analytical approach suggested by Claesson and Javed [40] for the 10th year of operation (Fig. 14). The analytical approach uses analytical radial solution (discussed in Section 3) for shorter times and a finite line source solution for longer times. The FRF model can predict the fluid temperature variation for the 10th year of operation with an acceptable accuracy (Figs. 14 and 15). The difference between the simulated fluid temperatures remains small, more precisely, it never exceeds 0.7 °C. It should, however, be stressed that the main goal of this rational FRF model is to help with the design of a controller that optimizes the performance of a GSHP system. Thus, a perfect simulation of the longterm heat transfer is not required. The long-term heat transfer model may be improved by combining the rational FRF model with a thermal response function, based on the finite line source solution. The impulse response that was obtained from the rational FRF model corresponds to the derivative of a traditional step-response function (e.g. g-function of [15]). In other words, the integral of the impulse response corresponds to the thermal response function. This means that the FRF model can be reformulated as a thermal response function for the short-term heat transfer. This short-term thermal response function can then be combined with a thermal response for the long-term heat transfer. This is, however, beyond the scope of this article, since the focus here is on obtaining a short-term heat transfer model for designing a controller for GSHP systems. 7. Conclusions This paper proposed a new frequency domain method to model the heat transfer in a borehole heat exchanger. The proposed method obtained a heat transfer model from in situ measurements. The measurements were taken using a thermal response test performed with multiple short-duration heat pulses. The measurement data obtained from the innovatively performed thermal response test contained more information about the short-term heat transfer than a traditional test. The proposed method was compared to a state-of-the-art analytical model. The models were compared by evaluating their ability to predict the fluid temperature variations caused by different heat injection and extraction rates. The comparison was based on experimental data. The proposed method delivered a more accurate model of the short-term heat transfer in the borehole than the analytical method. The higher quality of short-term heat transfer model obtained from the proposed method makes it better suited for efficient controller design.

Acknowledgments The Belgian part of this work is sponsored by the Vrije Universiteit Brussel, dept. ELEC, Brussels, Belgium, Fund for Scientific Research (FWO-Vlaanderen), the Flemish Government (Methusalem), and the Belgian Federal Government (IUAP VI/4). The Swedish part of this work is sponsored by Chalmers University of Technology, Gothenburg, Sweden and the Swedish Energy Agency through their Program EFFSYS+.

References [1] Q. Gao, M. Li, M. Yu, Experiment and simulation of temperature characteristics of intermittently-controlled ground heat exchanges, Renewable Energy 35 (6) (2010) 1169–1174. [2] T. Salque, D. Marchio, P. Riederer, Neural predictive control for single-speed ground source heat pumps connected to a floor heating system for typical french dwelling, Build. Serv. Eng. Res. Technol. (2013), http://dx.doi.org/ 10.1177/0143624413480370. . [3] W. Gang, J. Wang, Predictive ANN models of ground heat exchanger for the control of hybrid ground source heat pump systems, Appl. Energy (2013). . [4] C. Verhelst, Model predictive control of ground coupled heat pump systems in office buildings, Ph.D. Thesis, Katholieke Universiteit Leuven, 2012. [5] M. Beck, M. de Paly, J. Hecht-Méndez, P. Bayer, A. Zell, Evaluation of the performance of evolutionary algorithms for optimization of low-enthalpy geothermal heating plants, in: Proceedings of the 14th International Conference on Genetic and Evolutionary Computation Conference, GECCO ’12, ACM, New York, NY, USA, 2012, pp. 1047–1054, http://dx.doi.org/10.1145/ 2330163.2330309. [6] M. de Paly, J. Hecht-Méndez, M. Beck, P. Blum, A. Zell, P. Bayer, Optimization of energy extraction for closed shallow geothermal systems using linear programming, Geothermics 43 (2012) 57–65. [7] L.R. Ingersoll, O.J. Zobel, A.C. Ingersoll, Heat Conduction with Engineering, Geological, and Other Applications, McGraw-Hill, New York, 1954. [8] F. De Ridder, M. Diehl, G. Mulder, J. Desmedt, J. Van Bael, An optimal control algorithm for borehole thermal energy storage systems, Energy Build. 43 (10) (2011) 2918–2925. [9] S. Javed, J. Claesson, New analytical and numerical solutions for the short-term analysis of vertical ground heat exchangers, ASHRAE Trans. 117 (1) (2011) 3– 12. [10] P. Eskilson, J. Claesson, Simulation model for thermally interacting heat extraction boreholes, Numer. Heat Transfer 13 (2) (1988) 149–165. [11] C. Yavuzturk, J.D. Spitler, A short time step response factor model for vertical ground loop heat exchangers, ASHRAE Trans. 105 (2) (1999) 475–485. [12] R. Al-Khoury, P.G. Bonnier, R.B.J. Brinkgreve, Efficient finite element formulation for geothermal heating systems. Part I: Steady state, Int. J. Numer. Methods Eng. 63 (7) (2005) 988–1013. [13] R. Al-Khoury, P.G. Bonnier, Efficient finite element formulation for geothermal heating systems. Part II: Transient, Int. J. Numer. Methods Eng. 67 (5) (2006) 725–745. [14] M. He, S. Rees, L. Shao, Simulation of a domestic ground source heat pump system using a three-dimensional numerical borehole heat exchanger model, J. Build. Perform. Simul. 4 (2) (2011) 141–155. [15] P. Eskilson, Thermal analysis of heat extraction boreholes, Ph.D. Thesis, University of Lund, Department of Mathematical Physics, Lund, Sweden, 1987. [16] C. Yavuzturk, J. Spitler, S. Rees, A transient two-dimensional finite volume model for the simulation of vertical u-tube ground heat exchangers, ASHRAE Trans. 105 (2) (1999) 465–474. [17] H. Carslaw, J. Jaeger, Conduction of Heat in Solids, vol. 1, Clarendon Press, Oxford, 1947. [18] L. Lamarche, B. Beauchamp, A new contribution to the finite line-source model for geothermal boreholes, Energy Build. 39 (2) (2007) 188–198. [19] L. Lamarche, B. Beauchamp, New solutions for the short-time analysis of geothermal vertical boreholes, Int. J. Heat Mass Transfer 50 (7-8) (2007) 1408– 1419.

G. Monteyne et al. / International Journal of Heat and Mass Transfer 69 (2014) 129–139 [20] G. Bandyopadhyay, W. Gosnold, M. Mann, Analytical and semi-analytical solutions for short-time transient response of ground heat exchangers, Energy Build. 40 (10) (2008) 1816–1824. [21] J. Claesson, S. Javed, An analytical method to calculate borehole fluid temperatures for time-scales from minutes to decades, ASHRAE Trans. 117 (2) (2011) 279–288. [22] M. De Carli, M. Tonon, A. Zarrella, R. Zecchin, A computational capacity resistance model (CaRM) for vertical ground-coupled heat exchangers, Renewable Energy 35 (7) (2010) 1537–1550. [23] A. Zarrella, M. Scarpa, M.D. Carli, Short time step analysis of vertical groundcoupled heat exchangers: the approach of CaRM, Renewable Energy 36 (9) (2011) 2357–2367. [24] D. Bauer, W. Heidemann, H. Müller-Steinhagen, H.-J.G. Diersch, Thermal resistance and capacity models for borehole heat exchangers, Int. J. Energy Res. 35 (4) (2011) 312–320. [25] P. Pasquier, D. Marcotte, Short-term simulation of ground heat exchanger with an improved TRCM, Renewable Energy 46 (2012) 92–99. [26] J. Allwright, G. Papavasiliou, On linear programming and robust modelpredictive control using impulse-responses, Syst. Control Lett. 18 (2) (1992) 159–164. [27] R. Pintelon, G. Vandersteen, J. Schoukens, Y. Rolain, Improved (non)parametric identification of dynamic systems excited by periodic signals – the multivariate case, Mech. Syst. Signal Process. 25 (8) (2011) 2892– 2922. [28] R. Pintelon, J. Schoukens, G. Vandersteen, Frequency domain system identification using arbitrary signals, IEEE Trans. Automat. Control 42 (1997) 1717–1720.

139

[29] R. Pintelon, J. Schoukens, System Identification, second ed., A Frequency Domain Approach, Wiley-IEEE Press, Hoboken, USA, 2012. [30] I. Kollár, Frequency Domain System Identificaton Toolbox (FDIDENT) for Matlab, Gamax Ltd, Budapest, 2004–2009. [31] R. Shibata, Asymptotically efficient selection of the order of the model for estimating parameters of a linear process, Ann. Stat. (1980) 147–164. [32] P. Mogensen, Fluid to duct wall heat transfer in duct system heat storage, in: Proc. Int. Conf. on Subsurface Heat Storage in Theory and Practice, Stockholm, Sweden, 1983, pp. 652–657. [33] S. Gehlin, Thermal response test: method development and evaluation, Ph.D. Thesis, LuleåUniversity of Technology, Sweden, 2002. [34] S. Javed, Thermal modelling and evaluation of borehole heat transfer, Ph.D. Thesis, Chalmers University of Technology, Sweden, 2012. [35] I. Petráš, Fractional-Order Nonlinear Systems: Modeling, Analysis and Simulation, Springer Verlag, 2011. [36] J. Weideman, Computation of the complex error function, SIAM J. Numer. Anal. (1994) 1497–1518. [37] V. Faddeeva, N.M. Tables of Values of the Function pffiffi Terent’ev, Rz 2 2 wðzÞ ¼ ez ð1 þ 2i= ðpÞ 0 et dtÞ for Complex Argument, Pergamon Press, 1961. [38] ASHRAE, HVAC Applications, American Society of Heating, Atlanta, USA, 2007. [39] J. Spitler, L. Xing, J. Cullin, D. Fisher, J. Shonder, P. Im, Residential ground source heat pump systems utilizing foundation heat exchangers, in: Proc. of Clima 2010, Antalya, Turkey, 2010. [40] J. Claesson, S. Javed, A load-aggregation method to calculate extraction temperatures of borehole heat exchangers, ASHRAE Trans. 118 (2012) 530– 539.