COMSOL 6.4 - Geothermal Heating from a Pond Loop

Ponds and lakes can serve as thermal reservoirs in geothermal heating applications. In this example, fluid circulates underwater through polyethylene piping in a closed system. The pipes are coiled in a slinky shape and mounted onto sleds. The Nonisothermal Pipe Flow interface sets up and solves the equations for the temperature and fluid flow in the pipe system, where the geometry is represented by lines in 3D.

Figure 1: A sled carrying pipe coils shown before the system is submerged.

Model Definition

Geometry

High density polyethylene pipe (20 mm diameter) is rolled into sixteen coils. Groups of eight coils are mounted on two sleds. Each coil has a radius of 1 m and a length of approximately 75 m. The coil groups are connected to feed and return piping with a diameter of 50 mm (see Figure 2). The coil groups are 2.4 m in height an sit at the bottom of a pond that is 6 m deep. The total length of the piping is 1446 m.

Figure 2: Polyethylene pipe system. Elevation above the pond bottom is indicated. Feed and return piping (gray) is 50 mm in diameter while coils (black) are 20 mm in diameter. The pipes are insulated above the pond surface.

The heat exchange between pond water and pipe fluid depends on the temperature difference between the two. A slow current in the pond makes the heat transfer more effective than water at rest. The pond is warmer closer to the surface, as shown by the temperature data in Table 1 below.

Table 1: Pond Temperature.
Elevation (m)	Temperature (K)
0	284
2	288
4	291
6	293

It is easy to set up a function in the software with linear interpolation between points so that the varying pond temperature can be taken into account in the simulation.

Flow Equations

The continuity and momentum equations below describe the stationary flow inside the pipe system:

(1)

Above, A (SI unit: m2) is the cross section area of the pipe, ρ (SI unit: kg/m3) is the density, u (SI unit: m/s) is the fluid velocity, and p (SI unit: N/m2) is the pressure and F (SI unit: N/m3) is a volume force, like gravity.

Gravity can be included explicitly in the model, but since the variation in density is negligible, and the model is not pressure driven, the only effect of including gravity is a change in the total pressure level. It is therefore common modeling practice to exclude gravity by setting F=0 and interpret the pressure variable as the reduced pressure pr = p − ρg(z0 − z), where z0 is the datum level of the free liquid surface. This reduces the model complexity and yields the same results.

The right-hand side of Equation 1 describes the pressure drop due to internal viscous shear. The term contains the Darcy friction factor, fD, which is a function of the Reynolds number and the surface roughness divided by the hydraulic pipe diameter, e/dh. The Nonisothermal Pipe Flow interface provides a library of built-in expressions for the Darcy friction factor, fD.

This example uses the Churchill relation (Ref. 1) that is valid for laminar flow, turbulent flow, and the transitional region in between. The Churchill equation is also a function of the fluid properties, through the Reynolds number:

The physical properties of water as function of temperature are directly available from the software’s built-in material library.

Heat Transfer Equations

The energy equation for the pipeline flow is

(2)

where Cp (SI unit: J/(kg·K)) is the heat capacity at constant pressure, T is the temperature (SI unit: K), and k (SI unit: W/(m·K)) is the thermal conductivity. The second term on the right-hand side of Equation 2 corresponds to friction heat dissipated due viscous shear. Qwall (SI unit: W/m) is a source/sink term due to heat exchange with the surroundings through the pipe wall:

(3)

where Z (m) is the wetted perimeter of the pipe, h (W/(m2·K)) an overall heat transfer coefficient and Text (K) the external temperature outside of the pipe.

The overall heat transfer coefficient includes contribution from internal film resistance, wall resistance, and external film resistance. For a circular pipe, under assumption that the heat transfer through the wall is quasi static and that the temperature is equal around the circumference of the pipe, an effective hZ in Equation 3 is given by

where rn is the outer radius of wall n, hint, and hext are the film heat transfer coefficients on the inside and outside of the tube, and kn is the thermal conductivity of wall n.

The film resistance inside the pipe is given by