Dynamic Heat Transfer Modeling and Validation of Super-Long Flexible Thermosyphons for Shallow Geothermal Applications

Abstract

In comparison to borehole heat exchangers that rely on forced convection, super-long thermosyphons offer a more efficient approach to extracting shallow geothermal energy. This work conducted field tests on a super-long flexible thermosyphon (SFTS) to evaluate its heat transfer characteristics. The tests investigated the effects of cooling water temperature and the inclination angle of the condenser on the start-up characteristics and steady-state heat transfer performance. Based on the field test results, the study proposed a dynamic heat transfer modeling method for SFTSs using the equivalent thermal conductivity (ETC) model. Furthermore, a full-scale 3D CFD model for geothermal extraction via SFTS was developed, taking into account weather conditions and groundwater advection. The modeling validation showed that the simulation results aligned well with the temperature and heat transfer power variations observed in the field tests when the empirical coefficient in the ETC model was specified as 2. This work offers a semi-empirical dynamic heat transfer modeling method for geothermal thermosyphons, which can be readily incorporated into the overall simulation of a geothermal system that integrates thermosyphons.

1. Introduction

As global warming and environmental challenges intensify, developing and utilizing renewable energy has become crucial in addressing climate change [1]. Among various renewable energy sources, geothermal energy has garnered significant attention for its clean, sustainable nature [2]. In particular, shallow geothermal energy [3] stands out due to its abundant resources and low extraction costs, making it a key focus of contemporary energy research and application.

Traditionally, shallow geothermal energy is extracted using borehole heat exchangers (BHEs) coupled with heat pumps, which circulate a heat carrier fluid through buried pipe loops driven by circulation pumps [4,5]. However, this approach incurs significant energy consumption, with the circulation pump accounting for 10–20% of total system energy use.

The thermosyphon [6,7] is a passive heat transfer unit with highly equivalent thermal conductivity, offer a promising alternative for shallow geothermal energy utilization. By leveraging phase change heat transfer and passive fluid circulation, thermosyphons efficiently transfer heat without additional energy consumption, drastically reducing operational energy costs and maintenance demand [8]. Their potential spans diverse applications, including pavement de-icing [9,10,11] and spacing heating [12,13,14].

For the pavement de-icing, the use of thermosyphons is free of energy consumption. Zhang et al. [9] demonstrated the feasibility of a thermosyphon anti-freezing system for the airport runways at Beijing Daxing International Airport. The runways could be maintained free of ice and snow when the air temperature was above −4 °C. Zorn et al. [10] developed an CO₂-thermosyphon system to keep an asphalt pavement from icing in Bad Waldsee, Germany. The experimental results showed that the surface was kept above 0 °C under the test conditions. In addition, Yu et al. [11] integrated a multi-loop CO₂-thermosyphon with a bridge deck in Texas, USA, and the thermosyphon was connected with a ground heat exchanger to provide heat for de-icing. It was experimentally proven that the bridge deck surface was free of ice and snow as long as the air temperature was above −6.2 °C. Further, ground source heat pumps (GSHPs) with the incorporation of thermosyphons in borehole heat exchangers have been experimentally and theoretically investigated for providing heating to buildings [12]. The recent investigations revealed that the thermosyphon–GSHP system has a higher overall efficiency when compared with the traditional U-type or coil-type GSHP system. The long-term numerical analysis by Lim et al. [13] showed that a total energy consumption reduction of up to 10.3% was reached with the thermosyphon–GSHP system than with the traditional one based on the weather data in Seoul, South Korea. Additionally, super-long thermosyphon coupled GSHP systems applied in deep geothermal energy extraction for space heating in Taiyuan, China, have demonstrated significant application potential [14].

