Numerical simulation of ground heat exchangers based on dynamic thermal boundary conditions in solid zone

Accepted Manuscript Numerical simulation of ground heat exchangers based on dynamic thermal boundary conditions in solid zone Jun Wang, Enshen Long, Wen Qin PII:

S1359-4311(13)00326-8

DOI:

10.1016/j.applthermaleng.2013.04.055

Reference:

ATE 4789

To appear in:

Applied Thermal Engineering

Received Date: 14 December 2012 Revised Date:

7 March 2013

Accepted Date: 27 April 2013

Please cite this article as: J. Wang, E. Long, W. Qin, Numerical simulation of ground heat exchangers based on dynamic thermal boundary conditions in solid zone, Applied Thermal Engineering (2013), doi: 10.1016/j.applthermaleng.2013.04.055. This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. 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.

ACCEPTED MANUSCRIPT

Numerical simulation of ground heat exchangers based on dynamic thermal boundary conditions in solid zone Jun Wanga,*, Enshen Longa, Wen Qinb

b

College of Architecture and Environment, Sichuan University, Chengdu 610065, China

RI PT

a

Liaoning Urban and Rural Construction and Planning Design Institute, Shenyang 116000, China *Corresponding author: [email protected]

Abstract: Accurate calculation and prediction of solid zone (pipe, grout and soil) temperature

SC

variation caused by ground heat exchanger heat release or absorption are needed for evaluating

performance of ground source heat pump system. For this purpose, according to energy

M AN U

conservation, one numerical simulation model with dynamic thermal boundary conditions is established and verified with experimental data. As to single borehole and array boreholes, temperature variations in soil zone are assessed with the present model. Meanwhile, limitations of constant temperature and adiabatic boundary conditions are found and can be overcome with dynamic thermal boundary conditions. This study achieves improvement of numerical solution on heat transfer of ground heat exchanger and provides one favorable method for analyzing the

TE D

influence of ground heat exchanger operation on solid zone temperature. Keywords: ground source heat pump; ground heat exchanger; numerical simulation; dynamic

EP

thermal boundary condition

1. Introduction

AC C

Imbalance between heat release to soil in summer and heat absorption from soil in winter for ground source heat pump (GSHP) system can lead to soil temperature change, which is described as unbalance level. The GSHP may have multiple operation modes [1] and its long-term running can give rise to obvious deviation of soil temperature from original state, as a result of the imbalance. Moreover, this deviation can cause the GSHP efficiency to decline. Therefore, soil temperature recoverability is one important foundation to determine stable operation of the GSHP, which depends on heat accumulation in soil zone and reflects at annual change of the soil temperature. At present, there are two kinds of models for describing heat transfer of ground heat 1

ACCEPTED MANUSCRIPT

exchanger (GHE) and calculating soil temperature [2], including analytical solutions and numerical solutions. The analytical solutions, e.g. line source model [3-5] and cylindrical source model [6], are simple and convenient for application and have been improved time and again in

RI PT

previous research, although they behave a bit roughness to judge soil temperature response as a result of some assumptions. Bandyopadhyay et al. [7] established Laplace domain solutions for

equivalent single core of U-tube in grouted boreholes, from which average fluid temperature and borehole boundary temperature had been obtained by using Gaver-Stehfest numerical inversion

SC

algorithm. Yang et al. [8] proposed one updated two-region vertical U-tube GHE analytical model for system dynamic simulation of ground coupled heat pump. Li et al. [9] developed analytical solutions for continuous line and cylindrical-surface sources based on Jaeger’s instantaneous

M AN U

line-source solution and formed a new approach to model heat transfer of the GHE with composite media and complex geometry. Next, compared with analytical solutions, numerical solutions have advantages to fully describe mechanism and complexity of the GHE heat transfer, which involve finite difference method [10-13], finite element method [14] and finite volume method [15]. In recent years, some new or improved numerical simulation models have been established, e.g.

TE D

Florides et al. [16] developed a numerical model for simulation of energy flows and temperature changes in and around a borehole heat exchanger when fluid circulated through a U-tube. Of course, the numerical solutions may lead to computational cost rise. In order to reduce number of nodes and computation time, Bauer et al. [17] presented one three-dimensional (3D) numerical

EP

simulation model by using simplified thermal resistance and capacity models, which considered transient effects of heat and mass transports inside borehole.

AC C

