Effect of depth and fluid flow rate on estimate for borehole thermal resistance of single U-pipe borehole heat exchanger

Journal Pre-proof Effect of depth and fluid flow rate on estimate for borehole thermal resistance of single U-pipe borehole heat exchanger

Changxing Zhang, Xinjie Wang, Pengkun Sun, Xiangqiang Kong, Shicai Sun PII:

S0960-1481(19)31527-7

DOI:

https://doi.org/10.1016/j.renene.2019.10.036

Reference:

RENE 12405

To appear in:

Renewable Energy

Received Date:

22 May 2019

Accepted Date:

10 October 2019

Please cite this article as: Changxing Zhang, Xinjie Wang, Pengkun Sun, Xiangqiang Kong, Shicai Sun, Effect of depth and fluid flow rate on estimate for borehole thermal resistance of single U-pipe borehole heat exchanger, Renewable Energy (2019), https://doi.org/10.1016/j.renene.2019.10.036

This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2019 Published by Elsevier.

Journal Pre-proof

Effect of depth and fluid flow rate on estimate for borehole thermal resistance of single U-pipe borehole heat exchanger Changxing Zhanga,*, Xinjie Wanga, Pengkun Suna, Xiangqiang Kongb, Shicai Suna aShandong

Key Laboratory of Civil Engineering Disaster Prevention and Mitigation, Shandong

University of Science and Technology, Qingdao 266590, PR China bCollege

of Mechanical and Electronic Engineering, Shandong University of Science and

Technology, Qingdao 266590, PR China Abstract: Accurate estimates for ground thermal parameters and borehole thermal resistance are important to improve the design of borehole heat exchangers (BHEs) in ground-coupled heat pump systems(GCHPs). In order to improve the estimating accuracy of borehole thermal resistance, this paper presents a simple analytical method for evaluating the actual averaged –over-the –depth mean fluid temperature (MFT) in the U-pipe of BHE to calculate borehole thermal resistance Rb . Furthermore, the effects of borehole depth and volumetric flow rate on the calculating RMSE distribution between borehole thermal resistance Rb and effective borehole thermal resistance Rb* are investigated. The conclusion shows that the relative deviation between the two borehole thermal resistances corresponding to the volumetric flow rate 1.5e-4m3/s increases from 4.2% to 29.7% when borehole depth changes from 50m to 200m. Finally, the borehole depths corresponding to different volumetric flow rate are optimized to find the boundary line where Rb are nearly equal to

Rb* in the operating time, and the impacts of grout thermal conductivity and heat rate per unit depth of BHE on the boundary line are quantitatively analyzed. Volumetric flow rate has more effect on

Rb with the higher grout thermal conductivity, the relative error between the two Rb corresponding to V=1.5e-4m3/s and V=3e-4m3/s is 10.8% for kg =2.3W/(m.℃). The effect of heat rate per unit depth of BHE on Rb is very limited, the relative error between the two Rb corresponding to ql  50 W/m and ql  80 W/m is 5.4% under V=1.5e-4m3/s, and it is only 3.7% under V=3e-4m3/s. Keywords: Vertical temperature distribution; Mean fluid temperature; Borehole thermal resistance; Borehole depth; Volumetric flow rate *Corresponding author at: QianWanGang Road 579#, QingDao, ShanDong, China, 266590. Tel./Fax: +86 532 86057593. E-mail address: [email protected] (C.X. Zhang).

1

Journal Pre-proof 1. Introduction With ground as a heat source/sink, ground coupled heat pump system (GCHPs) plays an important role in reducing carbon emissions associated with building heating and cooling, and borehole heat exchangers (BHEs) are the critical component for improving energy efficiency and decreasing environmental impact. The performance of BHEs strongly depends on the ground thermal properties and borehole thermal resistance which is widely used to analyze the heat transfer between the fluid in the U-pipe and the borehole wall. If borehole thermal resistance can not be evaluated accurately, an overprediction causes an overestimation of the required BHE size in the design of GCHPs, which results in a noticeable economic impact. On the contrary, an underprediction of borehole thermal resistance leads to worse operating performance of GCHPs, which even affects the reliability of GCHPs in the practical application. With the development of measurement and analysis technologies, in-situ thermal response test (TRT) has been recognized as the effective tool to obtain accurate values of ground thermal properties and borehole thermal resistanceError! Reference source not found.Error! Reference source not found., and the international standard for TRT was also published in 2015 as EN ISO 17628[3] in which the frame work of TRT was presented. For the in-situ TRT, the inlet and outlet fluid temperatures of BHE are usually recorded by the data logger with the setting sampling rate, and the arithmetic average fluid temperature (AFT) is used to calculate the effective borehole thermal resistance in the analytical model of BHE where the heat transfer rate and borehole wall temperature are assumed uniform along the length, that is, Rb* 

Tav  Tb (Tout  Tin ) / 2  Tb  ql ql

(1)

In fact, the borehole thermal resistance should be defined based on the actual averaged –over-the – depth mean fluid temperature(MFT), which can be calculated by H

H

T  T (  (Tfu  Tb )dz  0 (Tfd  Tb )dz ) / (2 H )  Tb Rb  m b  0 ql ql

2

(2)

