FEFLOW | Tackling Heat Transport Modeling

MODELING TASK DESCRIPTION

Besides groundwater flow calculations, FEFLOW is a suitable modeling tool for any type of heat transport simulations. It is used for the design of open and closed-loop heating systems, different types of thermal energy storages (ATES, BTES, MTES) or entire BHE arrays with hundreds of borehole heat exchangers (BHEs). Heat transport simulations, however, are numerically complex and require specific settings, mesh configurations and boundary conditions to achieve high-performance simulations and reliable results. The following article provides some useful hints and tricks. 

SOLUTION

Step #1: Ideal element size
A successful heat transport model starts with an appropriate model discretization. Most models contain linear, subvertical heat extraction/injection wells or borehole heat exchangers. Both types require a specific mesh configuration honoring the concept of a virtual radius r_virtual (Fig. 1). This radius better represents anomalies around boreholes with radius r_b, if boundary conditions are applied to a series of singular vertical nodes (mesh singularity). The aim is to design the mesh in a way that the virtual borehole radius r_virtual best matches the physical radius r_b. In FEFLOW this concept is called ideal element size and is implemented with direct estimation of the nodal distance Δ. See Fig. 1.

Fig. 1. – Spatial discretization (n=6) around a BHE ‘well’ node.  

According to the above equation the nodal distance Δ is dependent on the virtual radius r_virtual and the number of nodes n around the well. Forcing r_virtual to be equal to r_b, the user can apply the following equation to calculate the ideal element size based on the number of adjacent nodes n є {4,6,8}. 

FEFLOW features the implementation of ideal element size nodes creating satellite points around each BHE well. After selecting all BHE nodes with equal borehole radius r_b, satellite points can be created with a specific button in the Mesh Editor toolbar, see Fig. 2. 

Step #2: Appropriate numerical settings
Transport simulations are more sensitive to numerical instabilities compared to flow simulations. The following Problem Settings (Fig. 3) apply to all kind of heat transport simulations and are recommended as initial settings. Advanced users can of course alter these settings during model calibration, but the suggested configuration has proved optimal results in various thermal applications. 

a) use First-order accurate (FE/BE) predictor-corrector scheme 

b) define a small error tolerance of 10e-04 (or even smaller)

c) choose the Maximum error norm instead of the Euclidean (default)

d) reduce the maximum number of iterations to 1

The latter point (d) reduces the computation time significantly. Heat transport simulations of BHE arrays can be computationally costly and especially if model convergence is not secured, the user benefits from shorter computation cycles per time-step.

If the model still suffers long computation times or missing convergence the user can apply stricter solver settings or try any of the parallelized solvers (e.g. SAMG, PARDISO). 

Step #3: Types of thermal boundary conditions

The energy balance of a geothermal model of the shallow or deep earth consists of multiple heat flows which all can be considered defining specific boundary conditions.

Advective heat
Groundwater flow or recharge create an advective heat flow which is a bulk transport of thermal energy contained in a moving fluid. It is typically implemented as a constant (or transient) temperature BC representing an average (or time-variant) groundwater temperature. It is typically assigned to all inflow nodes and/or to the top slice for groundwater recharge, respectively.
 

Convective heat
In arid climates, the convective heat exchange between air and ground surface may be of higher importance than the heat from recharge. Such heat exchange is typically implemented with a temperature BC of the ambient air temperature. Similarly, the thermal input from solar radiation can be considered as an additional heat flux via the top slice. It is recommended to use a Cauchy-type boundary condition defining a reference temperature and heat transfer-rates. 

Geothermal heat flux
For models with deep boreholes of several hundreds of meters, the geothermal heat flux from the earth itself might be considered. It is a constant heat flux from the bottom of the model which is global yet showing regional changes.  See Fig. 4. 

In FEFLOW, geothermal heat fluxes are realized with a Neumann-type flux boundary condition given in energy unit per square meter and time [J/m²/d] or power per square meter [e.g. mW/m²]. Please be aware of the negative sign convention for boundary conditions.

Consideration of the geothermal gradient
In deep geothermal applications, in addition, increasing temperatures with depth can play an important role for the heat extraction potential. In such cases, the geothermal gradient should be implemented both as initial condition (Fig. 5) and as boundary condition for the advective heat flux. 
 

The easiest way to assign vertically variable temperatures is to use the expression editor (Fig. 6). Applying the following formula calculates the depth-dependent temperature based on a minimum and maximum temperature at the upper and lower elevation of the profile while assuming a linear temperature gradient. 

T = temperature in [K]
z = elevation in [m]

The temperature initial condition is assigned to all nodes of the model domain, selecting the entire domain and applying the above expression. Boundary conditions are applied only to the inflow nodes using the process variable ‘Temperature’. 

Please remember to apply a geothermal gradient together with a Dirichlet type BC on the uppermost slice and a Neumann type flux BC at the bottom slice to account for the geothermal heat flux (see above).