In practice, both models and their thermal boundary conditions have impacts on calculation and prediction of soil temperature. Constant temperature or adiabatic is usually used as thermal boundary condition in solid zone (pipe, grout and soil) for numerical simulation [18-20]. However, for a long-term running or a large-scale application of the GHE, temperature or heat flux at the boundary of solid zone will change, especially at far-end surface of the solid zone. In this case, according to energy conservation among pipe heat transfer rate (Qin), heat consumption or release in solid zone (Qsink), boundary heat transfer rates at top surface (Qtop), bottom surface (Qbottom) and far-end surface (Qfar-end) of solid zone, shown in Fig. 1, thermal boundary condition of constant temperature or adiabatic may bring about inaccurate calculation of soil temperature and even 2

ACCEPTED MANUSCRIPT

make it gravely deviate from its actual value. Therefore, in order to accurately calculate or predict soil temperature change for a long-term operation or a large-scale application of the GHE, it’s necessary to improve numerical solution with the aid of new thermal boundary condition in the

RI PT

solid zone, which can express actual heat transfer of the GHE. The objective of this present work is to develop dynamic thermal boundary conditions in solid zone, establish one numerical simulation model based on energy conservation and achieve improvement of numerical solution, from which accurate calculation and prediction of solid zone

SC

temperature variation can be acquired, especially for a long-term running or a large-scale

application of the GHE. The established numerical solution is verified with experimental data and then used to evaluate soil temperature variation for single borehole and array boreholes,

M AN U

respectively.

2. Numerical simulation model of GHE

Mass, momentum, turbulence and energy conservation of water flow and the energy conservation of solid zone (pipe, grout and soil) are taken into consideration for numerical simulation of solid–fluid conjugated heat transfer of the GHE, which is shown in Fig. 1. As the

TE D

Reynolds number (Re) of water flow in GHE reaches 8,000-16,000, it can be recognized as fully turbulent pipe flow. The turbulence of water flow and heat transfer through the pipe wall are described by the standard k-ε model and the standard wall functions [21].The corresponding mathematic model of GHE heat transfer is as the following equation:

EP

∂ ( ρφ ) + div( ρU φ ) = div(Γφ gradφ ) + Sφ (1) ∂t where ф is the common variable; ρ is the water or solid region density; U is the water

AC C

velocity; Γф is the diffusion coefficient; Sф is the source term. Table 1 gives the expressions of

variables, diffusion coefficients and source terms for the water, pipe, grout and soil. In Table 1, ui is water velocity component at the direction i; Tf, Tp, Tg, Ts are temperature of the water, pipe,

grout and soil, respectively; µ and µt are molecular and turbulent viscosities, respectively; α is the

turbulent Prandtl number; λ denotes thermal conductivity; p is pressure of the water; g represents gravitational acceleration; Gk is the generation of turbulence kinetic energy due to the mean velocity gradients; Gb is the generation of turbulence kinetic energy due to buoyancy. In consideration of one kind of common initial condition that there exists non-uniform

3

ACCEPTED MANUSCRIPT

temperature distribution in the water, pipe, grout and soil, the initial condition is given as the following equations:

RI PT

Tf (i, j , k , 0) = F f (i, j , k )  Tp (i, j , k , 0) = F p (i, j , k )   (2) Tg (i, j , k , 0) = F g (i, j , k )  Ts (i, j , k , 0) = F s (i, j , k ) 

Next, according to energy conservation of the water, pipe, grout and soil, thermal boundary

SC

conditions should conform to requirements:

M AN U

Qin + Qsink + Qtop + Qbottom + Qfar-end = 0   Qtop = qtop Atop   (3) Qbottom = qbottom Abottom   Qfar-end = qfar-end Afar-end 

where qtop, qbottom and qfar-end are heat flux at top surface, bottom surface of the solid zone and far-end surface of the soil zone, respectively; Atop, Abottom and Afar-end are areas for corresponding positions.

TE D

Meanwhile, as impacts of the GHE running on temperature at bottom surface of the solid zone and far-end surface of the soil zone happen after a critical time ( tc ), running time ( t ) of the GHE needs to be divided into two stage. When t ≤ tc, bottom surface of the solid zone and far-end surface of the soil zone can be assumed adiabatic, qbottom = 0 and qfar-end = 0; when t > tc, qbottom and

EP

qfar-end are calculated as the following method that forms dynamic thermal boundary conditions in the solid zone, including the soil zone.

AC C

For t > tc, a scale factor ξ (0 ≤ ξ ≤ 1), acquired from thermal response test analysis, is employed to get the calculating methods of qbottom and qfar-end, expressions of which are as the following equations:

qbottom =

(Q + Qsink + Qtop ) Qbottom = −ξ in Abottom Abottom

(4)

qfar-end =

(Q + Qsink + Qtop ) Qfar-end = −(1 − ξ ) in Afar-end Afar-end

(5)