Journal Pre-proof The computations in the above two equations are important for proper design of efficient GCHPs. Underestimating the borehole thermal resistance leads to higher system energy consumption and reduced system capacity, on the contrary, overestimating the borehole thermal resistance results in increased BHE size and cost. By comparing Eq.(1) and Eq.(2), the relations between the two thermal resistances were described by Hellstrom [4]under the boundary conditions(uniform heat transfer rate and borehole wall temperature), it is observed that how to do with the mean fluid temperature is crucial for calculating borehole thermal resistance and achieving the optimal lengths of BHEs by use of the simulation-based design tool. In the Eq.(2), the distribution of temperatures along the two legs of BHE is considered, which explain what causes the difference between two thermal resistances. Furthermore, it is short-circuiting heat transfer between the upward-flowing and downward-flowing legs of BHE makes AFT deviate from MFT. Based on a full 3D numerical model, Marcotte and Pasquier[5] proposed the p-linear average value of inlet and outlet fluid temperature with p  1 , which closely fits MFT calculated by the 3D numerical model of BHE. Moreover, the conclusion showed that AFT led to an overprediction of borehole thermal resistance in the TRT, and resulted in an overestimation of 17.5% for the required BHE size in the design. Beier[6] developed an analytical model of the actual vertical temperature profile in the ground loop for the late-time period of the in-situ TRT. The soil thermal conductivity and borehole thermal resistance can be estimated without the mean temperature approximation. Based on the analytical model, the error in the total thermal resistance was less than 5% if the volume flow rate was maintained above a threshold, which increased linearly with the total length of the borehole. The analysis also indicated a p-linear average method gives smaller errors than AFT. Du and Chen [7]presented a p-linear dimensionless average fluid temperature to estimate borehole thermal resistance of BHE, and it is compared with theoretical dimensionless fluid temperature calculated by quasi-three-dimensional models for both single and double U-pipes[8]. Results showed that the dimensionless logarithmic mean temperature for p  0 and the dimensionless geometric mean temperature for p   1 should respectively be adopted to 2

3

Journal Pre-proof reasonably estimate the thermal resistance of single and double U-pipes boreholes. Considering the variation of the p value with time, Zhang et.al [9]used p(t)-linear average method for evaluating the ground thermal parameters as well as the p values at the different sampling times. It is concluded that the proposed method led to a 6.31% reduction of the borehole thermal resistance when compared to the AFT, and the maximum relative error for the p(t)-linear average method is less than 3% for all the typical practical cases. Spitler et. al [10]proposed a tool for computing Rb* and

Rb of single U-pipe BHE using Excel/VBA. Sample results and validation against experimental measurements from several boreholes were presented. The effect of thermal short-circuiting is proportional to the ratio of the pipe-to-pipe conductance to the thermal capacitance of the working fluid and is more important as borehole depth significantly exceeds 100 m. Javed and Spitler[11] reviewed ten different published methods for calculating borehole thermal resistance of single Upipe BHE, and compared 216 different cases that bracketed all or almost all real-world BHEs. The sensitivity of the borehole thermal resistance to various parameters was also discussed. Beier and Spitler [12]defined a weighted average temperature (WAT) combined with 1D radial model to account for the variations in temperature with depth. The proposed method was verified with measured data from TRTs on boreholes with single and double U-pipes, as well pipe-in-pipe (coaxial) boreholes, it gave more accurate results than AFT without requiring computationally intensive 3D models. The weighting factor f corresponding to a specific 1D model can be calculated by forcing the 1D model to agree with WAT from the quasi-3D model at the same time. That is f 

Twa  Tout Tin  Tout

(3)

The weighting factors are time-varying and account for short-circuiting at transit time. From Eq.(3), WAT correspond to AFT by setting f =0.5 for all times. Based on the transient weighting factor to calculate an average circulating fluid temperature along the borehole, Beier et. al [13]applied the proposed model to analyze multi-flow rate TRTs (MFR-TRT). The weighting factor model matched the measured inlet and outlet temperatures over a wide range of flow rates where the

4

Journal Pre-proof Reynolds number varies from 2400 to 33,300, it showed the advantage in calculating the fluid inlet and outlet temperatures at the smaller flow rates. Zanchini and Jahanbin [14]used Comsol [15]to study the effects of the temperature distribution on borehole thermal resistance of double U-pipes BHE whose inlet fluid temperature was set constant. The results of 3D simulations showed that the

Rb* was much higher than Rb . The percent difference can exceed 28% for high thermal conductivity grout and low flow rate. Higher percent differences would occur for longer BHEs with the same grout conductivity and flow rate. Furthermore, they presented the correlation approach for determining MFT of double U-pipes BHE, and provided the tables of a dimensionless coefficient that allowed an immediate evaluation of MFT in any working condition, which gave better estimations of borehole thermal resistance in TRT and more accurate evaluation of the outlet temperature in dynamic simulations of GCHPs[16]. On the basis of the simulated results of BHE using three levels and four factors (depth, shank spacing, grout thermal conductivity and volume flow rate), Zanchini and Jahanbin [17]determined the expression of the dimensionless coefficient by suitable best fitting method, which can calculate the difference between MFT and AFT in different working condition. For a specific installation of GCHPs, the knowledge of ground thermal properties is a prerequisite for correct design of BHE, and ground thermal conductivity and borehole thermal resistance are the essential parameters which are identified by solving inverse heat transfer problem using the experimental data from TRT. Different from estimating ground thermal conductivity, borehole thermal resistance is usually decided by the differences between mean fluid temperature in U-pipe and borehole wall temperature. Therefore, how to evaluating mean fluid temperature is critical to improve the estimating accuracy of borehole thermal resistance. Considering the effect of the vertical temperature distribution in U-pipe of BHE, this paper proposed a simple analytical method for evaluating the actual averaged –over-the –depth mean fluid temperature(MFT) in single U-pipe BHE by using variable step 4th-order and 5th-order Runge-Kutta-Felhberg algorithm (RKFA)[18] to numerically solve energy equilibrium equations. The effects of borehole depth and

