COMSOL 6.4 - Borehole Heat Exchanger with Pipe Flow

Renewable energies are a growing industry and the geothermal energy branch is a hot topic of active research. Over the past few decades, different techniques were established to extract geothermal heat from shallow to deep subsurface levels. The closed-loop borehole heat exchanger (BHE) is a standard approach for lower- and mid-depth applications.

This model shows how to compute an array of borehole heat exchangers (BHEs) for shallow geothermal energy production. The BHEs are simplified as lines using the Nonisothermal Pipe Flow interface. The array is embedded into a layered subsurface model with groundwater flow in one of the layers.

Model Definition

The following model solves for the heat transport around a shallow geothermal installation embedded in a geological domain. The domain is separated into various parts, representing different geological layers with their particular properties. The thermal influence of seasonal temperature changes at the surface is taken into account using the COMSOL Multiphysics internal ambient weather database (in this case meteorological weather data from Berlin Tempelhof, Germany).

The geometry contains a three-by-three borehole heat exchanger (BHE) array that is 150 m deep and located in layered bedrock. Between 10 m and 50 m below the surface is an aquifer where groundwater flow occurs, causing horizontal convective heat transport. The BHEs are approximated by lines. Water of 10°C is pumped into the pipes at a volumetric flow rate of 6 l/min, corresponding to 10−4 m3/s. Figure 1 shows the geometry as implemented in the model.

Figure 1: Geometry of the BHE model. The boreholes whose temperatures have been further investigated and displayed in the Results section are marked in blue. The black arrows illustrate the aquifer flow field.

Although the boreholes look just like lines, each borehole is modeled as a bent pipe with inlet and outlet at the top.

The heat transfer in the whole model domain is computed using the Heat Transfer in Porous Media interface. Within the aquifer layer, the aquifer flow velocity is implemented as a convective velocity field u in the heat equation

(1)

The heat transfer between the BHEs and the surrounding bedrock is modeled using the Pipe Wall Heat Transfer multiphysics coupling node.

Results and Discussion

Figure 2: A multislice plot of the temperature after 1 year of simulated time. The borehole temperatures are displayed in a line plot and the arrow field shows the aquifer flow velocity.

Figure 2 and Figure 3 are examples of predefined ways to display the 3D temperature distribution in COMSOL Multiphysics. The top surface temperature changes seasonally between about 0°C and 22°C according to the meteorological weather data from Berlin Tempelhof, Germany, which leads to a layered temperature structure near the surface as shown in Figure 2. Both figures are dominated by the initial temperature profile, which shows a temperature increase with increasing depth. However, the cooling effect of the BHEs and the elongated heat pattern in the aquifer layer are clearly visible.

Figure 3: Isothermal layers around the BHEs after 1 year of simulated time.

Figure 4: Temperature profiles at the three different borehole positions relative to the aquifer flow field.

Shown in Figure 4 are the pipe-surrounding temperatures of the three BHEs in the middle of the array (as marked in Figure 1). The temperatures in the parts with downward flowing water are lower than in the parts with upward flowing water, as expected due to the heat exchange with the surrounding bedrock. The temperature gradient is much smaller than the initial one because the cooling effect of the BHEs is larger for larger temperature differences and varies slightly with different thermal properties of the different porous material layers.

Due to the thermal interaction between the heat sinks, the temperature of the middle BHE relative to the aquifer flow field (green line) is lower than that of the other two in the pipe sections with upward flow. An effect of the aquifer can hardly be seen in the temperature profiles.

Figure 5: Total line heat source integrated over all BHEs as a function of time.

Figure 5 shows the total line heat source with time, integrated over all BHEs. It decreases fast in the beginning of the simulation and more slowly when the heat transfer between pipes and porous media reaches quasi-equilibrium. Finally it seems to level out at about −3000 W.

Notes About the COMSOL Implementation

The boreholes are modeled as lines using the Nonisothermal Pipe Flow interface. The heat transfer between the pipes and the surrounding soil is modeled using the multiphysics coupling node. This node can only be activated on boundaries where a node is active. The node allows to define multiple wall layers — each defined in a subnode — with different material properties, which is used here to represent both the pipe walls and the grouting between the pipes and the borehole walls.

Reference

1. www.comsol.com/blogs/modeling-geothermal-processes-comsol-software

Pipe_Flow_Module/Heat_Transfer/borehole_heat_exchanger_pipe_flow

Modeling Instructions

From the menu, choose .

New

In the window, click

Model Wizard

In the window, click

In the tree, select > > .

Click .

In the tree, select > > .

Click .

Click

In the tree, select > .

Click

Geometry 1

The geometry sequence can be imported from an external file. As the pipe curvature at the bottom of the boreholes is very small compared to the borehole depth and the array size, the checkbox is cleared in the section of the node.

In the toolbar, click and choose .

Browse to the model’s Application Libraries folder and double-click the file borehole_heat_exchanger_pipe_flow_geom_sequence.mph.

In the toolbar, click

Click the

button in the toolbar.

In the window, under click .

In the window for , locate the section.

Clear the checkbox.

The geometry sequence is parametric. Add a few more parameters used setting up the physics.

Global Definitions

Parameters 1

In the window, under click .

In the window for , locate the section.

In the table, enter the following settings:


Name	Expression	Value	Description
pth	3[mm]	0.003 m	Pipe wall thickness
u_aquifer	40[m/a]	1.2675E-6 m/s	Groundwater flow velocity
vfr	6[dm^3/min]	1E-4 m³/s	Volumetric flow rate at pipe inlet