For the moment of t = tc, Qin, Qsink and qtop can be calculated according to numerical simulation results, then sum of Qbottom and Qfar-end can be obtained from equation (3). 4

ACCEPTED MANUSCRIPT

Correspondingly, values of qbottom and qfar-end can be derived from equation (4) and equation (5) for the moment of t = tc, respectively. While values of qbottom and qfar-end at the moment of t = tc are used as thermal boundary conditions for next time step numerical simulation, Qin, Qsink and qtop at t

RI PT

= tc + t can be acquired. Furthermore, the values of qbottom and qfar-end can be updated with above method at the moment of t = tc + t, then which are used as thermal boundary conditions for the next time step numerical simulation again. By means of iteration, qbottom and qfar-end can be updated

time and again and dynamic thermal boundary conditions are formed. Therefore, in order to obtain

Qin = cw mw [Tfin ( t −∆t ) − Tfout ( t −∆t ) ] Qsink = −[ ∫ ρ p ( Ωp

∂Tp

∂t

)(t −∆t ) dv +

∫

Ωg

qtop = −h(Ttop( t −∆t ) − Ta ( t −∆t ) )

ρg (

M AN U

firstly for the moment of t – t, which are as the following:

SC

the dynamic thermal boundary conditions in solid zone, Qin, Qsink and qtop should be determined

∂Tg

∂t

)( t −∆t ) dv +

∂Ts

∫ ρ ( ∂t ) s

( t −∆t )

(6)

dv] (7)

Ωs

(8)

In above equations, cw is the water specific heat capacity; mw denotes the water mass flow

TE D

rate; Tfin and Tfout are water temperature at inlet and outlet of the GHE, respectively; Ωp, Ωg and Ωs are corresponding volume at pipe, grout and soil; h is coefficient of convective heat transfer at the top surface of solid zone, Ttop is its temperature; Ta is outside air temperature.

EP

It’s worth pointing out that the role of tc is providing time node for changing constant thermal boundary conditions to dynamic thermal boundary conditions. As equation (3) accords with

AC C

energy conservation, the constant thermal boundary conditions can be seen as the special case of the established dynamic thermal boundary conditions. Therefore, the value of tc has no effect on the simulation results when it’s less than or equal to its physical truth value and more than or equal to one time step ( t) chosen for numerical simulation. If the physical truth value of tc can’t

be determined accurately in advance by experiment, tc needs to be assigned a small value for

meeting above requirement. When t or 2 t is less than the physical truth value of tc, tc = t or tc = 2 t is one reasonable choose. As the thermal diffusion in solid zone is slow, if t is no more than one minute, t or 2 t will be much less than the physical truth value of tc for general case, 5

ACCEPTED MANUSCRIPT

which is equal to the time that pipe heat transfers to bottom surface of the solid zone and far-end surface of the soil zone. In addition, the assumption in the present model, that boundary heat flux are uniform

RI PT

distribution at top surface, bottom surface of the solid zone and far-end surface of the soil zone, doesn’t affect correctness of numerical simulation results to reflect actual heat transfer of the GHE

for a long period operation. Moreover, it makes the present model be convenient and practical for engineering application.

SC

3. Verification of the numerical simulation model

Experimental data in reference [22] are used for validation of the established numerical simulation model in this study. A 3D region is chosen and has dimensions of 5 m (width) × 5 m

M AN U

(length) × 80 m (depth), including zones of water, high density polyethylene (HDPE) pipe, grout and soil. Single U-shaped type is adopted in this GHE, with inner diameter of 25 mm and outer diameter of 32 mm. Table 2 gives physical properties of the water and solid materials. Initial temperature of all zones is 17.8 ℃. Water flow rate is 0.72 m3h-1 and its temperature (IWT) at the GHE inlet refers to the measured data. Turbulent intensity of the GHE inlet is calculated through

TE D

0.16Re-1/8. Boundary condition at the GHE outlet is single direction outflow. Heat flux at far-end surface of the soil zone is obtained according to the method provided in this study. Influence of heat flux at bottom surface of the solid zone can be ignored for running time of 96 hours and ξ =0. Meanwhile, heat flux at top surface of the solid zone is determined with equation (8), h equals to

EP

13.9 Wm-2℃-1 and outside air temperature (OAT) is derived from the measured value. The finite volume method is utilized to discretize equations of this present model in the

AC C

simulating region. Evaluation of the advection terms of Navier-Stokes equations, turbulent transport equations and energy equations is achieved by the QUICK scheme. The central differencing algorithm is used for two-order implicit scheme. The SIMPLEC algorithm is applied to the pressure-velocity coupling. Grid-independent solution is tested by using three levels of non-structured grid sizes, the final meshes used in this simulation are demonstrated in Fig. 2. Unsteady simulation proceeds to obtain variations of the GHE outlet water temperature (OWT) and heat exchange quantity per unit of GHE depth with running time. Fig. 3 gives the comparison of recorded data to simulated values for GHE outlet water temperature and heat exchange quantity per unit of GHE depth, the percentages of error in 6

