GHX
IDA ICE borehole model
Borehole model
Lars Eriksson and Patrik Skogqvist, EQUA Simulation AB, 2023
The model describes the heat exchange between the fluid in U-tubes, embedded in a filling material (“grout”) in each borehole, and the surrounding ground:
Fig. 1: Model ghx „ground heat exchanger“
A number of temperature fields are computed using finite difference technique and combined by superposition to generate the 3-dimensional field. For each borehole the following fields are computed:
- 1-dimensional heat transport equation in U-pipe liquid downward and upward with heat transfer to grout and ground.
- 1-dimensional heat equation in grout with heat transfer to liquid and ground.
- 2-dimensional heat equation in cylindrical coordinates around borehole with heat transfer to grout and liquid.
The boreholes utilize the same discretization in liquid, grout and surrounding ground but may have different inclination. This implies that the ground properties cannot vary with depth and that all boreholes have the same length.
In addition to this the undisturbed ground temperature field is needed, which is computed using the 1-dimensional heat equation with a heat transfer coefficient towards ambient temperature at the ground surface. A geothermal gradient can also be included in the undisturbed ground field.
Our model can utilize symmetry to speed up the computation. The degree of symmetry is governed by a user given parameter:
- mir = 0 : No symmetry utilized, boreholes may be placed in the whole x-y plane.
- mir = 1 : Mirror in x-axis, boreholes specified in half plane y ≥ 0.
- mir = 2 : Mirror in x- and y-axis, boreholes specified in quadrant y ≥ 0, x ≥0.
- mir = 3 : Mirror in x-, y-axis and 45°-axis, boreholes specified in sector 0 ≤ y ≤x
Description of the model for s Single Hole
The model equations for a single hole and its interaction with the ground is first described in this section. In
later sections is the interaction between the holes and the setup with groups of holes described.
Heat exchange fluid – filling material „grout“ – ground „earth“
The heat transfer from the brine into the filling material of the borehole is calculated, as well as the heat conduction from this filling material to the border of the borehole in the ground. The whole system ground-borehole-U-tube is divided into layers, in which temperature of borehole and fluid is assumed to be homogeneous:
Fig. 2: Heat exchange fluid – filling material („grout“) – ground material („earth“) in layer j of borehole i.
Fig. 3: Thermal masses (marked with points) and thermal connections (black lines) in borehole i.
The energy balance in the fluid is:
Tdi,j (°C) | Temperature of the down going fluid in node j of borehole i |
| Tdi,j-1 | = TIn,i for layer j=1 |
| Tui,j (°C) | Temperature of the up going fluid in node j of borehole i |
| Tui,j+1 | =Tdi,nzHole for node j = nzHole |
| mi (kg/s) | Mass flow of the fluid in borehole i |
| ρLiq (kg/m3) | Density of the fluid computed as a function of brine type, liquid freezing point and mean of incoming and leaving liquid temperatures |
| cp,Liq (J/(kg K)) | Specific heat capacity of the fluid computed as a function of brine type, liquid freezing point and mean of incoming and leaving liquid temperatures |
| KLiqGrout,i (W/K) | Heat transfer coefficient between fluid and filling material in borehole i |
| KLiqEarth,i (W/K) | Heat transfer coefficient between fluid and ground in borehole i |
| TGroutd,i,j (°C) | Grout temperature around up flow pipe(s) at layer j of borehole i |
| TGroutu,i,j (°C) | Grout temperature around down flow pipe(s) at layer j of borehole i |
| TReal,i,j (°C) | Ground temperature at the borehole i in node j |
| nzhole (- ) | Number of nodes in borehole |
The volume of each pipe cell is:
| nPipe (- ) | Number of U-tubes per borehole |
| rPipe (m) | Radius of tubes |
| zLay (m) | Thickness of layer |
The Heat transfer coefficient between fluid and filling material is:
| αi (W/(m2 K)) | Heat transfer coefficient from fluid to the inner surface of the U-tubes in borehole i |
| RPipeGrout (K/(W/m2)) | Surface specific heat resistance between fluid and filling material |
| RPipeEarth (K/(W/m2)) | Heat resistance between fluid and ground |
The heat transfer coefficient is:
| λLiq (W/(m K)) | Heat conductivity of the fluid |
| Nui (- ) | Nusselt number for borehole i |
The nusselt number is:
| Rei (- ) | Reynolds number for borehole i |
| NuLam,i (- ) | Nussel number for laminar flow |
| NuTurb,i (- ) | Nussel number for turbulent flow |
| NuMixed,i (- ) | Nussel number for mixed turbulent/laminar flow |
| a,b | Coefficients such that Nui becomes continuous |
| ci | Coefficients such that Nui has continuous derivatives at breaking points, Re = 2300 and Re = 3500. |
The Reynolds number is:
| vi (m/s) | Velocity of the fluid in borehole i |
| νLiq (m2/s) | Viscosity of the fluid |
The velocity of the fluid is:
The energy balance for the filling material is:
| Mcp,Grout 1 (J/K) | Absolute heat capacity of inner grout |
| Mcp,Grout 2 (J/K) | Absolute heat capacity of outer grout |
| KGroutEarth (W/K) | Heat conduction coefficient between grout and ground |
| TGrout,i,j (°C) | Temperature of outer grout in node j of borehole i |
| TReal u,i,j (°C) | Grout temperature around upflow pipe(s) at node j of borehole i |
| TReal d,i,j (°C) | Grout temperature around downflow pipe(s) at node j of borehole i |
| TReal,i,j (°C) | Ground temperature at the borehole i in node j |
The absolute heat capacities of the filling material is:
| ρGrout (kg/m3) | Density of the filling material |
| cp,Grout (J/(kg K)) | Mass specific heat capacity of the filling material |
| rHole (m) | Radius of the borehole |
| rPipe (m) | Radius of tubes |
| rGroutRing (m) | Radius of inner grout around tubes |
The heat conductivity coefficient from the core of the filling material to the border towards ground is:
| RGroutEarth (K/(W/m2) | Surface specific heat resistance between filling material and ground |
The heat conductivity coefficient between grout ring and ground is:
| RRingEarth (K/(W/m2) | Heat resistance between grout ring and ground |
The heat conductivity coefficient between grout ring and outer grout is:
| RGroutGrout (K/(W/m2) | Heat resistance between grout ring and grout. |
Pressure loss in the brine liquid loop
The pressure loss in the brine loop is calculated with the Darcy-Weisbach equation:
| Δp (Pa) | Difference in pressure |
| L (m) | Total length of pipe |
| F | the Darcy-Weisbach friction factor |
| ΣK | Sum of minor losses |
The Darcy-Weisbach friction factor is computed with the expanded form:
| ε (m) | Rougness height in pipe. |
Heat conductivity and storage in the ground
For the calculation of the heat conductivity and storage, the ground around a borehole is divided in a 2 dimensional grid. In axial direction, the grid corresponds to the layers of the borehole, except that it can be continued below. In radial direction, the ground is divided into ground rings.
Fig. 4: Computational mesh around borehole
The energy balance for the ground ring elements is:
| Mcp,Earth,j,k (J/K) | Absolute heat capacity associated to nodes (j, k) |
| DTi,j,k (°C) | Temperature change due to conditions in hole i at node (j, k) |
| KrEarth,k (W/K) | Heat conductivity coefficient in radial direction between radial node k and k+1 at axial nodes j |
| KzEarth,k (W/K) | Heat conductivity coefficient in axial direction between axial nodes j and j+1 at radial nodes k |
if k=1, the term with KrEarth,k-1 is replaced by 0 if above or below hole or else by
| nz1 | Number of nodes in extended computational domain above ground surface |
if k=nRing+1, the term with KrEarth,k is replaced by 0, i.e. insulation
| DTi,j-1,k | = 0 for j=1 |
| DTi,j+1,k | = 0 for j = nLayTot1 |
| TReal,i,j-nz1 (°C) | Ground temperature at the borehole i in node j |
| nLayTot1 | Number of nodes in extended computational domain |
The absolute heat capacity associated to a node is:
and
| rk (m) | = Radial distance to nodes k |
| Gridfac (-) | Enlargement factor for radial discretization |
| rk-1 | = rk for k = 1 |
| rk+1 | rk for k = nRing+1 |
| ρEarth (kg/m3) | Density of the ground |
| cp,Earth (J/(kg K)) | Specific heat capacity of the ground |
| zLay (m) | Thickness of layers |
| βk | Computed factor for symmetry of mass inside and outside of node |
| β0 | = 0 and βnRing+1 = 1 |
The heat conductivity coefficient in radial direction is:
| λLay (W/(m K)) | Heat conductivity of ground |
The heat conductivity coefficient in axial direction is:
Undisturbed ground field
The temperature field of the undisturbed ground (without borehole) is computed by discretizing the one dimensional transient heat equation in the same way as in the axial ground field around the boreholes. In that way we can reuse coefficients already defined:
| DTref,0 (°C) | Temperature change in undisturbed ground surface relative mean ground temperature |
| DTref,j (°C) | Temperature change in undisturbed ground node j relative mean ground temperature |
| DTref,nLayt+1 | = 0 i.e. Dirichlet boundary condition at bottom |
| ρsurf (kg/m3) | Density of surface layer |
| cp,surf (J/kg K) | Specific heat of surface layer |
| Lc (m) | Depth of surface layer (Length of casting) |
| TSurf (°C) | Ambient temperature |
| TMean(°C) | Mean ambient temperature used as mean ground temperature |
| λA (W/K) | Heat transfer at ground surface
|
| λSurf (W/( m2K)) | Heat transfer coefficient at ground surface |
A geothermal gradient may also be included giving the undisturbed ground field
| Tundist,j(°C) | Undisturbed ground temperature at vertical node j |
| geoGrad (°C/m) | Geothermal gradient |
Superposition of ground fields
The actual ground temperature at the boreholes, TReal,i,j, used in the equations above is computed using superposition.
| TReal,i,j (°C) | Actual ground temperature at borehole i and node j |
| DTi,j+nz1,1 (°C) | Temperature change due to conditions in hole i to node j |
| DTSumi,j (°C) | Superposition of all other borehole 2-D fields, their symmetry images, all borehole 2-D field images mirrored in ground surface and undisturbed ground temperature field |
With more than one borehole, the thermal behavior in the ground around each borehole is calculated as described for one single borehole. The influence of all boreholes to the ground temperature is summed up by superposition.
For this, the temperature of the undisturbed ground (without borehole) TRefj is calculated. Like this, the temperatures of the ground ring elements can be divided in:
| Ti,j,k (°C) | Temperature or the ground ring element k in layer j of borehole i |
| TRef,j (°C) | Temperature of the undisturbed ground in layer j |
| DTi,j,k (°C) | Ground temperature rise in ground ring element k in layer j of borehole i due to the same borehole i |
| DTSumi,j (°C) | Ground temperature rise in layer j at borehole i due to all other boreholes |
In the energy balance of the ground ring elements, this leads to:
This approximation can be done, as TRef,j and DTSumi,j do not vary significantly in time. The energy balance gets:
| DTi,j-1,k | = 0 for j = 1 |
| DTi,j+1,k | = 0 for j = nLayTot1 |
For calculation of TRef,j, a new energy balance is needed:
| TRef,j-1 | = TSurf for j = 1 |
| TRef,j+1 | = TDeep for j = nLayTot1 |
It remains the equation for DTSumi,j:
| Disti,h (m) | Distance between borehole i and borehole h |
Fig. 5: Place of borehole h in relation to borehole i.