Despite the low maintenance advantages of thermosyphons in shallow geothermal systems, their super-long length leads to high fabrication, transportation and installation costs, limiting commercialization [15]. To address this, we previously proposed a super-long flexible thermosyphon (SFTS) made by metal corrugated pipes [15]. Field tests of a 32 m-long SFTS revealed that, in addition to its flexibility, the corrugated pipe prevents the formation of a high liquid column in the evaporator, which might cause a stagnant zone at the bottom. This design improves heat transfer in the evaporator of SFTS, suggesting strong application potential.

To design a shallow geothermal system integrated with thermosyphons, an accurate prediction method for the heat transfer performance in the thermosyphons is necessary and essential. Many numerical efforts have been made in this regard. Hartmann et al. [16] developed a transient 2D numerical model to investigate the effects of grout and pipe materials on the heat extraction of an 80 m-long geothermal thermosyphon. Ebeling et al. [17] numerically studied the thermal performance of a 368 m-long geothermal thermosyphon using CO₂ as the working fluid. In addition, Wang et al. [18] developed a transient CFD model based on the VOF (Volume-of-Fluid) multiphase method to reveal the phase change and flow behaviors in a geothermal thermosyphon. However, these models involve complex phase change heat transfer and fluid flow, making the modeling steps difficult and computationally expensive, hindering their application for overall geothermal system performance analysis.

Alternatively, semi-empirical methods, such as the thermal-resistance network model, are commonly used for steady-state heat transfer analysis of heat pipes or thermosyphons. Ozsoy et al. [19] developed a thermal-resistance model for geothermal thermosyphons to investigate the thermal performance under various types of working fluids, evaporator lengths and ground temperatures. In our previous works [20], the system performance of an ice and snow melting system coupled with thermosyphon was numerically investigated for varying arrangement intervals, evaporator lengths and air conditions, in which the thermosyphons were also modeled by the thermal-resistance model.

Another similar semi-empirical method treated the heat pipes or thermosyphons as thermal superconductors, referred to as the equivalent thermal conductivity (ETC) model, which was often used for the analysis of heat pipe integrated systems as well. Dillig et al. [21] numerically studied the thermal effects of using a planar high-temperature heat pipe as the interconnector in a solid oxide cell stacks, and the heat pipe interconnector in the CFD model was modeled as a solid with constant thermal conductivity, up to 15,000 W/(m·K). These semi-empirical methods for thermosyphons or heat pipes are easier to implement for overall system analysis, but they generally assume steady-state operation, where heat transfer performance remains constant under specific conditions.

However, in a geothermal utilization system, the operating state of thermosyphons varied with weather or hydrogeological conditions. Dynamic heat transfer modeling of thermosyphons based on a simple semi-empirical method is more convenient and useful for the geothermal system design and analysis. However, there is limited associated research work.

In this work, field tests were conducted on a 32 m-long SFTS to investigate the heat transfer characteristics. According to the heat transfer characteristics of SFTSs, we proposed a dynamic heat transfer modeling method based on the ETC model. And then, a full-scale 3D CFD model for geothermal extraction via an SFTS considering weather conditions and groundwater advection was developed. The simulation results were validated against field test data for both the start-up and steady-state heat transfer. This work offers a simple semi-empirical modeling method, easily applicable for the overall analysis of geothermal system integrating thermosyphons.

2. Experimental Analysis of the SFTS

2.1. Experimental Setup and Test Procedures

A 32 m-long flexible thermosyphon fabricated with 304 stainless steel, using ammonia as the working fluid, was field tested in this work. It was made up of a smooth pipe (2 m in length) and a corrugated pipe (30 m in length), as depicted in Figure 1a. The smooth pipe was set as the condenser, with an outer and inner diameters of 22 and 19 mm, respectively. The corrugated pipe was inserted into a borehole as the evaporator, while its near-surface part was wrapped by thermally insulated material with a length of 2 m. The outer and inner diameters of corrugated pipe were 25 and 18 mm, respectively. In addition, a 1 mm-thick woven metal wire was wrapped around the outer wall of the corrugated pipe to enhance its mechanical strength. The detailed structure of the corrugated pipe is also illustrated in Figure 1a. The filling ratio of ammonia was 37%, which was defined as the ratio of the liquid ammonia volume to the total volume of the corrugated pipe.