ACCEPTED MANUSCRIPT

simulation results are less than 5%, which shows their well consistency and indicates that the present model and its dynamic thermal boundary conditions can simulate heat transfer of the GHE accurately.

RI PT

4. Numerical evaluation of ground temperature variation 4.1. Single borehole

GHE of single borehole is chosen for analyzing influence generated by three kinds of thermal

boundary conditions on soil temperature distribution and variation, including constant temperature

SC

(far-end surface temperature of the soil zone Tfar-end = 17.8 ℃), adiabatic (qfar-end = 0), dynamic

thermal boundary conditions. As difference among mean water temperature at different GHE depth and temperature gradient along ground depth are small for many GSHP projects, vertical

M AN U

heat transfer of the GHE is neglected and quarter of one 2D circular region with radius r = 3m is selected for simulating analysis, shown in Fig. 4. Meanwhile, diameter of the borehole is 0.12 m. Symmetrical boundary conditions are used at two interior faces of the region and chosen thermal boundary condition is applied to its exterior face. Moreover, specified heat flux is fed to surface of the borehole, corresponding to heat exchange quantity per unit depth of 20 W/m (low design load)

TE D

or 60 W/m (high design load) for cooling season and - 20 W/m (low design load) or - 60 W/m (high design load) for heating season. Everyday running to standby time ratios (RSR) of 24/0 and 8/16 are considered as two operating modes of the GSHP system in cooling and heating seasons. In addition, physical properties, initial temperature of soil zone and numerical method are the

EP

same as Section 3. The operation of GHE contains four months cooling season and two months recovery season, four months heating season and two months recovery season.

AC C

Fig. 5 shows soil temperature variations at point A1 (r = 0.6 m) and point B1 (r = 1.5 m) with running time of the GHE for low design load of 20 W/m and RSR = 24/0 in cooling season. It can be observed that difference of point A1 temperature acquired by these three kinds of thermal

boundary conditions is negligible, as which is far from the far-end boundary. However, obvious difference of point B1 temperature occurs when different thermal boundary conditions are applied.

After 71 days operation, the value of point B1 temperature under adiabatic boundary condition begins to be higher than those for other two kinds of thermal boundary conditions. After 11 days in recovery season, the point B1 temperature calculated with Tfar-end = 17.8 ℃ starts to be less than that corresponding to dynamic qfar-end. At the end of recovery season, the values of point B1 7

ACCEPTED MANUSCRIPT

temperature for thermal boundary conditions of dynamic qfar-end and Tfar-end = 17.8 ℃ are 0.86% and 2.78% lower than that obtained by adiabatic boundary condition, respectively. When design load rises to 60 W/m and RSR = 24/0 in cooling season, evolution of soil

RI PT

temperature at point A1 and point B1 with running time of the GHE are given in Fig. 6. The difference of point A1 temperature derived from application of different thermal boundary

conditions still can be ignored. But deviation of point B1 temperature got by adiabatic boundary condition from those calculated with other two kinds of thermal boundary conditions emerges on

SC

the sixty-third day ahead. Moreover, on the one hundred and twenty-third day, point B1

temperature under dynamic qfar-end begins to be higher than that for Tfar-end = 17.8 ℃. In addition, at the end of recovery season, values of point B1 temperature got from thermal boundary conditions

M AN U

of dynamic qfar-end and Tfar-end = 17.8 ℃ are 2.38% and 7.59% lower than that for adiabatic boundary condition, respectively. Therefore, increase of design load may lead difference among values of soil temperature calculated with these three kinds of thermal boundary conditions to be more obvious.

For low design load of - 20 W/m and RSR = 24/0 in heating season, variations of soil

TE D

temperature at point A1 and point B1 with running time of the GHE are provided in Fig. 7. It can be found that obvious temperature difference is only appearing point B1 for these three kinds of thermal boundary conditions. Meanwhile, the temperature at point B1 under adiabatic boundary condition starts to be lower than those derived from other two kinds of thermal boundary

EP

conditions after 63 days operation. In addition, after 5 days in recovery season, the point B1 temperature corresponding to Tfar-end = 17.8

starts to be higer than that calculated with dynamic

AC C

qfar-end. At the end of recovery season, the point B1 temperature obtained from thermal boundary conditions of dynamic qfar-end and Tfar-end = 17.8

are 0.98% and 3.06% higher than that for