5

Journal Pre-proof volumetric flow rate on estimate for borehole thermal resistance are analyzed by calculating the RMSE distribution between Rb* and Rb which are corresponding to the different mean fluid temperatures in U-pipe. Finally, the effects of grout thermal conductivity and heat rate per unit depth on the boundary line of the two borehole thermal resistances are studied and compared. 2. The simple analytical method for evaluating MFT For BHE, mathematical model is crucial for evaluating MFT in the U-pipe. Compared with numerical model, the analytical models show the advantages that they are fast and easy to use in the estimation on ground thermal parameters of TRT and operating simulation of GCHPs. Cylindrical source model (CSM) extends a line source to a cylindrical source with a constant radius and compares favourably well with 3D numerical models, so it is employed in the analysis on heat transfer in the outside part of BHE. As for the inside part, quasi-three-dimensional model is welcomed because of its unique advantages in the process of deducing vertical fluid temperature distribution along the borehole. 2.1 Cylindrical source model In the cylindrical source model(CSM), the borehole is assumed as an infinite cylinder surrounded by homogeneous medium with constant properties, i.e. the ground. It also assumes that the heat transfer between the borehole and soil with perfect contact is of pure heat conduction. CSM simulates heat conduction in a soil mass subjected to a constant heat flow rate. The difference, however, is that the contact area between BHE and the soil is along the surface area of the borehole, i.e. at r  rb , the borehole radius. Fig.1 shows the single U-pipe BHE diagram. The governing partial differential equation is described in the cylindrical coordinate system, as  1 T  2T 1 T  2  ; rb  r      t r  r  r  T (r , 0)  Ts ; T (r  , t )  T ; t  0 s   q (t ) T  l   ks r r  rb 2 r  6

(4)

Journal Pre-proof The analytical solution of CSM with a constant heat rate is described as[19]: 2

q 1  e  x Fo  1  )Y ( x)  J ( x)Y (rx  )]dx T (r , Fo)  T0  l 2  2 2 [ J o (rx 1 1 o ks  0 x [ J1 ( x)  Y12 ( x)]  

(5)

G ( r , Fo )

Where： r  r / rb ; Fo 

ks t ; J and Y are the Bessel functions of the first and second kind. rb2 Cs

In Eq. (5), the expression G (r , Fo) is a function of time and distance from the borehole center, the temperature at the borehole wall, Tb , can be obtained by setting r  1 . Philippe et al. [20]compared the simulated result from CSM with the numerical results with the finite different solution of CSM given by Lee and Lam[21], and it showed that there is a quite good matching between the two results. CSM is also widely applied in the performance simulation of GCHPs and the parameter estimation in TRT[22]. soil grout material b

borehole wall

1

2

HDPE pipe

a) b) Fig.1 Single U-pipe BHE diagram: a) Cross-section; b) Isometric view

2.2 The quasi-three-dimensional model for BHE For BHE, it is necessary for evaluating MFT to accurately predict the vertical fluid temperature distribution in the U-pipe. Bauer et.al [23]presented the 3D numerical model for BHE which included the thermal capacities (the fluid inside the tubes/the grouting material), the fluid temperature distribution along the borehole depth was calculated at different times. Using a multiblock mesh to represent each component of BHE in three-dimensions and applying a finite volume numerical method, Rees and He [24]developed 3D BHE model for studying conduction and fluid

7

Journal Pre-proof circulation processes over both short and long timescales. Ozudogru et.al [25]proposed 3D numerical model uses 1D linear elements for simulating the flow and heat transfer inside the pipes, and the validation of the model was carried out by comparing the numerical results with the results obtained from the analytical model. Furthermore, the method to estimate the steady state thermal resistances in the borehole/energy pile was presented in order to calculate the fluid temperatures analytically. Compared the computationally time intensive numerical model, the quasi-threedimensional model for BHE is an efficient analytical approach for calculating the vertical fluid temperature distribution in the U-pipe along the BHE depth, and the accuracy has been verified by the related experiment and BHE model[25][26]. In the quasi-three-dimensional model, the fluid temperature and pipe wall temperature vary in the axial direction, and the impact of the short-circuiting among U-pipe legs of BHE can be assessed. As seen the single U-pipe BHE in Fig.1, pipe’1’ and pipe’2’ are downward and upward fluid channel respectively, the thermal resistance between pipes in the borehole can be expressed as follows:  1  R11  2 kg    1  R12  2 kg 

 r ln  b   rpo

 kg  ks  rb2  D 2   ln(rpo / rpi ) 1 ln       2 2 kp 2 rpi hf    kg  ks  rb

(6)

  rb  kg  ks  rb2  D 2   ln  ln    2   2 D  kg  ks  rb  