The schematic of the experimental setup is shown in Figure 1b. To evaluate the dynamic heat transfer performance of the SFTS, a shell heat exchanger was installed at the condenser to facilitate continuous water cooling. The temperature of cooling water was regulated by a thermostatic bath. A peristaltic pump was used to maintain the flow rate of circulated water. A photograph of the experimental setup is provided in Figure 2. Ten T-type thermocouples were employed to monitor the temperature variations along the SFTS during tests, with T₁–T₆ at the evaporator, T₇ at the adiabatic section, T₈–T₁₀ at the condenser, as seen in Figure 1c. Additionally, two T-type thermocouples were employed to monitor the inlet and outlet temperatures of the circulating water, labelled as T_in and T_out. The temperature data were recorded using a data logger at a frequency of 1 Hz. The details of the devices used in experiments are summarized in

Table 1. Details of employed devices.

Items	Manufacturer	Type	Parameters
Peristaltic pump	Shenchen, Shanghai, China	V6	0.07–6000 mL/min; flow rate accuracy: 0.5%
Thermostatic bath	Teer, Zhengzhou, China	DY20/20	−20 to 99 °C (±0.1 °C); volume: 20 L; maximum power: 3000 W
Data logger	Anbai, Changzhou, China	AT4532	32-channel; resolution: 0.1 °C; scanning rate: 32 channel/s
Thermocouples	Shenhua, Beijing, China	T-type	−100–300 °C; accuracy: ±0.1 °C

.

In this work, we investigated the effects of cooling temperature (5.5–10.5 °C) and the inclination angle of the condenser (0–90°) on the heat transfer characteristics of SFTSs. The inclination angle was defined as the angle between the condenser axis and the ground surface. During each test, the inclination angle of the condenser and the flow rate of the cooling water v_cooling were maintained invariable, while the cooling water temperature T_cooling in the thermostatic bath was changed from 5.5 to 10.5 °C in increments of 1 °C. v_cooling was controlled at 900 mL/min for all tests in this work. It should also be noted that the flow rate and temperature of the cooling water were calibrated with a precision of 5% and 0.1 °C, respectively.

2.2. Data Reduction and Error Analysis

The practical heat transfer power Q_p of the SFTS was calculated based on the temperature rise of cooling water in the cooling shell heat exchanger.

$$Q_{\mathrm{p}}=c_p\rho v_{\mathrm{cooling}}(T_{\mathrm{out}}-T_{\mathrm{in}}-\Delta T_{\mathrm{ambient}})$$

(1)

where c_p denotes the specific heat of cooling water, ρ is the density of water and v_cooling represents the flow rate of cooling water. T_in and T_out are the water temperature at the outlet and inlet of the cooling shell heat exchanger, respectively. ΔT_ambient results from the thermal influence of ambience, which had been evaluated in our previous work [15].

The steady-state thermal performance of the SFTS under various conditions was assessed by the total thermal resistance R_t, using the following equation:

$$R_{\mathrm{t}}={\displaystyle \frac{T_{\mathrm{e}, \mathrm{ave}}-T_{\mathrm{c}, \mathrm{ave}}}{Q_{\mathrm{p}}}}$$

(2)

where T_{e, ave} and T_{c, ave} are the time-average value of the wall temperature at the evaporator and condenser sections under steady state for about 10 min, respectively.

In addition, the variation of heat transfer ability of the SFTS during start-up was estimated by an equivalent thermal conductivity Kt, defined in Equation (3).

$$Kt={\displaystyle \frac{Q_{\mathrm{p}}{L}_{eff}}{A(T_{\mathrm{e}}-T_{\mathrm{c}})}}$$

(3)