adiabatic boundary condition, respectively. Fig. 8 gives variations of soil temperature at point A1 and point B1 with running time of the

GHE for high design load of - 60 W/m and RSR = 24/0 in heating season. Main characteristics are similar with those shown in Fig. 7. But time is advance to the 58th day for the temperature at point B1 under adiabatic boundary condition being less than those obtained from other two kinds of thermal boundary conditions. Moreover, the point B1 temperature at dynamic qfar-end begins to be 8

ACCEPTED MANUSCRIPT

lower than that calculated with Tfar-end = 17.8

on the one hundred and twenty-fourth day. In

addition, the point B1 temperature derived from thermal boundary conditions of dynamic qfar-end and Tfar-end = 17.8

are 3.23% and 9.29% more than that for adiabatic boundary condition at the

RI PT

end of recovery season, respectively. So the thermal boundary conditions have more obvious influence on soil temperature in heating season than those in cooling season.

For RSR = 8/16 in cooling season, Fig. 9 and Fig. 10 give evolutions of soil temperature at point A1 and point B1 with running time of the GHE for low design load of 20 W/m and high

SC

design load of 60 W/m, respectively. In addition to some similar characteristics with those depicted in Fig. 5 and Fig. 6, a few new features can be found for intermittent operation of the

M AN U

GHE, e.g., time for difference among values of soil temperature obtained by different thermal boundary conditions appearing is delayed at point B1. Moreover, when recovery season ends as to the low design load of 20 W/m, the values of point B1 temperature for dynamic qfar-end and Tfar-end = 17.8 ℃ are 0.20% and 0.90% lower than that got from adiabatic boundary condition, respectively. On the other hand, for high design load of 60 W/m, soil temperature at point A1 under Tfar-end = 17.8

TE D

starts to be less than those at the other two kinds of thermal boundary

conditions on the second day in cooling season. Meanwhile, the soil temperature of point A1 at adiabatic boundary condition is higher than that for dynamic qfar-end after the forty-ninth day of cooling season. In addition, at the end of recovery season, values of point B1 temperature for

EP

dynamic qfar-end and Tfar-end = 17.8

are 1.15% and 3.60% lower than that obtained from adiabatic

boundary condition, respectively. Therefore, impacts of thermal boundary conditions on soil

AC C

temperature become weak for intermittent operation of the GHE. For RSR = 8/16 in heating season, Fig. 11 and Fig. 12 show soil temperature variation at

point A1 and point B1 with running time of the GHE for low design load of - 20 W/m and high design load of - 60 W/m, respectively. Compared with the results at continuous operation of the GHE, thermal boundary conditions also produce weaker influence on soil temperature under intermittent operation of the GHE. For the low design load of - 20 W/m, the B1 temperature calculated with dynamic qfar-end and Tfar-end = 17.8

are 0.25% and 1.02% higher than that under

adiabatic boundary condition at the end of recovery season, respectively. Meanwhile, the 9

ACCEPTED MANUSCRIPT

corresponding percentages reach 1.3% and 3.96% for the high design load of - 60 W/m. Fig. 13 presents temperature fields of the soil zone for high design load of 60 W/m in cooling season. It’s obvious that high temperature occurs in the vicinity of borehole. Moreover, when RSR

0.38% higher than those for dynamic qfar-end and Tfar-end = 17.8

RI PT

is 24/0, mean temperature of the soil zone under adiabatic boundary condition are 0.14% and , respectively. As for RSR = 8/16,

the above corresponding percentages become 0.07% and 0.20%. On the other hand, for high

design load of - 60 W/m in heating season, temperature distributions of the soil zone are given in

SC

Fig. 14. Under this circumstance, low temperature is observed in the vicinity of borehole. Meanwhile, for RSR = 24/0, mean temperature of the soil zone derived from dynamic qfar-end and are 0.15% and 0.46% more than that under adiabatic boundary condition. When

M AN U

Tfar-end = 17.8

RSR is 8/16, the above corresponding percentages are 0.08% and 0.22%. The above results indicate that soil temperature change calculated with adiabatic boundary condition is affected by heat release or absorption of the GHE at the highest level, compared with the other two kinds of thermal boundary conditions. 4.2. Array boreholes

TE D

In order to model and revalue soil temperature variation for array boreholes, the GHE applied in one practical project is chosen as analysis object. Fig. 15 shows annual cooling and heating loads distribution of the project, positive and negative values represent cooling load and heating

EP

load, respectively. Annual accumulative cooling load of the project is 882.6 MWh and its annual accumulative heating load equals 804.9 MWh, so the thermal imbalance ratio is 8.8%. Each

AC C