Materials

Define the different geological layers with their particular properties and use the selections provided by the feature in the geometry sequence.

Add Material from Library

In the toolbar, click

and choose .

Add Material

Go to the window.

In the tree, select > .

Click the button in the window toolbar.

In the toolbar, click

to close the window.

Materials

Water, liquid (mat1)

In the window for , locate the section.

From the list, choose .

Porous Material: Holocene Sediments

In the window, right-click and choose > .

In the window for , type Porous Material: Holocene Sediments in the text field.

Locate the section. From the list, choose .

Locate the section. Click

Fluid 1 (pmat1.fluid1)

In the window, click .

In the window for , locate the section.

From the list, choose .

Solid 1 (pmat1.solid1)

In the window, click .

In the window for , locate the section.

In the θs text field, type 0.7.

Locate the section. In the table, enter the following settings:


Property	Variable	Value	Unit	Property group
Density	rho	1600	kg/m³	Basic
Heat capacity at constant pressure	Cp	2000	J/(kg·K)	Basic
Thermal conductivity	k_iso ; kii = k_iso, kij = 0	1.2	W/(m·K)	Basic

Porous Material: Pleistocene Sands

In the window, under right-click and choose > .

In the window for , type Porous Material: Pleistocene Sands in the text field.

Locate the section. From the list, choose .

Locate the section. Click

Fluid 1 (pmat2.fluid1)

In the window, click .

In the window for , locate the section.

From the list, choose .

Solid 1 (pmat2.solid1)

In the window, click .

In the window for , locate the section.

In the θs text field, type 0.75.

Locate the section. In the table, enter the following settings:


Property	Variable	Value	Unit	Property group
Density	rho	1700	kg/m³	Basic
Heat capacity at constant pressure	Cp	800	J/(kg·K)	Basic
Thermal conductivity	k_iso ; kii = k_iso, kij = 0	1.3	W/(m·K)	Basic

Porous Material: Pleistocene Glacial Till

In the window, under right-click and choose > .

In the window for , locate the section.

From the list, choose .

In the text field, type Porous Material: Pleistocene Glacial Till.

Locate the section. Click

Fluid 1 (pmat3.fluid1)

In the window, click .

In the window for , locate the section.

From the list, choose .

Solid 1 (pmat3.solid1)

In the window, click .

In the window for , locate the section.

In the θs text field, type 0.85.

Locate the section. In the table, enter the following settings:


Property	Variable	Value	Unit	Property group
Density	rho	2100	kg/m³	Basic
Heat capacity at constant pressure	Cp	900	J/(kg·K)	Basic
Thermal conductivity	k_iso ; kii = k_iso, kij = 0	1.8	W/(m·K)	Basic

Porous Material: Tertiary Sands

In the window, under right-click and choose > .

In the window for , type Porous Material: Tertiary Sands in the text field.

Locate the section. From the list, choose .

Locate the section. Click

Fluid 1 (pmat4.fluid1)

In the window, click .

In the window for , locate the section.

From the list, choose .

Solid 1 (pmat4.solid1)

In the window, click .

In the window for , locate the section.

In the θs text field, type 0.8.

Locate the section. In the table, enter the following settings:


Property	Variable	Value	Unit	Property group
Density	rho	1900	kg/m³	Basic
Heat capacity at constant pressure	Cp	850	J/(kg·K)	Basic
Thermal conductivity	k_iso ; kii = k_iso, kij = 0	1.4	W/(m·K)	Basic

The node indicates that there are still missing material properties, which is the permeability required for Darcy’s Law. The permeability in the lower two layers is about two orders of magnitude smaller than in the upper two. Therefore add the permeability only for the upper layers.

Nonisothermal Pipe Flow (nipfl)

Now start to set up the physics. Start with the Nonisothermal Pipe Flow interface which should only be active on the pipe edges.

In the window, under click .

In the window for , locate the section.

From the list, choose .

Pipe Properties 1

In the window, under > click .

In the window for , locate the section.

From the list, choose .

In the di text field, type 2*(r_pipe-pth).

Temperature 1

In the window, click .

In the window for , locate the section.

In the Tin text field, type 283.15[K].

Initial Values 1

In the window, click .

In the window for , locate the section.

In the T text field, type 283.15[K].

In the p text field, type 1[atm].

Pressure 1

In the window, click .

In the window for , locate the section.

In the p0 text field, type 1[atm].

Inlet 1

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Locate the section. From the list, choose .

In the qv,0 text field, type vfr.

Heat Outflow 1

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Wall Heat Transfer 1

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

As the settings are influenced by the multiphysics coupling between heat transfer in pipes and heat transfer in the surrounding porous medium, add the corresponding multiphysics node before continuing.

Multiphysics

Pipe Wall Heat Transfer 1 (pwhtc1)

In the toolbar, click

and choose > .

Nonisothermal Pipe Flow (nipfl)

Wall Heat Transfer 1

In the window, under > click .

Wall Layer 1

In the toolbar, click

and choose .

In the window for , locate the section.

From the k list, choose .

In the text field, type 0.38[W/m/K].

From the Δw list, choose .

In the text field, type pth. Now you have added the thermal conductivity of the pipe material and the pipe wall thickness. In a second layer consider the grouting material which fills the space between pipes and borehole walls.

Wall Heat Transfer 1

In the window, click .

Wall Layer 2

In the toolbar, click