Estimation of borehole heat exchanger thermal properties

. An important issue in the design of installations with borehole heat exchangers (BHE) is knowledge of the thermo-physical properties of the grouting material and soil. The objectives of this work is to estimate the thermal conductivity and specific heat of the grouting material and the soil. The estimated values were determined on the basis of two numerical heat transport models in the borehole and using the results of the small sandbox experiment. The design of experiment technique (DOE) and the response surface methodology (RSM) were used to achieve the aim of this work. Numerical calculations were carried out with the use of the new finite element with multiple degrees of freedom (MDF) and a quasi-3D model in ANSYS package. The estimation requires minimization of four output parameters (quality criteria), and therefore it is a multiple-objectives optimization problem. Based on the DOE, RSM and multiple-objective optimization the material properties of BHE and soil was determined. For both models the discrepancies Y1 P1 are below 4 %, Y2 P1 are below 11.5%, Y1 P2 are below 6.5 % and Y2 P2 are below 11.5 %. The average discrepancies below 5 % based on the verification measurements with different operating parameters was obtained. It is noteworthy that for the two different models the same values of the estimated parameters were obtained.


Introduction
The ground source heat pump system (GSHP) is a commonly utilized renewable technology for heating and cooling of buildings [9]. The borehole thermal energy storage systems (BTES) also take on importance. To design such systems, the knowledge of the thermal properties of ground and BHE is crucial. Developments during the 1980s include design software which helped overcome under sizing problems. Current design software helps to obtain the needed inputs, especially the thermal conductivity of the ground. In paper by Mogensen [7] a method by which thermal conductivity of the ground and borehole resistance could be obtained in situ was presented. We now call it thermal response test (TRT). Further development of TRT was conduct in USA at the Oklahoma State University and in Sweden at Luleå University of Technology [9]. During TRT, one approach is to use active control of the heat input to maintain relatively uniform heat input. Alternatively, TRTs without uniform heat input may be analyzed with the parameter estimation approach [1]. In E3S Web of Conferences 46, 00017 (2018) https://doi.org/10.1051/e3sconf/20184600017 3 rd International Conference on Energy and Environmental Protection work [9] authors mentioned that developments in test rig design have included: small suitcase-sized units; telemetry; temperature measurements using optical fiber; downhole measurements using wireless sensors and wired sensors; and a low-powered test rig using a heating cable and no pumping. The analysis procedures also were developed. In the analysis procedures the groundwater flow [4] and ambient losses are included [8]. Improved analysis of tests with variable heat flow rates with [4,12,13] or without parameter estimation are conducted, improved error analyses [4] and proposals to compute the mean temperature more accurately [3]. To estimate the thermal conductivities and diffusivities of soil and grouting material, both analytical and numerical models have been used. Most of papers report that, estimation procedures based on analytical models are simpler, faster, and more robust than those based on numerical methods [5]. On the other hand, procedures using numerical models, either one-or multi-dimensional, can estimate the heat capacity of grout; and the field testing time may be short. However, they require more input parameters and are more sensitive to the accuracy of the input parameters than algorithms based on analytical models [2,6].
The objective of this work, is to conduct estimation of the thermal conductivity and specific heat of the grouting material and the soil. This research was conducted using the finite element method (FEM), DOE technique and RSM. The ANSYS software with an implemented new MDF [11] and q3D numerical model of heat and mass transport in a BHE was used. This models has been verified via experimental basis and compared to other literature BHE models in work [10]. The estimation procedure was conducted according presented in figure 1 algorithm.