borehole has depth of 70 m and 0.11 m diameter. For total area (80 m × 80 m) of the array boreholes with horizontal spacing of 5 m, vertical heat transfer of the GHE is ignored and one 2D corner of 40 m × 40 m is used to figure out the soil zone temperature distribution and variation, presented in Fig. 16. Its two interior faces apply symmetrical boundary conditions and chosen thermal boundary conditions are used to the exterior faces. Meanwhile, specified heat flux is fed to surface of each borehole, which is determined according to heat exchange quantity per unit depth corresponding to cooling and heating loads. Density, specific heat and thermal conductivity of the soil are 1352 kgm-3, 2632 Jkg-1K-1 and 1.8 Wm-1K-1, respectively, initial temperature of which is 14.1

. In addition, the applied numerical method is the same as Section 3. 10

ACCEPTED MANUSCRIPT

When first operation of the GHE is used for cooling, soil temperature variations in three years at point A2 (x = 1 m, y = 1 m) and point B2 (x = 19 m, y = 19 m) are shown in Fig. 17. It can be found that, except the last 102 days, values of the soil temperature at point B2 behave mostly

RI PT

identical for these three kinds of thermal boundary conditions, as which is far from far-end boundary and impacted by alteration of thermal boundary conditions slightly. Meanwhile, the soil temperature rise at point B2 reaches 0.43

for dynamic qfar-end and Tfar-end = 14.1

and 0.49

for adiabatic (qfar-end = 0) at end of the third year. On the other hand, the soil temperature at point

SC

A2 obtained by qfar-end = 0 increases continually with running time, rise of which attains 0.89

after three years operation. When dynamic qfar-end is applied, fluctuating change of the soil

M AN U

temperature at point A2 caused by cooling and heating loads emerges. Moreover, the soil temperature at point A2 calculated with Tfar-end = 14.1

varies flatly during the GHE running. In

addition, there is not obvious rise of the soil temperature at point A2 for dynamic qfar-end and Tfar-end = 14.1

.

As first operation of the GHE is applied for heating, characteristics of soil temperature

TE D

change at point A2 and point B2 are similar to the above case, shown in Fig. 18. But it should be noticed that deviation of the soil temperature at point A2 got from qfar-end = 0 from values of that acquired by dynamic qfar-end and Tfar-end = 14.1

decreases to a low level. Meanwhile, at end of , a bit higher

EP

the third year, rise of the soil temperature at point B2 under qfar-end = 0 reaches 0.65 than that in above case.

AC C

Fig. 19 gives temperature fields of the soil zone for these two kinds of operating mode, corresponding to results at end of the GHE three years running. For the case of cooling first, it can be found that the mean temperature of soil zone for adiabatic boundary condition is 0.09% and 0.11% higher than those at thermal boundary conditions of dynamic qfar-end and Tfar-end = 14.1

,

respectively. But for heating first, the former is 0.07% and 0.06% less than the latter two kinds of results, correspondingly. Again, it’s demonstrating that adiabatic boundary condition may make prediction be on the high side for influence of the GHE heat release or absorption on soil zone temperature. Meanwhile, constant temperature boundary condition may underestimate the

11

ACCEPTED MANUSCRIPT

influence. Therefore, the dynamic thermal boundary conditions established in this study are needed to obtain reasonable results of calculation and prediction. As a matter of fact, it’s the root cause that accurate temperature gradient variation at the boundary of soil zone can’t be acquired

RI PT

with adiabatic and constant temperature used as thermal boundary conditions, especially for a long-term operation or a large scale application of the GHE. 5. Conclusions

In order to obtain accurate calculation and prediction of solid zone (pipe, grout and soil)

SC

temperature variation caused by the GHE heat release or absorption, one numerical simulation

model with dynamic thermal boundary conditions is established in this study. The model is verified with experimental data and then used to evaluate soil temperature variations for single

M AN U

borehole and array boreholes. As to simulating results derived from thermal boundary conditions of constant temperature, adiabatic and dynamic thermal boundary conditions, the first one generates underestimation for influence of the GHE operation on solid zone temperature. Meanwhile, overestimation of this influence appears for the second one. In addition, a long-term operation or a large scale application of the GHE may make this phenomenon be more obvious.

TE D

Therefore, the present dynamic thermal boundary conditions are useful for overcoming limitations of constant temperature and adiabatic boundary conditions and acquiring reasonable results of solid zone temperature. 6. Acknowledgements

EP

The authors gratefully acknowledgement the financial support from China Postdoctoral Science Foundation under Grant No. 2012M511930 and the National Natural Science Foundation

AC C

of China under Grant No. 51178282. References