where L_eff and A are the effective length and the cross-sectional area of the evaporator of the SFTS, respectively. T_e and T_c are the average wall temperature at the evaporator and condenser, respectively. In this work, L_eff is the total length of the SFTS, and A is the cross-sectional area of evaporator [22].

The relative errors of heat transfer power and thermal resistance can be evaluated by

$$E(Q_{\mathrm{p}})={\displaystyle \frac{C_P\rho \sqrt{v_{\mathrm{cooling}}^2{\mathrm{d}}^2(\Delta T_{\mathrm{out}-\mathrm{in}})+\Delta T_{\mathrm{out}-\mathrm{in}}^2{\mathrm{d}}^2{v}_{\mathrm{cooling}}}}{C_P\rho v_{\mathrm{cooling}}\Delta T_{\mathrm{out}-\mathrm{in}}}}=\sqrt{{\displaystyle \frac{{\mathrm{d}}^2(\Delta T_{\mathrm{out}-\mathrm{in}})}{\Delta T_{\mathrm{out}-\mathrm{in}}^2}}+{\displaystyle \frac{{\mathrm{d}}^2{v}_{\mathrm{cooling}}}{v_{\mathrm{cooling}}^2}}}$$

(4)

$$E(R)=\sqrt{{\displaystyle \frac{{\mathrm{d}}^2(\Delta T)}{\Delta T^2}}+{\displaystyle \frac{{\mathrm{d}}^2{Q}_{\mathrm{p}}}{{Q_{\mathrm{p}}}^2}}}$$

(5)

where ${\mathrm{d}(∆T}_{\mathrm{out}–\mathrm{in}})$ denotes the measurement error of temperature difference between the outlet and the inlet cooling water, ${\mathrm{d}v}_{\mathrm{cooling}}$ is the measurement error of the flow rate of the cooling water, $\mathrm{d}\left(∆T\right)$ represents the measurement error of the temperature difference between the wall surfaces at different sections and $\mathrm{d}{Q}_{\mathrm{p}}$ is the absolute error of Q_p. Under the given testing conditions, the estimated values of E(Q_p) and E(R) are less than 5% and 10%, respectively.

2.3. Heat Transfer Characteristics of the SFTS

To reveal the steady-state heat transfer characteristics of the STFS, Figure 3 shows the variations of T_e, T_c and Q_p during the test with the condenser inclination angle set to 0° and a v_cooling of 900 mL/min. The initial temperature of the cooling water fed to the shell heat exchanger at the condenser was 5.5 °C, and it was incrementally increased by 1 °C after a semi-steady state was reached. As Figure 3 shows, T_e, T_c and Q_p all exhibited a sharp reduction at the initial stage of start-up, but they leveled out soon. When T_cooling was 5.5 °C, the practical heat transfer power of the SFTS at the semi-steady state was approximately 230 W, with T_e and T_c reaching nearly 13.7 and 11.5 °C, respectively. As T_cooling increased, Q_p showed a significant reduction. It decreased to as low as 90 W when T_cooling was 10.5 °C. However, both T_e and T_c increased with the increasing T_cooling. T_c increased by 1.2 °C when T_cooling was increased to 10.5 °C, while T_e merely increased by 0.4 °C. It implied that the difference between T_e and T_c also reduced with an increased T_cooling. The above discussion revealed that the heat transfer power of the SFTS and the temperature difference between the evaporator and condenser T_e–T_c at a steady state are linearly related. Moreover, the slight temperature increase of T_e observed in the test suggested that the heat extraction rate from underground was below the maximum heat recovery rate around the borehole.

