PDF

Borehole Heat Exchanger
Introduction
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 line heat sinks with a uniform extraction rate. The array is embedded into a layered subsurface model with groundwater flow in one of the layers.
Model Definition
This model has been set up for comparison with the model using the Pipe Flow Module found in the application library under Pipe_Flow_Module/Heat_Transfer/borehole_heat_exchanger_pipe_flow and is also presented in a COMSOL Blog post (Ref. 1). While that model uses the Nonisothermal Pipe Flow interface, the current model shows an alternative approach using line heat sources.
The model-geometry and boundary-condition setups are similar to those of the model version with pipe flow (Figure 1). The main difference is that here the boreholes rather than the pipes are represented, and therefore only straight lines instead of u-shaped curves at the bottom are implemented.
The borehole heat exchangers (BHEs) are arranged in a three-by-three array that is 150 m deep and located in layered bedrock. Between 10 m and 50 m is an aquifer where groundwater flow occurs, causing horizontal convective heat transport.
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.
The heat transfer in the whole model domain is computed using the Heat Transfer in Porous Media interface. In addition, the fluid flow in porous media is calculated within the top layer and the aquifer layer using the Darcy’s Law interface.
The borehole heat exchangers are modeled as line heat sources which extract energy from the model. The heat rate is set to 3000 W, which corresponds to the value that is calculated in the model using the Pipe Flow Module.
Results and Discussion
Figure 2: Temperature distribution after 1 year of simulated time within and around the BHEs.
Figure 2 and Figure 3 are examples of predefined ways to display the 3D temperature distribution in COMSOL Multiphysics. Both figures are dominated by the initial temperature profile, which shows a temperature increase with increasing depth. The top surface temperature changes seasonally between about 0°C and 22°C according to the meteorological weather data from Berlin, Germany, which leads to a layered temperature structure near the surface as shown in Figure 2. Figure 3 shows how the heat pattern elongates in the direction of flow due to convection in the region of the aquifer.
Figure 3: Isothermal layers around the BHEs after 1 year of simulated time.
When comparing Figure 2 and Figure 3 with the corresponding figures of the model using the Pipe Flow module, there is hardly any difference in the global temperature field (Figure 3). In Figure 2, the temperature of the pipe-surrounding bedrock is shown while for the version using pipe flow, the temperatures within the pipes are displayed. Therefore, the influence of the external temperature at the top is visible in this model version. However, the pipe flow model version shows a stronger cooling effect near the boreholes in the lower half of the model and a smaller cooling effect closer to the surface for the model version using pipe flow, as the heat extraction rate in the latter case can vary with height whereas here it is provided as a constant.
Figure 4: Temperature profiles at the three different borehole positions relative to the aquifer flow field.
Figure 4 displays the borehole wall temperatures of the three BHEs in the middle of the array (as marked in Figure 1). Due to the thermal interaction between the heat sinks, the temperature of middle BHE relative to the aquifer flow field (green line) is lower than the other two. Only in the region of the aquifer, the BHE further downstream than the other two (red line) is additionally affected by the heat exchange occurring upstream of it due to convection, resulting in a lower temperature.
Compared with the temperature profile of the model version using pipe flow, the temperature profiles here are dominated by the initial temperature profile. A vertically constant heat extraction only shifts the temperature profile to lower temperatures rather than influencing the temperature gradient. While in the model version using the Pipe Flow Module, the effect of the aquifer on the temperature profiles is exceeded by the effect of vertically varying energy extraction, it is clearly visible here as described above.
Notes About the COMSOL Implementation
The boreholes are modeled as line heat sources. To make it available to users that just have the Subsurface Flow Module or Porous Media Flow Module licensed, the line heat source option does not contain the radius of the line heat source. The option to enter a line heat source radius is available with the Heat Transfer Module and does actually affect the temperatures around the boreholes quantitatively by shifting the profiles about 0.5 K toward higher values. However, as described when interpreting Figure 4, for a vertically constant line heat source, the temperature profiles differ significantly from the ones modeled with the Pipe Flow Module in general, so the qualitative results of both the model and the model comparison stays the same.
Reference
1. www.comsol.com/blogs/modeling-geothermal-processes-comsol-software
Application Library path: Subsurface_Flow_Module/Heat_Transfer/borehole_heat_exchanger
Modeling Instructions
From the File menu, choose New.
New
In the New window, click  Model Wizard.
Model Wizard
1
In the Model Wizard window, click  3D.
2
In the Select Physics tree, select Fluid Flow > Porous Media and Subsurface Flow > Darcy’s Law (dl).
3
Click Add.
4
In the Select Physics tree, select Heat Transfer > Porous Media > Heat Transfer in Porous Media (ht).
5
Click Add.
6
Click  Study.
7
In the Select Study tree, select General Studies > Time Dependent.
8
Geometry 1
1
In the Geometry toolbar, click Insert Sequence and choose Insert Sequence.
2
3
In the Geometry toolbar, click  Build All.
The geometry sequence is parametric. Add a few more parameters used setting up the physics.
Global Definitions
Parameters 1
1
In the Model Builder window, under Global Definitions click Parameters 1.
2
In the Settings window for Parameters, locate the Parameters section.
3
Materials
Now set up the materials for each layer and use the selections provided by the Block feature in the geometry sequence.
Add Material from Library
In the Home toolbar, click  Windows and choose Add Material from Library.
Add Material
1
Go to the Add Material window.
2
In the tree, select Built-in > Water, liquid.
3
Click the Add to Component button in the window toolbar.
4
In the Home toolbar, click  Add Material to close the Add Material window.
Materials
Porous Material: Holocene Sediments
1
In the Model Builder window, under Component 1 (comp1) right-click Materials and choose More Materials > Porous Material.
2
In the Settings window for Porous Material, type Porous Material: Holocene Sediments in the Label text field.
3
Locate the Geometric Entity Selection section. From the Selection list, choose Top (Block 1).
4
Locate the Phase-Specific Properties section. Click  Add Required Phase Nodes.
Fluid 1 (pmat1.fluid1)
1
In the Model Builder window, click Fluid 1 (pmat1.fluid1).
2
In the Settings window for Fluid, locate the Fluid Properties section.
3
From the Material list, choose Water, liquid (mat1).
Solid 1 (pmat1.solid1)
1
In the Model Builder window, click Solid 1 (pmat1.solid1).
2
In the Settings window for Solid, locate the Solid Properties section.
3
In the θs text field, type 0.7.
4
Locate the Material Contents section. In the table, enter the following settings:
Porous Material: Pleistocene Sands
1
In the Model Builder window, under Component 1 (comp1) right-click Materials and choose More Materials > Porous Material.
2
In the Settings window for Porous Material, type Porous Material: Pleistocene Sands in the Label text field.
3
Locate the Geometric Entity Selection section. From the Selection list, choose Aquifer (Block 1).
4
Locate the Phase-Specific Properties section. Click  Add Required Phase Nodes.
Fluid 1 (pmat2.fluid1)
1
In the Model Builder window, click Fluid 1 (pmat2.fluid1).
2
In the Settings window for Fluid, locate the Fluid Properties section.
3
From the Material list, choose Water, liquid (mat1).
Solid 1 (pmat2.solid1)
1
In the Model Builder window, click Solid 1 (pmat2.solid1).
2
In the Settings window for Solid, locate the Solid Properties section.
3
In the θs text field, type 0.75.
4
Locate the Material Contents section. In the table, enter the following settings:
Porous Material: Pleistocene Glacial Till
1
In the Model Builder window, under Component 1 (comp1) right-click Materials and choose More Materials > Porous Material.
2
In the Settings window for Porous Material, locate the Geometric Entity Selection section.
3
From the Selection list, choose Glacial Till (Block 1).
4
In the Label text field, type Porous Material: Pleistocene Glacial Till.
5
Locate the Phase-Specific Properties section. Click  Add Required Phase Nodes.
Fluid 1 (pmat3.fluid1)
1
In the Model Builder window, click Fluid 1 (pmat3.fluid1).
2
In the Settings window for Fluid, locate the Fluid Properties section.
3
From the Material list, choose Water, liquid (mat1).
Solid 1 (pmat3.solid1)
1
In the Model Builder window, click Solid 1 (pmat3.solid1).
2
In the Settings window for Solid, locate the Solid Properties section.
3
In the θs text field, type 0.85.
4
Locate the Material Contents section. In the table, enter the following settings:
Porous Material: Tertiary Sands
1
In the Model Builder window, under Component 1 (comp1) right-click Materials and choose More Materials > Porous Material.
2
In the Settings window for Porous Material, type Porous Material: Tertiary Sands in the Label text field.
3
Locate the Geometric Entity Selection section. From the Selection list, choose Core (Block 1).
4
Locate the Phase-Specific Properties section. Click  Add Required Phase Nodes.
Fluid 1 (pmat4.fluid1)
1
In the Model Builder window, click Fluid 1 (pmat4.fluid1).
2
In the Settings window for Fluid, locate the Fluid Properties section.
3
From the Material list, choose Water, liquid (mat1).
Solid 1 (pmat4.solid1)
1
In the Model Builder window, click Solid 1 (pmat4.solid1).
2
In the Settings window for Solid, locate the Solid Properties section.
3
In the θs text field, type 0.8.
4
Locate the Material Contents section. In the table, enter the following settings:
The Material 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.
Porous Material: Holocene Sediments (pmat1)
1
In the Model Builder window, under Component 1 (comp1) > Materials click Porous Material: Holocene Sediments (pmat1).
2
In the Settings window for Porous Material, locate the Homogenized Properties section.
3
Porous Material: Pleistocene Sands (pmat2)
1
In the Model Builder window, click Porous Material: Pleistocene Sands (pmat2).
2
In the Settings window for Porous Material, locate the Homogenized Properties section.
3
Continue with setting up the physics. Add the ambient properties to provide meteorological data for the upper boundary conditions.
Definitions (comp1)
Ambient Properties 1 (ampr1)
1
In the Physics toolbar, click  Shared Properties and choose Ambient Properties.
2
In the Settings window for Ambient Properties, locate the Ambient Settings section.
3
From the Ambient data list, choose Meteorological data (ASHRAE 2021).
4
Locate the Location section. Click Set Weather Station.
5
In the Weather Station dialog, select Europe > Germany > BERLIN TEMPELHOF (103840) in the tree.
6
Darcy’s Law (dl)
Next, set up the physics for the Darcy’s Law interface. It is sufficient to activate it only in the domains with significant flow velocities (that is, the aquifer layer and top layer).
1
In the Model Builder window, under Component 1 (comp1) click Darcy’s Law (dl).
2
In the Settings window for Darcy’s Law, locate the Domain Selection section.
3
Click  Clear Selection.
4
5
Locate the Gravity Effects section. Select the Include gravity checkbox.
Assume a constant groundwater flow velocity in the aquifer layer. This can be done by using the Contributing Velocity feature.
Porous Medium 1
In the Model Builder window, under Component 1 (comp1) > Darcy’s Law (dl) click Porous Medium 1.
Contributing Velocity 1
1
In the Physics toolbar, click  Attributes and choose Contributing Velocity.
2
In the Settings window for Contributing Velocity, locate the Domain Selection section.
3
Click  Clear Selection.
4
5
Locate the Contributing Velocity section. Specify the uctrb vector as
Initial Values 1
1
In the Model Builder window, under Component 1 (comp1) > Darcy’s Law (dl) click Initial Values 1.
2
In the Settings window for Initial Values, locate the Initial Values section.
3
Click the Hydraulic head button.
Hydraulic Head 1
1
In the Physics toolbar, click  Boundaries and choose Hydraulic Head.
2
In the Settings window for Hydraulic Head, locate the Boundary Selection section.
3
From the Selection list, choose Exterior.
Precipitation 1
1
In the Physics toolbar, click  Boundaries and choose Precipitation.
2
3
In the Settings window for Precipitation, locate the Precipitation section.
4
From the P0 list, choose Ambient precipitation rate (ampr1).
5
Definitions (comp1)
Before setting up the heat transfer in porous media, add an analytic function to describe the initial temperature profile in the bedrock.
Analytic 1 (an1)
1
In the Definitions toolbar, click  Analytic.
2
In the Settings window for Analytic, type T0 in the Function name text field.
3
Locate the Definition section. In the Expression text field, type 283.15+3[K]/100[m]*abs(z).
4
In the Arguments text field, type z.
5
Locate the Units section. In the Function text field, type K.
6
7
Locate the Plot Parameters section. In the table, enter the following settings:
8
Heat Transfer in Porous Media (ht)
Now set up the physics of the Heat Transfer in Porous Media interface.
Porous Matrix 1
1
In the Model Builder window, under Component 1 (comp1) > Heat Transfer in Porous Media (ht) > Porous Medium 1 click Porous Matrix 1.
2
In the Settings window for Porous Matrix, locate the Matrix Properties section.
3
From the Define list, choose Solid phase properties.
Initial Values 1
1
In the Model Builder window, under Component 1 (comp1) > Heat Transfer in Porous Media (ht) click Initial Values 1.
2
In the Settings window for Initial Values, locate the Initial Values section.
3
In the T text field, type T0(z).
Porous Medium 2
1
In the Physics toolbar, click  Domains and choose Porous Medium.
2
Fluid 1
1
In the Model Builder window, click Fluid 1.
2
In the Settings window for Fluid, locate the Heat Convection section.
3
From the u list, choose Total Darcy velocity field (dl/porous1).
Porous Matrix 1
1
In the Model Builder window, click Porous Matrix 1.
2
In the Settings window for Porous Matrix, locate the Matrix Properties section.
3
From the Define list, choose Solid phase properties.
Temperature 1
1
In the Physics toolbar, click  Boundaries and choose Temperature.
2
3
In the Settings window for Temperature, locate the Temperature section.
4
From the T0 list, choose Ambient temperature (ampr1).
Temperature 2
1
In the Physics toolbar, click  Boundaries and choose Temperature.
2
3
In the Settings window for Temperature, locate the Temperature section.
4
In the T0 text field, type T0(z).
Open Boundary 1
1
In the Physics toolbar, click  Boundaries and choose Open Boundary.
2
3
In the Settings window for Open Boundary, locate the Upstream Properties section.
4
In the Tustr text field, type T0(z).
Line Heat Source 1
1
In the Physics toolbar, click  Edges and choose Line Heat Source.
2
In the Settings window for Line Heat Source, locate the Line Heat Source section.
3
From the Heat source list, choose Heat rate.
4
In the Pl text field, type -Pext.
5
Locate the Edge Selection section. From the Selection list, choose Array 1.
Mesh 1
The next step is to set up the mesh. Choose a fine mesh around the boreholes to guarantee that the heat exchange is correctly simulated. A swept mesh in the vertical direction is more efficient than a triangular mesh.
Free Triangular 1
1
In the Mesh toolbar, click  More Generators and choose Free Triangular.
2
Size
1
In the Model Builder window, click Size.
2
In the Settings window for Size, locate the Element Size section.
3
From the Predefined list, choose Extra fine.
Size 1
1
In the Model Builder window, right-click Free Triangular 1 and choose Size.
2
In the Settings window for Size, locate the Geometric Entity Selection section.
3
From the Geometric entity level list, choose Point.
4
5
Locate the Element Size section. Click the Custom button.
6
Locate the Element Size Parameters section.
7
Select the Maximum element size checkbox. In the associated text field, type 0.02.
8
Click  Build Selected.
Swept 1
In the Mesh toolbar, click  Swept.
Distribution 1
1
Right-click Swept 1 and choose Distribution.
2
In the Settings window for Distribution, locate the Domain Selection section.
3
Click  Clear Selection.
4
5
Locate the Distribution section. In the Number of elements text field, type 15.
Distribution 2
1
In the Model Builder window, right-click Swept 1 and choose Distribution.
2
In the Settings window for Distribution, locate the Domain Selection section.
3
Click  Clear Selection.
4
5
Locate the Distribution section. In the Number of elements text field, type 20.
6
Click  Build All.
Study 1
1
In the Model Builder window, click Study 1.
2
In the Settings window for Study, locate the Study Settings section.
3
Clear the Generate default plots checkbox.
Step 1: Time Dependent
1
In the Model Builder window, under Study 1 click Step 1: Time Dependent.
2
In the Settings window for Time Dependent, locate the Study Settings section.
3
From the Time unit list, choose a.
4
In the Output times text field, type range(0,1/24,1).
Solution 1 (sol1)
1
In the Study toolbar, click  Show Default Solver.
2
In the Model Builder window, expand the Solution 1 (sol1) node, then click Time-Dependent Solver 1.
3
In the Settings window for Time-Dependent Solver, click to expand the Time Stepping section.
4
From the Steps taken by solver list, choose Strict to ensure that the changes of the ambient temperature are represented properly.
Step 1: Time Dependent
In the Study toolbar, click  Compute.
Results
To create Figure 2 and Figure 3 open the Result Templates.
Result Templates
1
In the Results toolbar, click  Result Templates to open the Result Templates window.
2
Go to the Result Templates window.
3
In the tree, select Study 1/Solution 1 (sol1) > Heat Transfer in Porous Media > Temperature, Multislice (ht) and Study 1/Solution 1 (sol1) > Heat Transfer in Porous Media > Isothermal Contours (ht).
4
Click the Add Result Template button in the window toolbar.
5
In the Results toolbar, click  Result Templates to close the Result Templates window.
Results
Multislice 1
1
In the Model Builder window, expand the Temperature, Multislice (ht) node, then click Multislice 1.
2
In the Settings window for Multislice, locate the Multiplane Data section.
3
Find the X-planes subsection. From the Entry method list, choose Coordinates.
4
In the Coordinates text field, type lx.
5
Find the Y-planes subsection. From the Entry method list, choose Coordinates.
6
In the Coordinates text field, type ly.
7
Find the Z-planes subsection. In the Planes text field, type 5.
8
Locate the Expression section. From the Unit list, choose °C.
9
Click to expand the Range section. Select the Manual color range checkbox.
10
In the Maximum text field, type 14.6.
Temperature, Multislice (ht)
1
In the Model Builder window, click Temperature, Multislice (ht).
2
In the Settings window for 3D Plot Group, locate the Color Legend section.
3
Select the Show maximum and minimum values checkbox.
Arrow Volume 1
1
Right-click Temperature, Multislice (ht) and choose Arrow Volume.
2
In the Settings window for Arrow Volume, click Replace Expression in the upper-right corner of the Expression section. From the menu, choose Component 1 (comp1) > Darcy’s Law > Velocity and pressure > dl.u,dl.v,dl.w - Total Darcy velocity field.
3
Locate the Coloring and Style section. From the Color list, choose Black.
Line 1
1
Right-click Temperature, Multislice (ht) and choose Line to add the temperatures around the boreholes to the plot.
2
In the Settings window for Line, locate the Expression section.
3
In the Expression text field, type T2.
4
In the Unit field, type °C.
5
Click to expand the Inherit Style section. From the Plot list, choose Multislice 1.
6
Locate the Expression section. In the Expression text field, type T.
7
Locate the Coloring and Style section. From the Line type list, choose Tube.
Selection 1
1
Right-click Line 1 and choose Selection.
2
In the Settings window for Selection, locate the Selection section.
3
From the Selection list, choose Array 1.
4
In the Temperature, Multislice (ht) toolbar, click  Plot.
5
Click the  Go to Default View button in the Graphics toolbar.
6
Click the  Zoom Extents button in the Graphics toolbar.
Isothermal Contours (ht)
Here, specific levels are plotted to make it easier to compare the plot with the model version using the Pipe Flow module.
Isosurface 1
1
In the Model Builder window, expand the Isothermal Contours (ht) node, then click Isosurface 1.
2
In the Settings window for Isosurface, locate the Levels section.
3
From the Entry method list, choose Levels.
4
In the Levels text field, type 10 11 12 13 14 15.
5
In the Isothermal Contours (ht) toolbar, click  Plot.
6
Locate the Expression section. From the Unit list, choose °C.
7
In the Isothermal Contours (ht) toolbar, click  Plot.
Temperature Profile
Create Figure 4 by following the steps below.
1
In the Results toolbar, click  1D Plot Group.
2
In the Settings window for 1D Plot Group, locate the Data section.
3
From the Time selection list, choose From list.
4
In the Times (a) list box, select 1.
5
In the Label text field, type Temperature Profile.
Line Graph 1
1
Right-click Temperature Profile and choose Line Graph.
2
Click the  Go to XY View button in the Graphics toolbar.
3
Click the  Select Box button in the Graphics toolbar and span a rectangle around the point belonging to the most left borehole of those marked in Figure 1. This automatically marks all edges within and below the selection box. Now the edges 28-31 should be selected.
4
In the Settings window for Line Graph, locate the y-Axis Data section.
5
In the Expression text field, type z.
6
Locate the x-Axis Data section. From the Parameter list, choose Expression.
7
In the Expression text field, type T.
8
From the Unit list, choose °C.
9
Click the  Zoom Extents button in the Graphics toolbar.
10
Click to expand the Legends section. Select the Show legends checkbox.
11
From the Legends list, choose Manual.
12
Line Graph 2
1
Right-click Line Graph 1 and choose Duplicate.
2
In the Settings window for Line Graph, locate the Selection section.
3
Click  Clear Selection.
4
Select Edges 40–43 only. You could again use the Go to XY View and the Select Box buttons from the graphics toolbar to select the edges corresponding to the middle borehole.
5
Locate the Legends section. In the table, enter the following settings:
Line Graph 3
1
Right-click Line Graph 2 and choose Duplicate.
2
In the Settings window for Line Graph, locate the Selection section.
3
Click  Clear Selection.
4
5
Locate the Legends section. In the table, enter the following settings:
6
Click the  Go to Default View button in the Graphics toolbar.
7
In the Temperature Profile toolbar, click  Plot.
Because the ambient temperature varies strongly with time and is much higher than in the subsurface layers, define axis limits to get a better view of the region of interest.
Temperature Profile
1
In the Model Builder window, click Temperature Profile.
2
In the Settings window for 1D Plot Group, locate the Axis section.
3
Select the Manual axis limits checkbox.
4
In the x minimum text field, type 5.
5
In the x maximum text field, type 14.
6
In the y minimum text field, type -150.5.
7
In the y maximum text field, type -5.
8
In the Temperature Profile toolbar, click  Plot and compare with Figure 4.