Generally, the flow in the U-pipe is in the region of thermally developed flow, so hf in Eq.(6) is determined by the Dittus-Boelter correlation: hf  0.0115 Re0.8 Pr n kf / rpi , n = 0.4 for heating mode(wall hotter than the bulk fluid) and 0.33 for cooling mode(wall cooler than the bulk fluid)[27][28]. In order to improve the accuracy of calculating temperature profile in the U-pipe, water dynamic physical parameters including thermal conductivity, density and kinematic viscosity coefficient are applied in this study. The relationship between water dynamic physical parameters and transient temperature is regressed respectively as follows [29]: f (Tf )  (3  108  Tf 3  7  106  Tf 2  4  105  Tf  1)  1000 Tf  [12 52]

(7)

 f (Tf )  9.467215  104  (Tf  47)( 1.629311) / f (Tf ) Tf  [0 52]

(8)

8

Journal Pre-proof kf (Tf )  3.69  1011  (Tf  273.15)3  6.19  106  (Tf  273.15)2  5.11  103  (Tf  273.15)  0.369

Tf  [0 370]

(9)

Based on energy equilibrium equation, fluid temperature profile in pipes along the borehole depth can be expressed in a dimensionless form as: d 1  a1  b 2   dZ  d 2  b1  a 2  dZ  Boundary condition : 1 (0)  1, 1 (1)   2 (1);  

Where a 

R11 R12 H ; b ; 2 2 V f c R  R12 V f c R11  R122 H

temperature 1 

(10)

The parameters in Eq.(10), the dimensionless

2 11

Tf1  Tb T T z ;  2  f2 b ; and the dimensionless depth Z  ; Tin  Tb Tin  Tb H

The analytical solution of the Eq.(10) was derived by Zeng et.al using Laplace transform[8]. Yang et.al [30]proposed the updated method for solving the order differential equation groups so that BHE effectiveness can be directly determined. The complicated analytical equations and many intermediate variables are involved in the two solving process of obtaining the analytical solution of vertical temperature profile. Different from the two complicated solving process, the differential Eq.(10) is solved numerically using variable step 4th-order and 5th-order Runge-Kutta-Felhberg algorithm (RKFA) in Matlab, which makes the computation of temperature distribution in the Upipe fast and efficient. RKFA is an implicit and explicit iteration method for the numerical solution of ordinary differential equations derived from the Runge-Kutta method[18]. RKFA uses a fourth order and fifth order accurate method to provide two estimates of the next step in the solution, the difference in these results is related to the accuracy of the solution and can be used to decrease the integration step length if the error exceeds a specified accuracy. The clever feature of RKFA is that it reduces the number of evaluations of the force functions by using an integrator that uses some of the same values in both the 4th order and 5th order methods.

9

Journal Pre-proof 2.3 Evaluating the mean fluid temperature and borehole thermal resistance When the dimensionless differential Eq.(10) is solved corresponding to different time, the actual averaged–over-the–depth mean fluid temperature(MFT) and the arithmetic average fluid temperature (AFT) can be respectively described as: H

H

1

0

0

0

Tm  (  (Tf1  Tb )dz   (Tf2  Tb )dz ) / (2 H )  (Tin  Tb ) 

1 ( Z )   2 ( Z ) dZ 2

(11)

1   2 (0) (Tin  Tb )  Tb 2

Tav 

(12)

Simultaneously, the weighting factor corresponding to the same time in Eq.(3) can be conducted by setting Tm =Twa , that is:

f 



1

0

1 ( Z )   2 ( Z ) dZ   2 (0) Tb 2  1   2 (0) (Tin  Tb )(1   2 (0))

(13)

Based on the energy balance, the effective borehole thermal resistance

Rb* is calculated by

combining Eq.(1) as follows:

Rb* 

H  1   2 (0)    2V f c  1   2 (0) 

(14)

Combing Eq.(2) and Eq.(11), the borehole thermal resistance Rb is calculated as follows: 1

T  T (Tin  Tb ) 0 Rb  m b  ql

1 ( Z )   2 ( Z ) dZ  Tb 2 ql

(15)

3. Case study and discussion 3.1 Illustrative case In this paper, a case study is used to evaluate MFT based on the presented simple analytical method, and the difference between Rb* and Rb caused by different mean fluid temperature is analyzed and discussed. Heated by the constant heat rate per unit depth of BHE, a single U-pipe BHE grouted by the mixture of bentonite and sand is taken as an example for this study, and its properties are listed in Table.1. For the flow rate V=1.5e-4m3/s, the flow in the U-pipe is turbulent, 10

Journal Pre-proof MFT, AFT, Rb* and Rb are calculated every 5 minutes in the 50-hour operating time, and the water temperature profile in the single U-pipe BHE at different time is shown in Fig.2. As shown in Fig.2, the changes of downward water temperature along the depth is larger than that of upward fluid temperature, and the increasing tendencies of inlet and outlet water temperatures decrease with the increasing operating time. It is the rising ground temperature that reduces the heat transfer potential from fluid in pipe to the soil. Table 1 The physical properties of single U-pipe BHE Borehole

H

U-pipe

d　 b

(m)

(mm)

100

130

Material 　

HDPE

r　 pi

Heat input

D

rpo　

(mm) (mm) (mm)

13

16

60

Geologic

T0

ks 　

(℃)

(W/(m.℃))

15

2.6

Cs

(kJ/(m3.℃))

kp 　