Description of the sandbox experiment and computational domain
The parameters estimation were done using measurements from a small laboratory sandbox experiment. The small laboratory sandbox with dimensions of 0.46 m x 0.46 m x 2.92 m has been constructed with a borehole heat exchanger of length 2.5 m centered horizontally along the length of the box. The vertical cross section of the box is a square with sides of 0.46 m. In figure 2 the developed sandbox tank was presented. An aluminum pipe with an outer diameter of 0.080 m and 0.002 m wall thickness serves as the borehole wall. An copper pipe with an outer diameter of 0.018 m and 0.0011 m wall thickness serves as the U-tube. The distance between pipes axes is equal 0.046 m.  A testing unit to conduct the measurements (Fig. 3) was connected to the U-tube in the sandbox. Electric heating element supply heat to the circulating water. A pump circulates the water, and a flow meter was used to measure the volume flow rate of circulating water. The uncertainty of the measured flow rate is ±2 % of full range. The inlet and outlet fluid temperature measurements were conducted using PT100 temperature resistance sensors in 1/10B accuracy class. The BHE axis temperature in several location were measured using RTDs in 1/3B accuracy class. Locations of these RTDs are illustrated in figure 4. All of these RTDs were connected to a high-accuracy data acquisition system (NI PCI-6289 and NI E3S Web of Conferences 46, 00017 (2018) https://doi.org/10.1051/e3sconf/20184600017 3 rd International Conference on Energy and Environmental Protection SCXI-1503). The temperature of sand was measured using 110 of DS18B20 temperature sensors with uncertainty of ±0.5 K. Locations of these temperature sensors are illustrated in Figure 4.
To simplify the computational domain, the polypropylene housing of the sandbox tank was neglected (Fig. 2). The author conducted estimation based on two experiments (P1 and P2) with different operating parameters. The time of both experiment equal 71 hours.

Mathematical model
During estimation procedure author decided to use two numerical models of heat transport in borehole heat exchanger. The first one was a quasi 3D heat transport in BHE model [10,11] developed in ANSYS software and the second one was a MDF model presented in work [11].

The heat transport in the ground
Heat transport in the closest vicinity of a borehole heat exchanger is described by a function of spatial coordinates and time. This issue is also connected with the flow of fluid. It was assumed that the advection heat transport in the ground can be neglected (due to a lack of a water in sandbox). Also, evaporation and seepage of water has been neglected. In this case the heat transport in the rock mass was derived from Fourier's law and the conservation of energy, which is presented as partial differential equation of the transient heat conduction (Eq. 1).

Heat transport in the BHE using the q3D model
In this model the grout area is treated as a fully three-dimensional volume modelled by utilizing SOLID 70 type elements, whereas the fluid is modelled by one-dimensional FLUID 116 type elements in ANSYS software. Equation (Eq. 1) presents a partial differential equation of conduction in soil as well as in grout, and equations (Eq. 2) and (Eq. 3) present the heat transport process in a U-tube.

( )
Despite such dimensional simplification, real heat transfer within all the elements of the exchanger has been considered. Representing a fluid in the form of a one-dimensional finite element allows to maintain temperature variation along the vertical axis of this element while neglecting temperature variation in the radial direction. Heat fluxes that are normal to the contact surface of a U-tube pipe and the borehole along the vertical axis are fully considered. The quasi-3-D model was presented in detail in a paper [10].

Heat transport in the BHE using the MDF model
In this one-dimension finite element model, equations (Eq. 4 -8) presented heat transport in the inlet pipe, outlet pipe and grout g1, g2, g3, respectively.
The proposed formulation allows to present the elements of a borehole as a 1D finite element that describes the pipes and grout. The presented model allows to maintain temperature variation along this element's vertical axis while neglecting the temperature variation in the radial direction. Heat fluxes normal to the contact surface of the borehole and rock mass along the vertical axis are fully considered.

Initial and boundary conditions
In order to solve the presented partial differential equations require assuming the initial and boundary conditions which were measured during the experiment test. It was assumed that the initial condition was connected with the influence of an undisturbed environment in the sandbox. As the initial condition, in the whole computational domain, the average ground temperature Tp1,s,g,f(t=0)=294.5 K for experiment P1 and Tp2,s,g,f(t=0)=296.95 K for experiment P2, was assumed. To simulate external sandbox conditions, for the all face  (Fig. 2), the convection boundary condition was assumed as in equation (Eq. 9), where Ta,P1,P2=296.75 K is the three day average air temperature in the laboratory hall.
( ) he assumed temperature at the inlet to the borehole heat exchanger U-tube was equal to the temperature at the inlet of the U-tube during the experiment test Ti,inlet(t)=Tinlet,P1,P2(t).
During measurements the volume flow rate was constant

Grid independent test
Author did a grid independence test to ensure that the results would converge (Fig. 5). Solving problem is transient, so grid independence was checked for several time steps. The numerical simulation was done with 255 000 finite elements for q3D model and 130 000 elements for MDF model.
The Crank-Nicolson algorithm was used to solve the partial differential equations and an adaptive transient time step interval was assumed.

Results of the parameters estimation
The number of numerical experiments were reduced by using DOE. An extended central composite face centered (CCF) design in which factors adopt values on five levels were used. For four factors x1, x2, x3 and x4 there is a need to conduct 49 numerical experiments. The calculations were performed using "ZEUS" supercomputer at ACK Cyfronet Krakow. The table 3 presents the plan of the conducted simulations with the input and output parameters. For the output parameters, the author assumed two parameters described by the equations (Eq. 10 -12). Y1P1 and Y1P2 described the average discrepancy in a time range 0; 8 t h ∈ and Y2P1 and Y2P2 in a time range  The response points (Tab. 3) were approximated using the Kriging algorithm, which leads to the designation of the multidimensional response surface (Fig. 6 -8  https://doi.org/10.1051/e3sconf/20184600017 3 rd International Conference on Energy and Environmental Protection model, respectively. The response surface depends on two other design parameter, x3 and x4 which are equal   = 1.5, 2.5,   = 1000, 1800. The discrepancy, varies between 2% and 30%.
In figure 8 the presented multidimensional response surface was obtained from the design parameters, x3, x4 and the output parameters Y1P1 for q3D and MDF models. The response surface depends on two other design parameter, x1 and x2 which are equal   = 0.15, 0.2 0.25,   = 750, 950. The discrepancy, varies between 2% and 20%.

Fig. 8. The response surface for q3D and MDF models
To estimate the parameters the concept of Pareto dominance in multiple-objective optimization on predicted response surface was used. The objective of optimization was to minimize the output parameters (Y1P1, Y2P1, Y1P2, Y2P2) for both models. To conduct optimization procedure the screening algorithm in ANSYS 15 was used. The 1000 samples was set to generate and three candidate points was generated by optimization algorithm. In figures 9 and 10 the trade-off charts which shows the Pareto-dominant solutions was presented.  The estimated parameters were verified by experiment data and simulation with different operating parameters than previous tests (P1 and P2). The discrepancies were calculated according equations (Eq. 11 -12) and equal   = 3.87,   = 1.8 for q3D model and   = 4.16,   = 1.62 for MDF model. In figure 11 the outlet and inlet fluid temperature differences for experiment and both model with estimated parameters were shown. The discrepancy according equation (Eq. 10) was also presented (Fig. 11).   Fig. 11. The outlet and inlet fluid temperature differences (left) and discrepancy according eq. (10) (right) The average discrepancy based on verification measurement for both heat transport models are below 5 %.

Conclusion
In this work, the estimation of the ground and BHE thermal properties has been investigated. Two numerical model of heat transfer in BHE were used to conduct estimation procedure. The work begins with a general determination of the balance equation for the flow and transport of heat within each element of the small laboratory sandbox. The BHE exchanger is modelled as a quasi-three dimensional and as a one-dimensional finite element with multiple degrees of freedom. The DOE is then used to obtain a response surface and conduct the multiple-objective optimization analysis. From the simulations mentioned above, the following conclusions were obtained: • Based on the DOE, RSM and multiple-objective optimization the material properties of BHE and soil was determined. For both models the discrepancies Y1P1 are below 4 %, Y2P1 are below 11.5 %, Y1P2 are below 6.5 % and Y2P2 are below 11.5 %.
• For the heat transport models, the thermal properties variations with temperature and moisture content were neglected therefore, the calculated discrepancies below 11.5 % are considered satisfactory.
• The same optimal parameters values was determined based on two different heat transport models.
• The average discrepancies below 5 % based on the verification measurements with different operating parameters was obtained.
• The DOE, RSM and multiple-objective optimization as estimation procedure is an alternative to TRT.
As a further direction of research, it is proposed to conduct estimation procedure for BHE in full scale and compare estimated parameters with laboratory measurements.
This research was supported in part by PLGrid Infrastructure.