[1] Jalaluddin, A. Miyara, Thermal performance investigation of several types of vertical ground heat exchangers with different operation mode, Appl Therm Eng, 33-34 (2012) 167-174. [2] Y. Yuan, X. Cao, L. Sun, B Lei, N. Yu, Ground source heat pump system: A review of simulation in China, Renew Sust Energ Rev, 16 (2012) 6814-6822. [3] W. Kelvin, Mathematical and Physical Papers, Cambridge: Cambridge University Press, 1 (1882) 100-106. [4] L. Ingersoll, H. Plass, Theory of the ground heat pipe heat source for the heat pump, ASHRAE 12

ACCEPTED MANUSCRIPT

Trans, 20 (1948) 119-122. [5] L. Ingersoll, O. Zoeble, A. Ingersoll, Heating conduction with engineering, geological and other applications, New York: McGRaw-Hill, 1954.

RI PT

[6] S. Kavanaugh, Simulation and experimental verification of vertical ground-couple heat pump systems, PhD thesis, OSU, 1984.

[7] G. Bandyopadhyay, W. Gosnold, M. Mann, Analytical and semi-analytical solutions for short-time transient response of ground heat exchangers, Energy Build, 40 (2008) 1816-1824.

SC

[8] W. Yang, M. Shi, G. Liu, Z. Chen, A two-region simulation model of vertical U-tube ground heat exchanger and its experimental verification, Appl Energ, 86 (2009) 2005-2012.

[9] M. Li, A. Lai, New temperature response functions (G functions) for pile and borehole ground

M AN U

heat exchangers based on composite-medium line-source theory, Energy, 38 (2012) 255-263. [10] V. Mei, S. Fisher, Vertical concentric tube ground coupled heat exchanger, ASHRAE Trans, 89 (1983) 391-406.

[11] V. Mei, C. Emerson, New approach for analysis of ground-coil design for applied heat pump system, ASHRAE Trans, 91 (1985) 1216-1224.

TE D

[12] V. Mei, Theoretical heat pump ground coil analysis with variable ground far-field boundary conditions, AICHE J, 92 (1986) 22-25.

[13] T. Lei, Development of a computational model for a ground-coupled heat exchanger, ASHRAE Trans, 99 (1993) 149-159.

EP

[14] N. Muraya, Numerical modeling of the transient thermal interference of vertical U-tube heat exchanger, PhD thesis, TAMU, 1995.

AC C

[15] C. Yavuzturk, J. Spitler, Comparative study of operating and control strategies for hybrid ground source heat pump systems using a short time step simulation model, ASHRAE Trans, 106 (2000) 192-209.

[16] G. Florides, P. Christodoulides, P. Pouloupatis, An analysis of heat flow through a borehole heat exchanger validated model, Appl Energ, 92 (2012) 523-533. [17] D. Bauer, W. Heidemann, H. Diersch, Transient 3D analysis of borehole heat exchanger modeling, Geothermics, 40 (2011) 250-260. [18] Z. Gu, Y. Wu, Z. Tang, C. Ma, Numerical simulation of heat transfer of a U-tube ground heat exchanger, J Eng Thermophys-rus, 27 (2006) 313-315. 13

ACCEPTED MANUSCRIPT

[19] C. Guan, Z. Wan, P. Hu, Numerical analysis on multilayer geotechnical temperature field of GSHP buried tube, J Wuhan Inst Tech, 33 (2011) 42-44. [20] J. Ma, Z. Zheng, Simualtion of heat transfer capacity of single and double U-type heat

RI PT

exchangers, J HVAC, 42 (2012) 108-112. [21] B. Launder, D. Spalding, The numerical computation of turbulent flows, Comput Meth Appl Mech Eng, 3 (1974) 269–89.

[22] R. Zhang, X. Zhang, L. Dong, Impact of flow rate on ground heat exchanger performance

SC

according to experiment, 7th Proceedings of International Symposium on Heating, Ventilation and

AC C

EP

TE D

M AN U

Air Conditioning, (2011) 901-907.

14

ACCEPTED MANUSCRIPT

Highlights Numerical simulation model with dynamic thermal boundary conditions is established.




Soil zone temperature variations are assessed with this present model.




Limitations of constant temperature and adiabatic boundary conditions are overcome.




This study achieves improvement of numerical solution on ground heat exchanger.

AC C

EP

TE D

M AN U

SC

RI PT




ACCEPTED MANUSCRIPT

Figure Caption Fig. 1. Solid–fluid conjugated heat transfer of the GHE. Fig. 2. Meshes used for the numerical simulation ( 1,127,175 nodes for the triple U-shaped): (a)

RI PT