Previous studies have demonstrated that the heat transfer limit of a geothermal thermosyphon decreases significantly when its condenser is bent away from the vertical position [20]. However, the heat transfer power of the SFTS in this work remained well below its heat transfer limit. Consequently, the effects of inclination angle on the heat transfer performance at the steady state was investigated. Figure 4 shows the experimental results of the total thermal resistance R_t under varying T_cooling and inclination angles. It is shown in Figure 4 that the increase of T_cooling increased R_t extremely, implying a degradation in the overall heat transfer performance of the SFTS with an increase in T_cooling. Furthermore, as Figure 4 shows, the inclination angle had an insignificant influence on R_t at a lower T_cooling. When T_cooling was below 6.5 °C, R_t ranged from 0.02 to 0.03 K/W, with the lowest R_t observed at an inclination angle of 60°. However, as the inclination angle exceeded 60°, R_t increased, likely due to an increase in the shearing force of the two-phase flow at higher heat transfer power. Additionally, as T_cooling increased above 6.5 °C, the impact of inclination angle on R_t became more pronounced. At a T_cooling of 10.5 °C, the lowest R_t was approximately 0.031 K/W at an inclination angle of 90°, while the highest R_t was approximately 0.057 K/W at an inclination angle of 20°.

It follows that the arrangement of the condenser of an SFTS in practical applications with a lower temperature can be flexible, such as vertical placement for spacing heating and approximate horizontal placement for de-icing, anti-freezing [9] and lawn heating [23,24].

To further reveal the heat transfer characteristics of the SFTS during field tests, Figure 5 shows the variations of outer wall temperatures (T_e and T_c), Q_p and Kt during the start-up when the inclination angle was 0°, v_cooling was controlled at 900 mL/min and T_cooling at 5.5 °C. As Figure 5 shows, when the cooling water was fed into the shell heat exchanger at the condenser, a rapid thermal response of the SFTS was observed. At the initial stage of start-up, the temperatures at the condenser, adiabatic section and upper part of the evaporator (T₅–T₁₀) dropped sharply, while the low part of the evaporator (T₁–T₄) experienced only minor fluctuations, as seen in Figure 5a. Since the temperature variations leveled out soon, the overall temperature drop at the evaporator was minimal, but the temperature drop at the condenser was relatively high, as shown in Figure 5b. Initially, the heat transfer power of the SFTS was exceptionally high, peaking at approximately 550 W, as shown in Figure 5c. But it rapidly decreased to below 300 W in 150 s and then became stable. Since the SFTS was not normally operated at the initial stage of start-up and the geothermal energy could not be extracted via SFTS effectively, it appears that the high heat transfer power at that time was physically unreasonable. However, this phenomenon can be explained by the heat release from the heat capacity of the SFTS shell and the latent heat in vapor. Due to the high heat transfer power, the equivalent thermal conductivity Kt was also nominally much higher, as shown in Figure 5d. But after 150 s, when the phase change and flow circulation of the working fluid was formed in the SFTS chamber, Kt leveled out rapidly.

In addition, a comparison of Figure 5b,d reveals that when Kt began to become stable at about 150 s, the temperature difference between T_e and T_c also nearly reached a turning point. This observation suggested that the operation state of the SFTS during start-up could be roughly judged by ΔT_e–c, namely T_e–T_c, which can serve as a valuable parameter for modeling the dynamic heat transfer process in the SFTS.

3. CFD Modeling Method and Validation

3.1. Treatment Method of Equivalent Thermal Conductivity in an SFTS

The above discussion of heat transfer characteristics of the SFTS suggested that the heat transfer power of the SFTS Q_p and the temperature difference between the evaporator and condenser ΔT_e–c were linearly dependent. In addition, the operation state of the SFTS during the start-up, that is, whether or not it reached an effective heat transfer state, could be judged by the variation of ΔT_e–c. Accordingly, a simple CFD modeling method was developed in this work for the dynamic heat transfer simulation of an SFTS. In this modeling method, the complicated phase change and fluid flow in the SFTS was neglected, and the SFTS was treated as a homogeneous thermal superconductor just as the previous work [21]. However, the equivalent thermal conductivity of the SFTS was treated as a function of ΔT_e–c. When the SFTS reached an effective operating state, namely ΔT_e–c was higher than a specific value ΔT₀, the relationship of Kt and ΔT_e–c was determined by a fitting equation deduced from the field test results at the steady state, as shown in Figure 6. When ΔT_e–c was lower than ΔT₀, the SFTS did not work well, and thus Kt was given a low but constant value. To make sure the numerical solution was stable and converged, when ΔT_e–c was below ΔT₀, the constant value of Kt was also calculated by the fitting equation, but ΔT_e–c was replaced by ΔT₀. In addition, as the diameter of the evaporator differed from that of the condenser, and the inner cross-sectional area of the evaporator varied with the axis, an empirical coefficient λ was introduced to compensate for prediction errors. And then, the treatment of equivalent thermal conductivity of an SFTS was obtained, as shown in Equation (6).