(W/(m. ℃))

0.38

q　 l

t

(W/m) (h)

80

50

Heat carrier

kg

(W/(m. ℃))

2403

2.3

c

V

(J/(kg ℃))

(m3/s)

4181

1.5e-4

Fig.2 The fluid temperature distribution in U-pipe 3.2 Weighted factor As shown in Eq.(3), the difference between WFT and AFT can be characterized by the variation of weighted factor at transit time. The weighted factors subjected to three different borehole depth are calculated and compared using Eq.(13) in this case. Fig.3 shows their variation with the operating time in the steady-flux period[12]. Although the three groups of weighted factors increase

11

Journal Pre-proof with the operating time, the weighted factors are all higher than 0.5 under H=50m, which illustrates that MFT is larger than AFT and the effects of inlet fluid temperature are underestimated in the AFT expression. However, the two groups of weighted factors are lower than 0.5 under H=100m and H=150m, which illustrates that MFT is lower than AFT and the effects of inlet fluid temperature are overestimated in the AFT expression. Consequently, the two borehole thermal resistances ( Rb* and Rb ) must be deviated due to the difference between AFT and MFT.

Fig.3 The variation of weighted factors with the operating time 3.3 Analysis on borehole thermal resistance of BHE In order to investigate the effects of MFT and AFT on different borehole thermal resistances, the relative differences between Rb* and Rb are described as R  ( Rb*  Rb ) / Rb , and they are shown in Fig.4 with the dynamic MFT and AFT. It is observed that AFT is obviously higher than MFT in the whole time, which causes Rb* larger than Rb . However, the relative difference R decreases slightly with the increasing operating time, and its changing range is about 4.65%~5.1%. According to Eqs.(14) and (15), it is necessary for BHE to further study how to reduce the differences between

Rb* and Rb by adjusting borehole depth and volumetric flow rate in BHE. In this case, the root mean squared error (RMSE) between Rb* and Rb are calculated by setting various independent values of (H, V) in the analytical method, RMSE is given by the following equation: RMSE 

1 n   ( Rb,i  Rb,* i )2 n i 1

12

(16)

Journal Pre-proof The variation ranges for borehole depth and volumetric flow rate are set to [50, 200] and [1.5, 3] respectively, e.g. H  [50, 200] (in m), V  [1.5,3] (in 10-4m3/s). RMSEs are screened by a grid search, and their distribution obtained with various independent values of (H, V) is shown in Fig.5. It is obviously seen that the narrow RMSE valley in the surf is corresponding to the field where MFTs are nearly equal to AFTs in the whole time. AFTs are higher than MFTs in the right field of RMSE valley where weighted factors f are all lower than 0.5 in WAT, and Rb* is overestimated using AFT. On the contrary, Rb* is underestimated in the left field where weighted factors f are all larger than 0.5 in WAT. It is seen that the effect of volumetric flow rate on RMSE in the 50m depth is limited by comparing with that of 200m depth, R corresponding to 50m depth only increases from 4.2% to 6.5% when volumetric flow rate changes from 1.5e-4m3/s to 3e-4m3/s. As illustrated in Fig.5, borehole depth plays an important role in the calculated deviation between Rb* and Rb under the lower volumetric flow rate, R corresponding to the volumetric flow rate 1.5e-4m3/s increases from 4.2% to 29.7% when borehole depth changes from 50m to 200m. In order to reduce R , it is necessary to obtain the optimal borehole depth corresponding to different volumetric flow rate, which can be realized by parameters optimization to find the lowest line in the bottom of narrow RMSE valley of Fig.5.

Fig.4 The relative differences between Rb* and Rb with the dynamic MFT and AFT

13

Journal Pre-proof

Fig.5 RMSE distribution corresponding to various independent values of (H, V) 3.4 Optimization on borehole depth to minimize RMSE The parameter optimization is implemented by using iterative algorithm and minimizing RMSE to find the optimal borehole depth corresponding to setting volumetric flow rate. The simulated annealing algorithm (SAA) is applied to optimize borehole depth using the presented analytical method. SAA is characterized by the need for a number of function evaluations, usually about three orders-of-magnitude greater than that commonly required for a single run of unimodal algorithms. More significantly, this algorithm allows inputting lower and upper bounds for optimizing parameters, which can produce more reliable results, independent of the quality of the initial guessed values. Details of this algorithm were reported in Ref. [31]. In the optimization process of this case, the optimization constraints of borehole depth are set H  [50, 200] , and its initially guessed value is 100m. Fig. 6 shows the dynamic optimization

process of borehole depth under the volumetric flow rate of 3e-4m3/s. As the annealing temperature in SAA always keeps falling, it is seen that RMSE fluctuates with the variation of borehole depth before 35 iterations, and their minor adjustments after 50 iterations is also observed. The optimal borehole depth (H =145m) is obtained until minimum objective function (RMSE =0.002665) at the 93th iteration is attained, which indicated that SAA is feasible and efficient in the optimization of borehole depth. Furthermore, the lowest line in the bottom of narrow RMSE valley corresponding to the smallest RMSEs under different volumetric flow rates can be found by optimizing the borehole 14