view of the total mesh and (b) top view of the mesh. Fig. 3. Comparison between measured data and simulated values of the GHE outlet water temperature and heat exchange quantity per unit of GHE depth.

Fig. 4. One 2D circular region selected for simulation: (a) geometry size and (b) grid division.

SC

Fig. 5. Soil temperature variations at point A1 and point B1 at low design load of 20 W/m and RSR = 24/0.

Fig. 6. Soil temperature variations at point A1 and point B1 at high design load of 60 W/m and

M AN U

RSR = 24/0.

Fig. 7. Soil temperature variations at point A1 and point B1 at low design load of -20 W/m and RSR = 24/0.

Fig. 8. Soil temperature variations at point A1 and point B1 at high design load of -60 W/m and RSR = 24/0.

= 8/16.

TE D

Fig. 9. Soil temperature variations at point A1 and point B1 at low design load of 20 W/m and RSR

Fig. 10. Soil temperature variations at point A1 and point B1 at high design load of 60 W/m and RSR = 8/16.

RSR = 8/16.

EP

Fig. 11. Soil temperature variations at point A1 and point B1 at low design load of -20 W/m and

AC C

Fig. 12. Soil temperature variations at point A1 and point B1 at high design load of -60 W/m and RSR = 8/16.

Fig. 13. Temperature fields of the soil zone at high design load of 60 W/m: (1) Tfar-end = 17.8℃ and

RSR = 24/0; (2) Adiabatic and RSR = 24/0; (3) Dynamic qfar-end and RSR = 24/0; (4) Tfar-end = 17.8℃ and RSR = 8/16; (5) Adiabatic and RSR = 8/16; (6) Dynamic qfar-end and RSR = 8/16. Fig. 14. Temperature fields of the soil zone at high design load of -60 W/m: (1) Tfar-end = 17.8℃ and RSR = 24/0; (2) Adiabatic and RSR = 24/0; (3) Dynamic qfar-end and RSR = 24/0; (4) Tfar-end =

ACCEPTED MANUSCRIPT

17.8℃ and RSR = 8/16; (5) Adiabatic and RSR = 8/16; (6) Dynamic qfar-end and RSR = 8/16. Fig. 15. Annual cooling and heating loads distribution. Fig. 16. The 2D region of array boreholes used for simulation: (a) geometry size and (b) grid

RI PT

division. Fig. 17. Soil temperature variations at point A2 and point B2 for cooling first. Fig. 18. Soil temperature variations at point A2 and point B2 for heating first.

Fig. 19. Temperature fields of the soil zone: (1) Tfar-end = 14.1℃ and cooling first; (2) Adiabatic

SC

and cooling first; (3) Dynamic qfar-end and cooling first; (4) Tfar-end = 14.1℃ and heating first; (5)

AC C

EP

TE D

M AN U

Adiabatic and heating first; (6) Dynamic qfar-end and heating first.

ACCEPTED MANUSCRIPT

Table 1 Expressions of ф, Γф and Sф.

φ

Γφ

Sφ

Mass conservation

1

0

0

Momentum conservation

ui

µ+µt

−

Energy conservation

Tf

(µ+µt)/αT

0

Turbulent energy

k

(µ+µt)/αk

Gk + Gb − ρε

ε

(µ+µt)/αε

ε

turbulent energy

k

Energy conservation

Tp

λp

0

Grout

Energy conservation

Tg

λg

0

Soil

Energy conservation

Ts

λs

0

(C1ε Gk + C3ε Gb − C2ε ρε )

AC C

EP

TE D

M AN U

Pipe

RI PT

Dissipation rate of

∂u j ∂p ∂ + [( µ + µ t ) ] + ρ gi ∂xi ∂x j ∂x j

SC

Water

Equations

ACCEPTED MANUSCRIPT

Table 2 Physical properties of water and solid materials [22]. Density

Specific heat -1

Thermal conductivity

( kgm )

( Jkg K )

( Wm-1K-1 )

Water

1000

4180

0.60

HDPE

950

2300

0.44

Grout

1860

1200

1.60

Soil

1925

1400

1.40

RI PT

-1

AC C

EP

TE D

M AN U

SC

-3

AC C

EP

TE D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

AC C

EP

TE D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

AC C

EP

TE D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

AC C

EP

TE D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

AC C

EP

TE D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

AC C

EP

TE D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

AC C

EP

TE D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

AC C

EP

TE D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

AC C

EP

TE D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

AC C

EP

TE D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

AC C

EP

TE D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

AC C

EP

TE D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

AC C

EP

TE D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

AC C

EP

TE D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

AC C

EP

TE D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

AC C

EP

TE D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

AC C

EP

TE D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

AC C

EP

TE D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

AC C

EP

TE D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT