The Dynamic 2D Temperature Integration

The Energy Balance

The temperature distribution in the well at a particular time is represented by one central temperature in each grid cell. Given the temperatures at a time, the objective of this computation is to compute the temperature of the next timestep.

The computation is based on applying the first law of thermodynamics (energy conservation) to each grid cell volume. The grid cell volume is the volume obtained when the grid cell rectangle in the $r-z$ plane is rotated around the well center. Similarly, the four boundaries of this ring like structure, are the areas we get when the four edges of the grid cell rectangle are rotated around the well center. With this in mind, the energy balance for a grid cell volume can be written as:

$\underbrace{\frac{\text{d}}{\,\text{d}t}\iiint u \left(T,p \right) \rho \left(T,p \right)\,\text{d}V}_\textrm{Accumulation term } = \underbrace{\iint {q}_{n} \,\text{d}A}_\textrm{Conductive term } + \underbrace{\iint{e}_{n}\,\text{d}A}_\textrm{Convective term} + \underbrace{\iiint{s}_{n}\,\text{d}V}_\textrm{Source term}$

or

where

• $\rho$ is the mass density

• $u$ is the specific internal energy (per mass)

• $q$ is the conductive energy flow

• $e$ is the convective energy flow

• $s$ is the volumetric energy source

For all grid cells except the ones in the flow path(s), there is no mass flow over the grid cell boundaries, and the last term in the equation is zero. These are all the grid cells containing a solid material (steel, cement etc.) and the fluid filled grid cells in the close annuli. The thermally induced free convection within these grid cells is treated as a modification of the thermal conductivity of the grid cell and not as a convective energy transport.

Energy is transported in the axial direction through convection in the flowpath(s) and are then transported in the radial and axial directions through conduction. Let us look at the three terms in the equation separately.

The Accumulation Term

The accumulation term can be written in terms of enthalpy, $h$, as:

$\frac {\text{d}}{\text{d}t} \iiint u \left(T,p \right) \rho \left (T,p \right )\,\text{d}V=\frac {\text{d}}{\text{d}t} \iiint \left(h \left(T,p \right) - \frac{p}{\rho} \right) \rho \left(T, p \right ) \,\text{d}V$

$T$ and $p$ are supposed to be constant within the grid cell and equal to the midpoint values. Furthermore, the density variation due to the temperature is ignored. A steady state pressure $p$ is also assumed. Then this accumulation term can be written as:

$\frac {\text{d}h}{\text{d}t} \rho V - \frac {\text{d}p}{\text{d}t} V = \frac{\text{d}h}{\text{d}T} \rho V \frac{\text{d}T}{\text{d}t} = {c}_{p} \rho V \frac{\text{d}T}{\text{d}t}.$

The extension to 3 phase flow is straightforward. The expression is a sum over the phases, where the volume fractions $\alpha_{i}$are computed in the 1D flow computation:

$\left(\sum_{i} {{\alpha}_{i}{c}_{p_i}{\rho}_{i}}\right) V \frac{\text{d}T}{\text{d}t}$

This term is inserted into the energy balance for each grid cell.

The Conductive Term

The heat conduction is based on Fourier's law for heat conduction: $\vec{q} = -k\nabla T$

where the conductive heat flow is proportional to the temperature gradient. The proportionality factor, $k$ , is the thermal conductivity for the material. This low is then integrated over a grid cell (Fig. 1). The integrated conductive heat flows passing the four boundaries, $Q_u$, $Q_b$, $Q_l$ and $Q_r$ can be expressed as:

$Q_r = C_r\left(T_r-T_c\right)$

$Q_l = C_l\left(T_l-T_c\right)$

$Q_u = C_u\left(T_u-T_c\right)$

$Q_b = C_b\left(T_b-T_c\right)$

Where $C_r$, $C_l$, $C_u$, $C_b$ are the conductances between the grid cells in the different directions, $T_r$, $T_l$, $T_u$, $T_b$ are the temperatures of the neighboring cells and $T_c$ is the temperature of the central cell.

Fig 1 The central grid cell has a midpoint temperature Tc. This grid cell has four neighboring cells (upper, bottom, left and right) with temperatures Tu, Tb, Tl and Tr respectively.

When doing this, it must be taken into account that in general the thermal conductivity $k$ and the temperature gradient are discontinuous over the grid cell boundary, while the heat flux is always continuous over the boundary. The estimation of the conductances $C_u$, $C_b$, $C_l$ and $C_r$ are done in three steps.

Step 1 - Computation of Internal Conductance in the Grid Cell

The intern effective conductance, between the center of the cell and the boundaries, are computed by applying a uniform heat flow field between the center of the grid cell and the boundaries (see figure below).

The conductances in the $z$direction are given by the thermal conductivity $k$, the area $A$ where the conductive flow goes through and the distance between the central point and the boundary of the cell:

${C}_{u}^{i} = {C}_{b}^{i} = \frac{k A} {dz/2}$

The expression for the horizontal conductances $C_r$ in a cylinder geometry is:

${C}_{r}^{i} = k \frac{2 \pi dz}{\ln \left(\frac{{{r}_{2}}}{{r}_{1}} \right)}$

where $r_1$ is the average radius of the cell and $r_2$ is the radius of the external boundary of the cell. A similar expression is obtained for $C_l^{i}$.

Step 2 - Modification of the conductance terms due to free and forced convection

In the flow path(s) fluids, either one phase or three phases, are transported. In the bulk of the flow there is normally strong turbulent mixing. Close to the wall there is a laminar boundary layer. The thickness of this boundary layer defines the heat transfer number between the bulk flow and the wall. We modify the horizontal conductance in the bulk flow by multiplying with a large number to simulate the turbulent mixing. The horizontal conductance connecting the bulk flow with the wall is modified in order to effectively give the correct heat transfer number.

Similarly, in the annuli filled with stationary fluid, there is thermally induced free convection due to the temperature gradients. This free convection gives added heat transfer over the annulus which effectively is seen as an increased thermal conductivity. This increase is expressed as a Nusselt number. We correct all the horizontal conductance over the annuli with the Nusselt number which is included in Appendix.

Step 3 - Computation of conductances based on the modified internal conductances

The figure 3 shows two internal conductances sharing the same boundary. By requiring continuity in the heat flow over the boundary, it is possible to eliminate the boundary temperature and to compute the effective conductance between the two central temperatures as:

$C_r = \frac{C_r^i C_l^{i+1}}{{C_r^i+C_l^{i+1}}}$

and the heat flow from the right to the left grid cell can then be written as:

$Q_r = C_r\left(T_r-T_c\right)$

This term is added to the energy balance for the left cell and subtracted from the energy balance for the right cell. In this way all the contributions from the heat conduction flows are included in the balances.

Fig 3 Computation of the effective conductance between the two central temperatures

The Convective Term

Let us consider $F$ the mass flow between a downstream grid cell and an upstream grid cell and look at the energy transfer during a short time interval $\delta_t$. The total internal energy transferred is:

$\Delta U = F u \, \delta_t$

where $u$ is the specific energy in the upstream grid cell. Furthermore, the mass flow displaces a volume

$\Delta V = \frac{F \, \delta_t}{\rho}$

Adding the internal energy transferred and the $p\,\text{d} V$ work done by the upstream volume on the downstream volume, it is obtained:

$\Delta U + p \Delta V = F u \, \delta_t + p \frac{F \, \delta_t}{\rho} = F \, \delta_t \left(u + \frac{p}{\rho} \right) = F \, \delta_t h$

In other words, by convicting the enthalpy instead of the internal energy, the $p\,\text{d}V$ work is automatically included.

Two additional energy flows must be included in the convection: the potential energy and the kinetic energy. If neglecting these energies, the equation will say that all the net $p\,\text{d}V$ work being done on a grid cell will be dissipated into internal energy. We know that this is not true. A major part of this work is used to lift the fluid out of the well, in other words to increase the potential energy of the fluid. Part of the potential energy is also converted into kinetic energy. Only the work associated with the frictional pressure drop is dissipated into heat. To account for this fact, we have to include the potential energy in the energy flow between the grid cell volumes. The convective energy flow over the boundary can be written as:

$F\, \left(h \left(T,p \right) - g d +0.5 v^2 \right)$

As $d$ is the True Vertical Depth, the equation is written with a minus sign. $T$ and $p$ should really be taken at the boundary. $v$ is the fluid velocity. For reasons of numerical stability, an upstream approximation is used, with the upstream values of $T$ and $p$:

$F\, \left(h \left({T}_{\text{up}} , {p}_{\text{up}} \right) - g d +0.5 v^2 \right)$

This expression is nonlinear in $T$ and is linearized before added to the balance equations. Ideally, it should be linearized around $T$, but since we do not know $T$ yet, the next best thing to do is to linearize it around the old value of $T$ which is $T_0$. Doing so one gets:

$F\left(h \left({T}_{\text{up}} , {p}_{\text{up}} \right) - g d +0.5 v^2 \right)= F \left(h \left ({T}_{\text{up}}^{0} , {p}_{\text{up}} \right) + \frac{\partial {h}}{\partial{T}} \left({T}_{\text{up}} - {T}_{\text{up}}^{0} \right) -g d +0.5 v^2 \right)$

This expression is linear in $T_{\text{up}}$, and can be added as a term in the downstream energy balance and subtracted from the upstream energy balance.

For three phase flow, the composition of the flow is constant in steady state. Therefore by using thermodynamics, the mixture properties of the flow can be computed ($h^{\text{mix}}$, $c_p^{\text{mix}}$ etc.).

In most applications the potential energy term can be neglected. However, here it is very important due to the huge height differences. In fact, the mechanical work needed to lift the fluid out of the well is typically of the order of the energy required to heat the fluid 10C. Failure to include the potential energy in the energy flows would therefore lead to a large overestimation of the temperature in the upper part of a production well.

The formulation used in this program is quite general. The phase transition heat is automatically included. Since it is based on an enthalpy flow formulation, and the variation of enthalpy with pressure is taken into account, the Joule-Thomson heating/cooling is also contained in the equations.

The Source Term

A source term is taken into account for the cementation (heat of hydration, produced during cement curing) and the drilling operations. A large part of the energy from the drilling motor is transformed into heat, generated by friction. Only a small part is actually used to cut the rock. For the other simulated operations, there is no source term.