Journal Pre-proof depths. Fig.7 shows the lowest line in the bottom of narrow RMSE valley in this case. It is observed that the RMSE is lower than 0.003 except for that of V=1.5e-4m3/s. The optimal lowest line is also a boundary line where MFTs are nearly equal to AFTs in the operating time, the deviation degree between the two temperatures can be characterized using the relative position between the boundary line and the combined working point of borehole depth and volumetric flow rate. Furthermore, RMSE between Rb* and Rb can be directly obtained once the position is fixed by setting borehole depth and volumetric flow rate in the Fig.7, and the definite corrections on Rb* can be arrived in order to use Rb to improve the accuracy of sizing BHE in practical applications because ground thermal parameters and Rb* can be easily identified using the conventional TRT[1][10].

Fig.6 The dynamic optimization process of borehole depth under V=3e-4m3/s

15

Journal Pre-proof

Fig.7 RMSE distribution and the lowest line in the bottom of narrow RMSE valley 4. Influencing factors on the boundary line and corresponding borehole thermal resistance 4.1 Effect of grout thermal conductivity on the boundary line For BHE, grout thermal resistance is the important part of borehole thermal resistance, which is used to calculate heat transfer between U-pipe and borehole wall as described in Eq.(6). As the critical parameter, the effect of grout thermal conductivity on the boundary line in this paper should be investigated. Using the presented analytical method and the related parameters in the case, the boundary lines and corresponding Rb subjected to the three grout conductivities are deduced by optimizing borehole depths under different volumetric flow rates, and they are illustrated in Fig.8. Combing with Fig.5 and Fig.7 , it is seen that AFTs are higher than MFTs at same time in the upper field of the boundary line, the area expands with increasing grout conductivities. However, the slope of the boundary line for kg =1 W/(m.℃) is larger than those of the two boundary lines, which indicates the increasing tendency of optimal borehole depth for kg =1 W/(m.℃) with rising volumetric flow rate is sharper than those for the other grout thermal conductivities. It is concluded that volumetric flow rate has more effect on optimal borehole depth in the lower grout thermal conductivity. As for borehole thermal resistance Rb corresponding to the boundary line, the obvious

16

Journal Pre-proof falling tendency can be observed with the rising volumetric flow rate, and Rb for kg =3.2W/(m.℃) is the lowest. Borehole thermal resistance Rb from kg =3.2W/(m.℃) to kg =1W/(m.℃) is increased by 73.9% under V=3e-4m3/s, and the increasing ratio is 65.7% under V=1.5e-4m3/s. Undoubtedly, the effect of grout thermal conductivity on Rb corresponding to the boundary line is greater than that of volumetric flow rate.

Fig.8 The boundary lines and corresponding Rb subjected to the three grout conductivities 4.2 Effect of heat rate per unit depth of BHE on borehole thermal resistance As illustrated in Eqs.(1), (2) and (5), it is necessary to analyze the effect of heat rate per unit depth of BHE on the boundary lines and corresponding Rb . Using the presented analytical method and the related parameters in the case, the boundary lines and corresponding Rb subjected to the four heat rate per unit depth are deduced by optimizing borehole depths under different volumetric flow rates, and they are illustrated in Fig.9. Combing with Fig.5 and Fig.7, it is seen that AFTs are lower than MFTs at same time in the nether region of the boundary line, the area reduces with increasing heat rate per unit depth. The optimal borehole depth from ql  50 W/m to ql  80 W/m is decreased by 18.2% under V=1.5e-4m3/s, and the decreasing ratio is 12.3% under V=3e-4m3/s. As for borehole thermal resistance Rb corresponding to the boundary line, the obvious falling tendency can be observed with the rising volumetric flow rate, especially, the falling slope for ql  50 W/m is

17

Journal Pre-proof the sharpest in the lower range of volumetric flow rates (1.5e-4~ 2e-4m3/s). Borehole thermal resistance Rb is decreased by 5.4% under V=1.5e-4m3/s when heat rate per unit depth of BHE is adjusted from ql  50 W/m to ql  80 W/m, and the decreasing ratio is 3.7% under V=3e-4m3/s.

Fig.9 The boundary lines and corresponding Rb subjected to the four heat rate per unit depth 5. Conclusions This paper proposes a simple method for evaluating the mean fluid temperature in the U-pipe of BHE by using the variable step 4th-order and 5th-order Runge-Kutta-Felhberg algorithm[18] to solve the dimensionless energy equilibrium equations in the pipes along vertical borehole depth, which is applied to improve the estimating accuracy of borehole thermal resistance. In order to investigate the effects of borehole depth and volumetric flow rate on borehole thermal resistance, RMSE distribution between Rb* and Rb are calculated corresponding to the two independent values using MFT and AFT as the transition parameters, and the lowest line in the bottom of narrow RMSE valley under different volumetric flow rates is obtained by optimizing the borehole depth based on SAA. Moreover, the impacts of grout thermal conductivity and heat rate per unit depth of BHE on the boundary line are quantitatively analyzed in the case. The following conclusions are drawn: 1) Based on CSM and quasi-three-dimensional model, the actual averaged –over-the –depth mean fluid temperature in the pipes of BHE can be calculated using the presented simple method, which provides an efficient tool for analyzing the errors of estimated borehole thermal

18

Journal Pre-proof resistance in TRT and predicting energy consumption of GCHPs by the hourly simulations in design software for BHE fields. 2) For the effect of depth and fluid flow rate on estimate for borehole thermal resistance, it is seen that the effect of volumetric flow rate on RMSE in the 50m depth is limited by comparing with that of 200m depth, R corresponding to 50m depth only increases from 4.2% to 6.5% when volumetric flow rate changes from 1.5e-4m3/s to 3e-4m3/s. However, borehole depth plays an important role in the calculated deviation between Rb* and Rb under the lower volumetric flow rate, R corresponding to the volumetric flow rate 1.5e-4m3/s increases from 4.2% to 29.7% when borehole depth changes from 50m to 200m. 3) The boundary line where Rb are nearly equal to Rb* in the operating time can be obtained by using SAA to optimize borehole depth under different volumetric flow rate. Using the RMSE distribution as the Fig.7, Rb can be obtained conveniently after ground thermal parameters and

Rb* were identified using the conventional TRT, which provides technical support for improving the accuracy of sizing BHE in the design of GCHPs. 4) For the boundary lines subjected to different grout thermal conductivities, volumetric flow rate has more effect on optimal borehole depth in the lower grout thermal conductivity, and less effect on Rb in the lower grout thermal conductivity. R corresponding to kg =1 W/(m.℃) decreases 6.4% when volumetric flow rate increases from 1.5e-4m3/s to 3e-4m3/s, and it is 10.8% for kg =2.3W/(m.℃). 5) For the boundary lines subjected to different heat rate per unit depth of BHE, lower heat rate per unit depth has more effect on Rb in the lower range of volumetric flow rates (1.5e-4~ 2e-4m3/s). As for the effect of heat rate on Rb under different volumetric flow rate, Rb is decreased by 5.4% under V=1.5e-4m3/s when heat rate per unit depth of BHE is adjusted from ql  50 W/m to

ql  80 W/m, and the decreasing ratio is 3.7% under V=3e-4m3/s.

19

Journal Pre-proof As a future work, the presented analytical method for evaluating MFT will be validated by the distributed TRT so that it is applied in estimation on the error between Rb* and Rb to improve the accuracy of sizing BHE in the design of GCHPs. At the same time, an estimation procedure for borehole thermal resistance Rb of BHE based on this method is being developed to improve the identified accuracy of ground thermal parameter in TRT.

Acknowledgements This work was financially supported by A Project of Shandong Province Higher Educational Science and Technology Program (J16LG06), the Overseas Visiting Scholars Project of Shandong University of Science and Technology, the Natural Science Foundation of Shandong Province (NO:ZR2019MEE116), the National Natural Science Foundation of China (51776115) and SDUST Research Fund (2015KYTD104). References [1] Zhang C, Guo Z, Liu Y, Cong X, and Peng D. A review on thermal response test of ground-

coupled heat pump systems[J]. Renewable and Sustainable Energy Reviews[J], 2014, 40(12): 851-867. [2] Lamarche L, Kajl S, Beauchamp B. A review of methods to evaluate borehole thermal

resistances in geothermal heat-pump systems[J]. Geothermics, 2010, 39(2): 187-200. [3] ISO C. 17628 Geotechnical Investigation and Testing—Geothermal Testing—Determination of

Thermal Conductivity of Soil and Rock Using a Borehole Heat Exchanger[J]. CEN: Brussels, Belgium, 2015. [4] Hellström G. Ground heat storage: thermal analyses of duct storage systems[D]. Doctor thesis,

Lunds University, Sweden, 1991. [5] Marcotte D, Pasquier P. On the estimation of thermal resistance in borehole thermal

conductivity test[J]. Renewable energy, 2008, 33(11): 2407-2415. [6] Beier R A. Vertical temperature profile in ground heat exchanger during in-situ test[J]. 20

Journal Pre-proof Renewable Energy, 2011, 36(5): 1578-1587 [7] Du C, Chen Y. An average fluid temperature to estimate borehole thermal resistance of ground

heat exchanger[J]. Renewable energy, 2011, 36(6): 1880-1885. [8] Zeng H, Diao N, Fang Z. Heat transfer analysis of boreholes in vertical ground heat

exchangers[J]. International Journal of Heat and Mass Transfer, 2003, 46(23): 4467-4481. [9] Zhang L, Zhang Q, Huang G, and Du Y. A p (t)-linear average method to estimate the thermal

parameters of the borehole heat exchangers for in situ thermal response test[J]. Applied energy, 2014, 131(10): 211-221. [10] Spitler J, Saqib J, and Rachel G. Calculation tool for effective borehole thermal resistance.

Proceedings of the 12th REHVA World Congress (Clima 2016), May 22-25, Aalborg, Denmark. 2016. [11] Javed S, Spitler J. Accuracy of borehole thermal resistance calculation methods for grouted

single U-tube ground heat exchangers[J]. Applied energy, 2017, 187(2): 790-806. [12] Beier R A, Spitler J D. Weighted average of inlet and outlet temperatures in borehole heat

exchangers[J]. Applied energy, 2016, 174(7): 118-129 [13] Beier R A, Mitchell M S, Spitler J D, and Javed S. Validation of borehole heat exchanger

models against multi-flow rate thermal response tests[J]. Geothermics, 2018, 71(1): 55-68. [14] Zanchini E, Jahanbin A. Effects of the temperature distribution on the thermal resistance of

double u-tube borehole heat exchangers[J]. Geothermics, 2018, 71(1): 46-54. [15] Pepper D W, Heinrich J C. The finite element method: basic concepts and applications[M].

