Flow in geothermal well
Computational determination of temperature and pressure profiles in a geothermal well
Introduction
The subject of interest in the GeoModel project is the use of low-temperature geothermal resources. For this reason, it is assumed that the fluid flow in the geothermal well is single-phase, specifically the flow of the fluid in its liquid state.
The rock formation is treated as a continuous, isotropic medium without internal heat sources and with constant thermophysical properties. It is assumed that there is no fluid flow within the rock formation – the only process that is considered is heat exchange, which occurs under unsteady conditions. Temperature changes over time in the rock formation adjacent to the well are described using explicit formulation of the Finite Difference Method (FDM).
Description of the geometry of the modelled area
The assumptions made allow to formulate the process in cylindrical coordinates (Figure 1).
Figure 1. Proximity zone of the well in cylindrical coordinates, discretization scheme of the rock formation and location of the well in space
Based on the assumptions, the temperature change over time will be described by the equation:
Equation 1
It is easy to notice that the chosen discretization scheme allows for accounting for different types of materials in each element of the space. This enables for flexible specification of the well construction, including the diameters and materials of the casing pipes, cement, and the rock formation itself. It also allows for the representation of the properties of the rock formation, which may vary with depth.
Heat exchange
Heat exchange in rock formation
According to the adopted assumptions, the description of unsteady heat exchange in the rock formation can be performed using the Finite Difference Method after discretizing the rock formation in space (Figure 1). By using symmetric finite difference quotients of the appropriate order to describe the spatial derivatives, the differential form of equation 1 can be written as follows (derivations and more information can be found in the literature (Wiśniewski and Wiśniewski, 1997):
Equation 2
The unknown in the above equation is ti,j,k+1 – temperature at node i,j in the next time step k+1:
Equation 3
Initial conditions
Establishing initial conditions is one of the requirements for the uniqueness of the solution of the Fourier equation, formulated in cylindrical coordinates by Equation 1. Initial conditions may refer to scenarios where the well begins to be exploited after a long period of inactivity. If the well has been inactive long enough for the temperature field in the rock formation to return to its natural state, the initial conditions may reflect the temperature profiling curve of the well (temperature distribution with depth) – if it is known. If the thermal profile is not known, the initial conditions can be defined according to the regional geothermal gradient. In both cases, the temperature in the computational grid nodes depends only on the depth.
Another example where it is necessary to define initial conditions is the desire to restore a disturbed temperature field in the rock formation, caused by, for example, previous exploitation of the well. In this case, applying initial conditions requires assigning temperatures to all nodes in the modelled area. It is most convenient to do this by importing an input file in the appropriate data format. This will be facilitated if the program that determines the temperature distribution in the modelled space saves the temperature distribution as one of the output files
Boundary conditions
Boundary conditions are the next requirement for the uniqueness of the solution to the unsteady heat transfer equation in the modelled area. Dirichlet boundary conditions, also known as first-type boundary conditions, are proposed. Their use provides significant stability to the numerical solution, and their implementation in the analyzed process involves setting temperatures at selected nodes of the model (Figure 2). It should be noted that nodes with the first-type boundary conditions should be sufficiently distant from the well so that thermal disturbances caused by fluid flow do not reach the nodes with fixed temperatures. In practice, for the analyzed process, this means setting conditions at the boundary of the cylindrical domain at distances of several tens of meters (~50 m), and for the bottom, distances of several to 20 meters. The suggested model extents are estimated based on (Pająk and Bujakowski, 2000).
Figure 2. Suggested scheme for determining boundary conditions
In the case of the near-surface zone, an additional layer of elements has been added where a first-kind boundary condition is set. This layer replaces the atmospheric air layer and reflects the heat exchange effects between the surface soil layer and the atmospheric air. The density and specific heat of the material replacing the air layer correspond to typical air parameters. However, the thermal conductivity coefficient is a substitutive value (significantly higher than the thermal conductivity coefficient of the stationary air layer), which has been determined to account for the convective heat exchange conditions occurring in the air layer near the Earth’s surface. Assuming that the average annual wind speed in Poland is approximately 3.24 m/s (Dygulska and Perlańska, 2015), the convective heat transfer coefficient can be estimated as (Recknagel, Sprenger, Hönmann, and Schramek, 1994).
Equation 4
If we assume that the heat exchange between the atmospheric air and the ground causes a temperature change in the air layer up to Δz = 2030 cm, and we introduce an effective heat conductivity coefficient λair‘ with such a value that the temperature difference between the ground surface and the air Δtsa, remains unchanged:
Equation 5
Hence, the value of the effective heat conductivity coefficient of air λair’, adopted in Figure 2, takes into account the convective heat exchange with the environment.
Heat exchange between the fluid and the inner wall of the well
The temperature distribution of the fluid flowing through the well can be described using the zone balance method. Assuming that mass balance is maintained in the well (the mass flow rate m ̇ remains constant throughout the exchanger), the temperature in the analyzed heat exchange zone n (Figure 1) changes linearly from the inlet temperature t(n)j−1 to the outlet temperature t(n)j. The thermal power exchanged between the fluid flowing through the well and the rock formation causes a change in the fluid temperature, according to the equation:
Equation 6
The above equation can be solved iteratively by determining the temperature of the fluid at the outlet of zone n. This allows for an accurate determination of the specific heat capacity of the fluid in the zone, which significantly impacts the effect of its heating or cooling.
The heat transfer occurring between the fluid and the inner side of the casing pipes is the convective heat transfer and is described by the equation:
Equation 7
The value of the temperature twsf(n)k on the inner surface of the casing pipe in zone n, at height j, at the current time (time step) k, can be determined using the equation:
Equation 8
The convection heat transfer coefficient can be determined using the Nusselt number Nu:
Equation 9
In the case of laminar flow in a pipe (Reynolds number < 2000), the convective heat transfer coefficient can be calculated using the equation 8 (Wiśniewski and Wiśniewski, 1997):
Equation 10
In the case of turbulent flows (Re > 2000), for water and other low-viscosity fluids (1 < Pr < 20), the following empirical relationship can be used (Wiśniewski and Wiśniewski, 1997):
Equation 11
Determination of the properties of the geothermal brine
In the case of wells extracting geothermal fluid, it should be considered that the fluid will not be pure water but rather a brine. From the perspective of modelling fluid parameters in a geothermal well, mineralization affects the specific heat, density and viscosity of the fluid. Parameters of the fluid, accounting for mineralization, can be determined based on literature. Equations for brines properties are proposed on the project’s website (IGSMiE PAN, 2024).
Flow resistance
The pressure distribution within the well can be determined by calculating the flow resistance of the geothermal fluid through each zone n of the well ∆p(n) (Recknagel, Sprenger, Hönmann, and Schramek, 1994):
Equation 12
The total pressure of the water at the end of the n-th zone of the well, at time step k, is determined by the following relationship:
Equation 13
The static pressure of the fluid at the outlet of zone n, at time step k denoted as ps(n)j,k – depends on fluid parameters such as viscosity, density, specific heat. It is determined by subtracting the dynamic pressure, the hydrostatic pressure drop ∆ps(n)k and the flow resistance ∆p(n)k from the total pressure at the beginning of the zone pc(n)j-1,k:
Equation 14
Literature:
Dygulska, A. & Perlańska, E. (2015). Mapa wietrzności Polski. Słupsk: Akademickie Centrum Czystej Energii (https://kierunkizamawiane.upsl.edu.pl/pliki/czystaenergia/raport2_II.pdf).
IGSMiE PAN. (2024, March 15). Geomodel. Retrieved from Parametry termofizyczne wody zasolonej: https://geomodel.pl/en/thermophysical-parameters-of-brines/
Pająk, L. & Bujakowski, W. (2000). Efektywność wykorzystania głębokich odwiertów wiertniczych jako wymienników ciepła wykorzystujących energię geotermiczną. Kraków: IGSMiE PAN.
Recknagel, H., Sprenger, E., Honmann, W., & Schramek, E. R. (1994). Poradnik Ogrzewanie i Klimatyzacja z uwzględnieniem chłodnictwa i zaopatrzenia w ciepłą wodę. Gdańsk: EWFE – Edition 1.
Wiśniwski, S. & Wiśniewski, T. S. (1997). Wymiana ciepła. Warszawa: Wydawnictwa Naukowe i Techniczne.
The project is financed under the Fund for Bilateral Relations through the European Economic Area Financial Mechanism (EEA FM) and the Norwegian Financial Mechanism (NFM) 2014-2021, programme “Environment, Energy and Climate Change”.
Project Operator: Ministry of Climate and Environment