$$K{t}_{\mathrm{CFD}}=\left\{\begin{array}{l}\lambda ·(1490196.04\Delta T_{\mathrm{e}-\mathrm{c}}-5602654.2),\Delta T_{\mathrm{e}-\mathrm{c}}> \Delta T_0\\ \lambda ·(1490196.04\Delta T_0-5602654.2),\Delta T_{\mathrm{e}-\mathrm{c}}\le \Delta T_0\end{array}\right.$$

(6)

where ΔT₀ is 4 °C in this work according to the variation of ΔT_e–c in Figure 5b.

3.2. CFD Modeling and Solution

3.2.1. Geometry and Meshing

We built a full-scale 3D computational geometry for the CFD simulation of the SFTS using the Design Modeler of ANSYS 19.0, as shown in Figure 7. The underground soil zone was a cylinder with a radius of 6 m and a depth of 30 m. The soil in the near-surface zone was treated as a homogeneous solid body, while the soil zone at depths of more than 2 m was assumed as homogeneous porous media to consider the influence of groundwater advection. The full-scale SFTS was also treated as a solid body coupled with the soil zones through the boundary conditions of coupled interfaces.

The discretization of the computational zone was performed by Meshing of ANSYS 19.0. After several iterations of debugging, a meshing structure was developed employing the tetrahedron method and the patch-conforming algorithm, with a maximum element size of 0.5 m. This configuration, shown in Figure 7c, resulted in a total of 7,148,104 elements, achieving an optimal balance between the computational cost and prediction accuracy.

3.2.2. Governing Equations

In the near-surface zone, merely the energy equation was solved to predict the temperature distribution. As the porous media model was applied in the soil zone at depths of more than 2 m, the corresponding governing equations are listed as follows [25]:

Continuity equation:

$${\displaystyle \frac{\partial }{\partial t}}\left(\epsilon {\rho }_{\mathrm{w}}\right)+\nabla ({\rho }_{\mathrm{w}}\overrightarrow{u})=0$$

(7)

where ε is soil porosity, ρ_w is the density of groundwater and $\overrightarrow{u}$ is the fluid velocity in soil.

Momentum equation:

$${\displaystyle \frac{\partial }{\partial t}}\left({\rho }_{\mathrm{w}}\overrightarrow{u}\right)+\nabla \cdot \left({\rho }_{\mathrm{w}}\overrightarrow{u}\overrightarrow{u}\right)=-\nabla p+\nabla \cdot \mu (\nabla \overrightarrow{u})+S_{\mathrm{i}}$$

(8)

where p represents the local pressure, μ is the dynamic viscosity and S_i is the source term.

In the porous media model, S_i is defined as [26]

$$S_{\mathrm{i}}=-\left({\displaystyle \frac{\mu }{\alpha }}{\overrightarrow{u}}_{\mathrm{i}}+C_2{\displaystyle \frac{1}{2}}\rho \left|u\right|{\overrightarrow{u}}_{\mathrm{i}}\right)$$

(9)

where α and C₂ are the permeability and inertial resistance factor, respectively. As the groundwater advection is of significantly low velocity, typically in laminar flow, C₂ can be considered to be zero here. Ignoring the convective acceleration and diffusion, the porous medium model is then reduced to Darcy’s law, which states that the pressure drop is proportional to velocity [15]

$$\nabla p=-{\displaystyle \frac{\mu }{\alpha }}\overrightarrow{u}$$

(10)

where the permeability α can be determined using the following expression:

$$\alpha ={\displaystyle \frac{{D_p}^2}{150}}{\displaystyle \frac{{\epsilon }^3}{{(1-\epsilon )}^2}}$$

(11)

where D_p is the mean particle diameter.

Energy equation:

$${\rho }_{\mathrm{s}}{c}_{\mathrm{s}}{\displaystyle \frac{\partial T_{\mathrm{s}}}{\partial t}}+{\rho }_{\mathrm{w}}{c}_{\mathrm{w}}\overrightarrow{u}\nabla T_{\mathrm{s}}=\nabla (k_{\mathrm{s}}\nabla T_{\mathrm{s}})$$

(12)

where ρ_sc_s is the effective volume-specific heat capacity of the soil, T_s is the temperature of soil, k_s represents the effective thermal conductivity of soil and ρ_wc_w is the volume-specific heat capacity of groundwater.

3.2.3. Boundary Conditions and Solution Strategy

The CFD modeling of the SFTS for shallow geothermal energy extraction in this work was carried out by ANSYS Fluent 19.0. The boundary conditions for the computational domain are illustrated in Figure 8. The inlet of groundwater was set to a velocity inlet boundary condition with the flow direction parallel to the z-axis, and the outlet of groundwater was set to an outflow. The temperature and advection of groundwater were set to constant values, 290.45 K and 0.75 m/day, respectively, according to the field tests [15]. The wall boundary condition was applied to all the other surfaces, while the thermal conditions were varied. For the bottom surface of the computational zone, the thermally adiabatic condition was applied; namely, the heat flux was zero. For the top surface of the computational zone, heat convection was imposed to consider the influence of air convection, with the heat transfer coefficient h_top and the free stream temperature as 13.3 W/(m²·K) and 284.95 K, respectively. The value of h_top was determined using Equation (13), an empirical correlation for calculating the total heat transfer coefficient on top surfaces [27], assuming an air velocity of 2 m/s. The free stream temperature corresponds to the air temperature measured in the testing area on the test day.

$$h_{\mathrm{top}}=5.7+3.8u_a$$

(13)

The local temperature on the side surface of the near-surface zone varied linearly with depths. The outer wall of condenser of the SFTS employed the heat convection boundary condition as well. The heat transfer coefficient at the condenser was determined using Equation (14), with the free stream temperature T_∞ corresponding to the average temperature of the cooling water at the inlet and outlet.

$$h_{\mathrm{c}}={\displaystyle \frac{Q}{A(T_{\mathrm{c}, \mathrm{ave}}-T_{\infty })}}$$

(14)

The modification of equivalent thermal conductivity in the SFTS with varying temperatures at the condenser and evaporator was performed by the User-defined Function (UDF). The corresponding UDF code was attached as an Supplementary Materials. The interface between the SFTS and the soil and among layers of soil was set to coupling heat transfer. The soil and groundwater parameters used for CFD simulation in this work are listed in

Table 2. Parameters used in CFD simulation.

	Item	Unit	Value
Porous media (soil)	Initial temperature	K	290.45
	Thermal conductivity	W/(m K)	2.1
	Specific heat capacity	J/(kg K)	1300
	Density	kg/m3	2100
	Porosity		20%
	Mean particle diameter	mm	0.01
Groundwater	Density	kg/m3	1000
	Conductivity	W/m·K	0.6
	Specific heat capacity	J/(kg K)	4200
	Flow velocity	m/day	0.75
	Inlet temperature	K	290.45

.

The pressure-based solver was used with the SIMPLE scheme for pressure–velocity coupling. For spatial discretization, the least squares cell-based scheme and the second-order scheme were employed for the gradient and pressure, respectively. While both the momentum and energy equations used the second-order upwind, each step was considered converged when the residuals of variables reduced to 10−3 for the continuity and momentum equations and 10−6 for the energy equation.

Prior to the transient solution, a steady-state solution was performed to obtain the initial conditions with the condenser of the SFTS specified as the thermally adiabatic condition. When the transient solution was carried out, the heat-convection condition of water cooling was imposed at the condenser. Initially, the transient simulation was calculated for 1 h with a fixed time step size of 10 s, corresponding to the conditions under which the thermostatic bath was controlled at 5.5 °C. Subsequently, just as in the field test, the cooling water temperature at the condenser varied from 5.5 to 10.5 °C in increments of 1 °C, and the corresponding heat transfer coefficient was modified in accordance with the test results. The modification of cooling conditions was merely performed when a semi-steady state was reached.

3.3. Modeling Validation

Three simulation cases were calculated with the empirical coefficient λ at 1, 2 and 3, respectively, to validate the model for both the start-up and steady state conditions.

Figure 9 shows the comparison between the experimental and simulation results at a steady state with T_cooling varying from 5.5 to 10.5 °C. The variation of practical heat transfer power Q_p is plotted in Figure 9a. As Figure 9a shows, the predicted trend of Q_p with T_cooling for different λ values closely aligned with the experimental results. However, Q_p was overpredicted when λ was set to 3. For T_cooling in the range of 5.5 to 7.5 °C, the CFD predictions of Q_p were in good agreement with the experimental data when λ was specified as 2, while, for a higher T_cooling, λ equal to 1 showed better prediction accuracy. In addition, Figure 9b compares the experimental and simulation results for the condenser temperature T_c. The variation trend of T_c with T_cooling was opposite to that of Q_p, but the prediction accuracy for various λ values was similar to that of Q_p. In most shallow-geothermal applications, an SFTS is operated at a low temperature at the condenser, typically near 5 °C for space heating and even lower for pavement de-icing. It is thus recommended that the empirical coefficient λ take 2.

The experimental and simulated variations of Q_p and T_c during the start-up are shown in Figure 10 under the conditions where T_cooling was controlled at 5.5 °C. As Figure 10 shows, in the initial 250 s, the variations of Q_p and T_c for various values of λ were not well consistent with the experimental results. This is due to the fact that the simplification assumptions made for the SFTS neglected the inner phase change and flow circulation dynamics, treating the equivalent thermal conductivity in the SFTS as a constant. However, when the simulation time was longer than 250 s, the predictions of Q_p and T_c showed good agreement with the experimental variations for λ, which was equal to 2. These findings indicate that the dynamic heat transfer modeling method developed for the SFTS in shallow geothermal applications is reliable and effective for capturing steady-state behavior.

4. Conclusions and Future Work

In this work, field tests of an SFTS on heat transfer characteristics were carried out. The effects of the cooling water temperature and inclination angle of the condenser on the start-up characteristics and steady-state heat transfer power were investigated. The results showed that the heat transfer power of the SFTS and the temperature difference T_e–T_c were linearly dependent at the steady state. The dependence of the inclination angle on heat transfer performance of the SFTS was negligible with a lower temperature at the condenser.

In addition, according to the field test results, a dynamic heat transfer modeling method of the SFTS based on the ETC model was proposed. A full-scale 3D CFD model for geothermal extraction via the SFTS was developed considering weather conditions and groundwater advection. Modeling validation was performed for both the start-up and steady state. It showed that the simulation results were in well agreement with the field test results when the empirical coefficient in the ETC model was specified as 2. It is proven that the modeling method proposed in this work can be further used for the overall simulation of a geothermal system with an integration of thermosyphons.

Future work will focus on applying the proposed model, which treats the SFTS as a homogeneous thermal superconductor, to predict the performance of geothermal systems designed for pavement ice and snow melting. Particular attention will be given to studying the effects of air conditions and the layout parameters of SFTS arrays.