Taylor & Francis, 2005. [16] Zanchini E, Jahanbin A. Correlations to determine the mean fluid temperature of double U-tube

borehole heat exchangers with a typical geometry[J]. Applied energy, 2017, 206(11): 14061415. [17] Zanchini E, Jahanbin A. Simple equations to evaluate the mean fluid temperature of double-U-

tube borehole heat exchangers[J]. Applied energy, 2018, 231(12): 320-330.

21

Journal Pre-proof [18] Butcher J C. Numerical methods for ordinary differential equations[M]. John Wiley & Sons,

2016. [19] Ingersoll L R, Zabel O J, Ingersoll A C. Heat conduction with engineering, geological, and

other applications[J]. 1954. [20] Philippe M, Bernier M, Marchio D. Validity ranges of three analytical solutions to heat transfer

in the vicinity of single boreholes[J]. Geothermics, 2009, 38(4): 407-413. [21] Lee C K, Lam H N. Computer simulation of borehole ground heat exchangers for geothermal

heat pump systems[J]. Renewable Energy, 2008, 33(6): 1286-1296. [22] Yu X, Zhang Y, Deng N, Wang J, and Wang D. Thermal response test and numerical analysis

based on two models for ground-source heat pump system[J]. Energy and Buildings, 2013, 66(11): 657-666. [23] Bauer D, Heidemann W, Diersch H J G. Transient 3D analysis of borehole heat exchanger

modeling[J]. Geothermics, 2011, 40(4): 250-260. [24] Rees S J, He M. A three-dimensional numerical model of borehole heat exchanger heat transfer

and fluid flow[J]. Geothermics, 2013, 46(4): 1-13. [25] Ozudogru T Y, Olgun C G, Senol A. 3D numerical modeling of vertical geothermal heat

exchangers[J]. Geothermics, 2014, 51(7): 312-324. [26] Zhang C, Chen P, Liu Y, Sun S, and Peng D. An improved evaluation method for thermal

performance of borehole heat exchanger[J]. Renewable energy, 2015, 77(5): 142-151. [27] Florides G A, Christodoulides P, Pouloupatis P. An analysis of heat flow through a borehole

heat exchanger validated model[J]. Applied energy, 2012, 92: 523-533. [28] Rosen M A, Koohi-Fayegh S. Geothermal energy: Sustainable heating and cooling using the

ground[M]. John Wiley & Sons, 2017. [29] ISI Steam Table in SI Units, SP26, the Indian Standards Institution, 1983 [30] Yang W, Zhang S, Chen Y. A dynamic simulation method of ground coupled heat pump

system based on borehole heat exchange effectiveness[J]. Energy and Buildings, 2014, 77(7):

22

Journal Pre-proof 17-27. [31] Corana A, Marchesi M, Martini C, Ridella S. Minimizing multimodal functions of continuous

variables with the “simulated annealing” algorithm Corrigenda for this article is available here[J]. ACM Transactions on Mathematical Software (TOMS), 1987, 13(3): 262-280.

23

Journal Pre-proof Nomenclature a

constant (-)

b

constant (-)

C

volumetric heat capacity(J/( m3 ℃))

c

fluid specific heat (J/(kg ℃))

D

half spacing of U-pipe shanks (mm)

d

borehole diameter (mm)

Fo

Fourier number (-)

f

weighted factor

G

G function (-)

H

borehole depth Length (m)

h

convection heat transfer coefficient between fluid and pipe inner wall ( W/(m2.℃))

i

time variable

J

Bessel function of the first kind

k

thermal conductivity (W/(m.℃))

n

time variable

Pr

Prandtl number(-)

ql

heat transfer rate per unit depth of BHE (W/m)

R

borehole thermal resistance ( (m.℃) /W)

R*

effective borehole thermal resistance ( (m.℃) /W)

Re

Reynolds number(-)

r

radius (mm)

T

temperature (℃)

t

operating time ( h)

V

volumetric flow rate (m3/s)

x

integration variable

Y

Bessel function of the second kind

Z

dimensionless coordinate

z

axial coordinate ( m)

Greek symbols



ground thermal diffusivity (m2 /h)



coefficient of fluid kinematic viscosity ( m2/s)



dimensionless fluid temperature



fluid density ( kg/m3)

Subscripts 0

initial 24

Journal Pre-proof 1,2

pipe sequence in borehole

Av

average

b

borehole or borehole wall

f

fluid

fd

fluid down

fu

fluid upward

g

grout

in

inlet fluid

m

mean

out

outlet fluid

p

pipe

pi

pipe inner

po

pipe out

s

soil

su

single U-pipe

wa

weighted average

Abbreviations AFT

arithmetic average fluid temperature

BHE

borehole heat exchanger

CSM

cylindrical source model

GCHPs

ground-coupled heat pump system

HDPE

high density polyethylene

MFR

multi-flow rate

MFT

averaged –over-the –depth mean fluid temperature

RKFA

Runge-Kutta-Felhberg algorithm

RMSE

root mean squared error

SAA TRT WAT

simulated annealing algorithm thermal response test weighted average fluid temperature

25

Journal Pre-proof

Highlights

 

A simple analytical method is proposed to evaluate the mean fluid temperature of BHE. Effects of depth and fluid flow rate on borehole thermal resistance are investigated.



The boundary line where Rb are nearly equal to Rb is found by parameter optimization.

*